© The Authors 2024

1 Introduction

The study of metal-poor massive O-type stars has received renewed interest in recent years. They dominate the rest-frame ultraviolet (UV) spectroscopic appearance of high-redshift (z) galaxies (Rix et al. 2004) and are the source of the ionizing radiation responsible for their rest-frame optical and UV nebular properties (Steidel et al. 2014; Lecroq et al. 2024). Early release observations with JWST (Pontoppidan et al. 2022) have already identified metal-poor, star-forming galaxies at z ≥ 7 (e.g., Schaerer et al. 2022; Arellano-Córdova et al. 2022; Arrabal Haro et al. 2023; Curti et al. 2023; Robertson et al. 2023; Trussler et al. 2023), highlighting the vital role of massive stars in the first gigayear and the demand for accurate knowledge of massive stars in metal-poor environments. in addition, metal-poor massive binaries are considered to be progenitors of black hole and neutron star mergers (e.g., Neijssel et al. 2019; Boco et al. 2021; Stevance et al. 2023), which are increasingly frequently detected by gravitational wave observatories (e.g., Abbott et al. 2021).

High-quality optical spectroscopy of O stars in the Magellanic Clouds – our closest metal-poor star-forming galaxies – was scarce until the advent of multiobject and integral field spectrographs on large ground-based telescopes greatly improved the samples (e.g., Evans et al. 2004a; 2006, 2011; Castro et al. 2018; Ramachandran et al. 2019). In the near future, the 1001 MC survey performed with 4MOST will provide a further order-of-magnitude increase in sample size (Cioni et al. 2019). In contrast, far-ultraviolet (FUV) spectroscopy of Magellanic Cloud OB stars – directly sampling P Cygni wind diagnostic lines – remains exceptionally scarce, and has often been limited to low-spectral-resolution observations (e.g., Walborn et al. 1995, 2000; Crowther et al. 2016; Rickard et al. 2022).

A significant new Hubble Space Telescope (HST) initiative, the Ultraviolet Legacy Library of Young Stars as Essential Standards (ULLYSES, Roman-Duval et al. 2020), seeks to address this deficiency through the acquisition of high-quality UV spectroscopy of several hundred OB stars in the Magellanic Clouds, each of which has also been observed in the optical range with the X-shooter spectrograph at the Very Large Telescope (VLT) in the framework of the XShooting ULLYSES (“XShootU”) initiative (Vink et al. 2023). Historically, detailed studies of individual O stars in the Magellanic Clouds involved application of planeparallel model atmospheres not assuming local thermodynamic equilibrium (“non-LTE”) to optical spectroscopy (Lanz et al. 1996; Heap et al. 2006), with stellar winds handled separately (Puls et al. 1996). Spherical, non-LTE model atmosphere codes incorporating the effects of metal line blanketing were subsequently developed, namely CMFGEN (Hillier & Miller 1998), WM-BASIC (Pauldrach et al. 2001), PoWR (Gräfener et al. 2002), and FASTWIND (Puls et al. 2005). Of these, all are capable of analyzing UV and optical spectroscopy, with the exception of WM-BASIC, whose focus is on UV spectroscopic studies (e.g., Garcia & Bianchi 2004).

These sophisticated model atmosphere codes have been applied to optical spectroscopic samples of OB stars in the Magellanic Clouds; see, for example, Bestenlehner et al. (2014), Ramachandran et al. (2018b), and Massey et al. (2009) and Rivero González et al. (2012) using CMFGEN, PoWR, and FASTWIND, respectively. Massey et al. (2013) made a rare comparison of two particular codes, finding broad agreement in the temperatures of O stars using CMFGEN and FASTWIND, although systematically lower gravities were obtained with the latter. Combined UV and optical spectroscopic studies of OB stars in the Clouds were rare until recently (Ramachandran et al. 2018a; Bouret et al. 2021; Hawcroft et al. 2021; Brands et al. 2022; Rickard et al. 2022), albeit with notable exceptions (Crowther et al. 2002; Hillier et al. 2003; Evans et al. 2004b; Bouret et al. 2013).

The era of ULLYSES/XShootU permits combined UV and optical studies of large samples of OB stars in the Magellanic Clouds, but it is critical to quantify any systematic differences between the model atmosphere codes and the various techniques employed by individual groups. This is the focus of the present study, where we analyze ULLYSES/XShootU spectroscopy of representative O stars in the Magellanic Clouds – showing prominent stellar winds – with different methods and provide a detailed comparison of the results. The paper is structured as follows: In Sect. 2, we present a summary of the UV and optical spectroscopic datasets. Section 3 outlines the various analysis techniques, with the comparison of results being presented in Sect. 4. A discussion of the implications of the derived stellar and wind parameters follows in Sect. 5 before we draw conclusions in Sect. 6. In the Appendix, we provide detailed information about the different codes and methods as well as large parameter and atomic data comparison tables.

2 Sample

To compare different analysis techniques, we selected three stars from the ULLYSES database that sample spectral types from early to late O-type, and luminosity classes from dwarfs to super-giants, with no a priori indication for binarity: We selected one O5 dwarf in the SMC, namely AzV 377, which has a fine-classification as O5 V((f)) (Evans et al. 2004a) following the additional “((f))” notation from Walborn (1971). This means that the He II 4686Å line is in absorption while the N III line complex at 4634–4640–4642 Å is seen in emission. The two other sample stars are located in the LMC and consist of the O9 supergiant Sk-66° 171 (classified as O9Ia by Fitzpatrick 1988) and the O7 star Sk -69° 50, which does not formally have a luminosity class and belongs – with its O7(n)(f)p fine classification (Walborn et al. 2010) – to the group of Ofnp stars (Walborn 1973). This group is marked as peculiar (“p”) and is characterized by broadened absorption lines (“n”) as well as the above-mentioned “f”-character, albeit with the involved N III lines portraying stronger emission in the case of Sk -69° 50 compared to AzV 377, hence the fine classification with single parenthesis compared to the double parenthesis designation for AzV 377. The SMC star (AzV 377) is the only one to have been analyzed previously with quantitative spectroscopy (Massey et al. 2004).

Optical and near-infrared (NIR) photometry is gathered from the ULLYSES project for each star, and is summarized in Table 1. As introduced in the first paper of the XShootU series (Vink et al. 2023), we adopt d_SMC = 62.44 kpc, corresponding to DM(SMC) = 18.98 mag (Graczyk et al. 2020), and d_LMC = 49.59 kpc, corresponding to DM(LMC) = 18.48 mag (Pietrzyński et al. 2019).

2.1 ULLYSES ultraviolet spectroscopy

The three stars have been observed with the Far Ultraviolet Spectroscopic Explorer (FUSE, Moos et al. 2000), providing spectroscopic coverage of λλ905–1187 Å (R ~ 20 000; for an OB atlas see Walborn et al. 2002). AzV 377 has been observed with HST in the FUV with the Cosmic Origins Spectrograph (COS, Green et al. 2012) using the G130M/1291 (λλ1132–1433 Å, R ~ 14 000) and G160M/1611 (λλ1419–1790 Å, R ~ 14 000) gratings, while the Space Telescope Imaging Spectrograph (STIS, Kimble et al. 1998) was used for Sk -69° 50 (added later to the ULLYSES dataset) and Sk-66° 171 with the E140M grating (λλ1143–1710 Å, R ~ 46000). For the latter star, the spectral coverage extends into the near-UV (NUV) due to an additional observation with STIS using the E230M/1978 grating (λλ1607–2366 Å, R ~ 30 000). Only the observations of Sk-66° 171, performed on January 28, 2022 (GO/DD 16365, PI Roman-Duval), were obtained within the DDT provided for the ULLYSES project, while the UV observations for AzV 377 were part of GO 15837 (PI Oskinova), performed on June 25, 2020, and the observations of Sk -69° 50 date back to October 11, 2011 and were part of GO 12218 (PI Massa).

2.2 XShootU optical spectroscopy

Optical, normalized X-shooter (Vernet et al. 2011) spectroscopy of each star from VLT/X-shooter was reduced and processed according to Sana et al. (2024, eDR1). In this work, we use the reduced data for the UBV (λλ3100–5600 Å, R ~ 6700) and VIS (λλ5600–10 240 Å, R ~ 11 400) arms. The data reduction for the NIR arm requires additional work and is therefore not yet available. For the purpose of our O star analysis, the broad wavelength coverage from the combined UV and optical spectra contains a sufficient number of spectral lines from different elements and ionization stages in order to avoid any deficiencies due to the absence of NIR data.

3 Analysis methods

The three targets in the present work are limited to the regime of O stars. All selected stars have noticeable winds that leave an imprint (i.e., diagnostic) in the spectrum and therefore mark ideal targets for our method comparison. Thus, we only use the expanding atmosphere codes CMFGEN, PoWR, and FAST-WIND. For O dwarfs with weak winds, plane-parallel model atmosphere codes, such as TLUSTY (Lanz & Hubeny 2003, Lanz & Hubeny 2007), would be suitable as well. Investigations of ULLYSES B-type stars will be presented in subsequent papers.

3.1 Common aspects

The three atmosphere codes utilized in this work – CMFGEN, FASTWIND, and PoWR – are 1D codes assuming spherical symmetry and a stationary outflow. Targeting hot stars, they account for an expanding, non-LTE environment by iteratively solving the equations of statistical equilibrium for a large set of levels of individual elements and ions together with the solution of the radiative transfer. For POWR and CMFGEN, the radiative transfer is completely solved in the co-moving frame (CMF). In the FASTWIND version (Sundqvist & Puls 2018) applied in this work, only a few “explicit” elements (cf. Sect. 3.2), plus the most important lines from the other elements, are treated in the CMF. For all other lines, a Sobolev approach with a pseudo-continuum irradiation – accounting for the combined line-opacity/emissivity from the metallic background – is used. Both approximations enable comparatively short computation times. After discussing common aspects and tools in this subsection, the following subsections introduce the different model atmosphere codes and their specific application methods. In the Appendix, we provide a more in-depth review of the physical treatments for all expanding atmosphere codes employed in this work (Appendix A). Detailed method descriptions sorted by the different aspects necessary (including aspects such as the determination of the projected rotational velocity or the bolo-metric luminosity) for quantitative spectral analysis are given in Appendix B. In the following, each individual method is denoted by a letter and a number with the letter denoting the initial of the employed atmosphere code.

3.1.1 Velocity and density structure

The models computed in this work are not dynamically consistent, but assume a prescribed velocity v(r) in the form of a so-called β-law $$v(r)=v_{\infty }{\left(1-\frac{R_0}{r}\right)}^{\beta },$$(1)

where v_∞ is the velocity for r → ∞, ß is a free parameter, and R₀ is a reference radius. In methods where a pre-calculated grid of models is used, ß is often fixed to a specific value; for example, ß = 0.8 motivated by extensions of the CAK (named after Castor et al. 1975) theory (e.g., Pauldrach et al. 1986). When individual models are calculated, the value of ß can instead be determined from combining constraints from UV and optical lines that are affected by the stellar wind. In the subsonic part, the models aim to obtain a (quasi-)hydrostatic stratification, with the detailed techniques and the connection to the supersonic ß-law differing between the codes. The solution techniques vary between the different codes. A notable difference with respect to the derived values of the surface gravity log g can arise due to different assumptions in the solution of the hydrostatic equation $$\frac{\text{d}P}{ \text{d}r}=-\rho (r)\left[g(r)-g_{\text{rad}}(r)\right],$$(2)

with P denoting the pressure and ρ the matter density. In addition to the required knowledge of the radiative acceleration g_rad(r), the solution of Eq. (2) demands an equation of state. In all model codes used in this work, this is the ideal gas equation, which we can write as P(r) = ρ(r)a²(r). For the speed a, there is the opportunity to not only include the (thermal) speed of sound, but also a turbulence term, such that $$a^2(r)=\frac{k_{\text{B}}T(r)}{\mu (r)m_{\text{H}}}+\frac{1}{2}{\xi }^2(r),$$(3)

with T(r) denoting the (electron) temperature, μ(r) the mean particle mass (including electrons), k_B Boltzmann’s constant, and m_H the mass of a hydrogen atom. From the codes used in this work, only PoWR has the option to include a nonzero term for a turbulent pressure described by a velocity ξ(r) when solving the hydrostatic equation (Sander et al. 2015). The value of ξ(r) in Eq. (3) can – but does not have to – be chosen similar to the microturbulence entering the formal integral. The use of ξ > 0 in the solution of the hydrostatic equation will lead to larger values for the determined log g. The difference can be estimated via $$\Delta (\log g)=\log \left(1+\frac{{\xi }^2\mu m_{\text{H}}}{2k_{\text{B}}T}\right).$$(4)

3.1.2 Wind inhomogenieties

All atmosphere codes take wind inhomogenieties (“clumping”) into account. Most of the applied methods only make use of the so-called “microclumping” approximation, whereby it is assumed that clumping is limited to small scales and the clumps themselves are optically thin. The medium between the clumps is assumed to be void. Defining an average density via the (stationary) equation of continuity $$\langle \rho \rangle =\frac{\dot{M}}{4\pi r^{2}v(r)},$$(5)

with $\dot{M}$ denoting the mass-loss rate, one can define a “clumping factor”: $$f_{\text{cl}}=\frac{\langle {\rho }^2\rangle }{{\langle \rho \rangle }^2}.$$(6)

Inside the clumps, the over-density relative to a smooth wind can be described by a factor D = ρ_cl/ρ_smooth with ρ_cl denoting the density inside the clumps and ρ_smooth denoting the density of a smooth wind with the same mass-loss rate. The mean density can further be expressed as 〈ρ〉 = f_Vρ_cl, with f_V describing the volume filling factor of the clumps.

In practice, the different atmosphere codes employ different quantities as their free parameter: FASTWIND uses f_cl, while CMFGEN requires f_V to be given, and PoWR has D as its free parameter. Fortunately, these values can easily be converted into each other for a void interclump medium and optically thin clumps, namely $$f_{\text{cl}}\equiv D\equiv f_{\text{V}}^{-1}.$$(7)

This relation does not hold for optically thick clumps or a non-void interclump medium, which is used in one of the employed methods (F3). The more detailed clumping implementations in the different codes, including the optically thick clumping approach in FASTWIND, are described in Appendix A.2. For f_cl = 1, a smooth (“unclumped”) wind situation is recovered in all cases.

3.1.3 Abundance notations

The input and output format for abundances differ between the atmosphere codes. While for example PoWR expects either mass fractions or absolute number fractions to be given, FASTWIND requires number ratios and CMFGEN can handle a mixture of number ratios and mass fractions. Fortunately, these quantities can easily be converted as long as information on all elements with major abundances is provided. From a given set of either absolute number fractions N(i) or number ratios relative to an element such as hydrogen N(i)/N(H), the (absolute) mass fraction X_i for an arbitrary element i can be determined via $$X_i=\frac{{\mathcal{A}}_i\frac{N(i)}{N(\text{H})}}{{\displaystyle \sum _{j=1}^{n_{\text{clem }}}{\mathcal{A}}_j}\frac{N(j)}{N(\text{H})}}=\frac{{\mathcal{A}}_{i}N(\text{i})}{{\displaystyle \sum _{j=1}^{n_{\text{clem }}}{\mathcal{A}}_j}N(\text{j})},$$(8)

with 𝒜_i; denoting the atomic weight of element i. If absolute mass fractions are provided, the absolute number fractions are given by $$N(i)=\frac{X_i}{{\mathcal{A}}_i{\displaystyle \sum _{j=1}^{n_{\text{clem }}}\frac{X_j}{{\mathcal{A}}_j}}}.$$(9)

In FASTWIND and CMFGEN, the He/H number ratio $$y_{\text{He}}=\frac{N(\text{He})}{N(\text{H})}$$(10)

is an input quantity, and is commonly also simply denoted Y in the literature. We refrain from using the latter notation and instead use y_He in this work to avoid any confusion with the common (X, Y, Z) notation, which refers to the mass fractions of hydrogen, helium, and all other elements, respectively, with these latter being commonly referred to as “metals” in astrophysics. This fraction of metals Z is referred to as metallicity.

A further common astrophysical abundance notation is $$ϵ(x)=\log \frac{x}{H}+12\equiv \log \frac{N(x)}{N(\text{H})}+12,$$(11)

with the first expression referring to the typical literature standard and the second being the equivalent in our notation accounting for the different format specifications. In the extragalactic community, Z is sometimes also used as a label for [O/H] = ϵ(O) − ϵ(O)_⊙ or [Fe/H]. When gas-phase abundances are measured from nebular lines, ϵ(O) is often treated as a proxy for Z or is even termed “metallicity”. We only use the term in its original meaning, referring to the total mass fraction of all elements beyond helium.

3.1.4 Rotation, macroturbulence, and lACOB-Broad

To determine the projected rotational and macroturbulent velocities, which broaden the spectral lines, a couple of methods use the iacob-broad¹ package described in Simón-Díaz & Herrero (2014, see also references therein). iacob-broad combines the Fourier transform method (based on the presence of zeros introduced by the transform function of the rotational profile) and the “goodness of fit” method (based on the best fit of a combination of rotational and radial-tangential macroturbulent profiles) to the observed spectral lines. In their package, Simón-Díaz & Herrero (2014) make three major assumptions: the stellar surface is spherical, the rotational and macroturbulent profiles are convolved with the emergent flux profiles (and not with the intensity profiles), and other broadening mechanisms (e.g., collisional and instrumental) are comparatively small. The last approximation implies that H and He lines should be avoided if possible, as they will be broadened by the linear Stark effect². The selection of lines is decided by the user and can be adjusted for each star depending on the available spectra and the strength of the individual lines therein.

The results obtained for the projected rotational velocity v sin i and the macroturbulence v_mac can be degenerate, depending on the spectral resolution, the available lines, and – if a Fourier transform method such as iacob-broad is used – the selection of the correct zero in Fourier space (see, e.g., the discussions Simón-Díaz & Herrero 2007, 2014, for more details). Therefore, some of the methods employed in this work choose to fix v_mac, while others keep it as a free parameter. The origin of macroturbulence in massive stars and the interpretation of the derived v_mac values are a topic of active research (e.g., Aerts et al. 2009; Sundqvist et al. 2013; Grassitelli et al. 2016; Debnath et al. 2024).

3.2 FASTWIND (F methods)

Developed with the intent to provide a computationally fast non-LTE scheme (Santolaya-Rey et al. 1997), FASTWIND focuses on providing models and synthetic spectra for OBA stars with winds that are not significantly optically thick in the (optical) continuum. The initial efforts of the code are documented in Santolaya-Rey et al. (1997) with subsequent improvements and extensions described in Puls et al. (2005); Rivero González et al. (2012); Carneiro et al. (2016); Sundqvist & Puls (2018). Unlike the other codes applied in this work, FASTWIND distinguishes between line and continuum transfer as well as “explicit” and “background” elements. Only the explicit elements have a flexible, user-supplied model atom³ and employ the CMF radiative transfer for their line transitions. The radiative transfer for the background elements is mainly performed with the Sobolev (1960) approximation, but the most important transitions can be done in the CMF as well (Puls et al. 2005). There is also a recent version of FASTWIND that can treat all elements in the CMF (Puls et al. 2020), but this version has so far only been employed to perform mass-loss predictions (e.g., Björklund et al. 2021) and is not used in this work.

The input reference radius for all FASTWIND models is the radius corresponding to an effective temperature for a Rosseland optical depth of τ_Ross = 2/3. In the literature, the corresponding temperature is commonly termed T_eff, while the radius is denoted R_*. However, this designation is not unique among the different atmosphere codes. While we use the established label T_eff to refer to the effective temperature at τ_Ross = 2/3, the corresponding radius is denoted R_2/3 throughout this paper in order to avoid any confusion between the different codes that use the label R_* for different radii.

To describe the strength of the wind, the mass-loss rate $\dot{M}$, terminal velocity v_∞, and clumping factor f_cl (assuming optically thin clumping with no interclump medium) can be combined into the “wind strength parameter”: $$Q_{\text{ws}}=\frac{\dot{M}\sqrt{f_{\text{cl}}}}{M_{\odot }{\text{yr}}^{-1}}{\left(\frac{\text{km\hspace{0.17em}}{\text{s}}^{-1}}{v_{\infty }}\frac{R_{\odot }}{R_{2/3}}\right)}^{3/2},$$(12)

which is a common input parameter for FASTWIND models, in particular when calculating grids of models. The quantity Q_ws was originally defined in Puls et al. (1996), later adjusted for clumping (e.g., Puls et al. 2008), and is also known as “optical depth invariant”. If the wind is optically thin, models with the same stellar parameters and the same Q_ws yield very similar spectra, allowing a reduction of the calculation effort for model grids. For optically thick winds, the v_∞ scaling changes and instead the “transformed radius” R_t (Schmutz et al. 1989) (or the “transformed mass-loss rate” ${\dot{M}}_{\text{t}}$ discussed in Sect. 5) are better scaling quantities in this regime (see, e.g., Bestenlehner et al. 2020). In the subsequent sections, we briefly introduce the general concepts of all methods employing FASTWIND.

3.2.1 F1 – Optical/IACOB-GBAT

The F1 method makes use of a grid-based automatic tool (GBAT) called IACOB-GBAT (Simón-Díaz et al. 2011), which was developed as part of the IACOB project (Simón-Díaz & Herrero 2014) and is regularly applied there (e.g., Holgado et al. 2018, 2020). IACOB-GBAT determines the goodness of a fit within a given grid of atmosphere models via a χ² criterion applied on a list of selected, normalized lines.

Only the optical H/He spectra are used in the F1 method, with H and He being the only explicit elements. A small grid of models was calculated for the SMC star, while a much larger grid is employed for the two LMC stars. The parameter range for the LMC model grid is listed in Table 2. The χ² calculation enables us to also estimate the uncertainties of the derived stellar (and wind) parameters. No wind clumping is included in any of the F1 models (f_cl = 1), meaning that any estimates of the mass-loss rate are only upper limits.

3.2.2 F2 and F3 – Kiwi-GA: Optical and optical + UV

The F2 and F3 methods use a similar approach (see below) to derive the stellar and wind parameters plus the He and CNO abundances (by means of H, He, C, N, O, Si, P as explicit elements), but differ in the usage of the underlying data. In F2, only the optical data are taken into account, while F3 uses both optical and UV data. For the optical spectra, the normalization of Sana et al. (2024) is adopted, but the data are renormalized where the continuum clearly lies above unity. F2 and F3 make use of FASTWIND (v10.6, Sundqvist & Puls 2018) with optically thick wind clumping (macroclumping), combined with a genetic algorithm (GA) called Kiwi-GA⁴ (Brands et al. 2022). Earlier forms of this method have been used in several analyses of massive star spectra (e.g., Mokiem et al. 2005; Tramper et al. 2014; Abdul-Masih et al. 2021; Brands et al. 2022). Genetic algorithms are based on the concepts of natural selection and “survival of the fittest”. First, an initial group of model input parameters is selected randomly from a given parameter space. The “fitness” of the resulting model spectra is tested against the data – a stellar spectrum – by computing a χ² value for a selection of normalized lines, thus deciding which parameters are selected for the next generation of models: parameters of the models with the lowest χ² value have the greatest chance of being selected. With the new parameters, but also random “mutations”⁵, models of the next generation are computed and their fitness is analyzed again. This process is repeated for 40–120 generations, after which the algorithm converges to a set of best-fit parameters. The χ² values further enable the calculation of uncertainties for the best-fitting model. More details about Kiwi-GA are given in Brands et al. (2022). For details regarding the uncertainty derivation, see Brands et al. (in prep, part of the XShootU series). Requiring no model grid, but instead the calculation of new models on the fly, the GA concept for spectral fitting has so far only been combined with FASTWIND atmosphere models due to their short computing times (15–45 minutes).

Technically, the F3 analysis is not performed independently, but builds up on F2. In F2, ß is fixed to unity, and also the clumping parameters are fixed (for details, see Appendix B.4).

F3 then allows us to vary ß and the full set of wind and clumping parameters in FASTWIND, but fixes the projected rotational velocity and the y_He ratio obtained in F2.

3.2.3 F4 – LMC optical model grid

The F4 method uses a grid-based approach, but is performed with a different set of models and a different pipeline than F1, though also minimizing the χ² (see Appendix B.2). The underlying model grid has dedicated LMC abundances and is calculated with FASTWIND v10.6, with H, He, C, N, O, and Si as explicit elements. The grid explores the T_eff, log g, $\dot{M}$, and y_He parameter space plus three CNO abundance combinations. The CNO abundance sets represent LMC baseline abundances (Vink et al. 2023) plus semi and fully processed CNO composition due to the CNO-cycle according to the 60 M_⊙ evolutionary track by Brott et al. (2011). In the grid, a smooth wind is assumed (i.e., clumping factor f_cl = 1), the wind velocity field uses a fixed value of ß = 1.0, and the microturbulence velocity is fixed to 10 km s⁻¹. In total, around 120 000 stellar models were computed. Including convolutions for rotation results in a total number of about 1 100 000 synthetic spectra.

The spectroscopic analysis in F4 is solely based on the optical VLT/X-shooter data (Sana et al. 2024, eDR1) and uses the spectral lines of H, He I-II, CII-IV, NII-V, and Si II-IV in the wavelength range of λλ3800–7100 Å. Further redward wavelengths have been ignored to avoid any impact of telluric lines on our results. Wavelengths of <3800 Å were omitted to avoid spurious effects from the normalization around the Balmer jump. For reproducing the line spectra, the normalized spectra provided by Sana et al. (2024, eDR1) were used without further renormalization. The uncertainties of the best-fitting model were derived with an empirical Bayesian approach and maximum a posterior approximations utilizing de-idealized models as described in detail in Bestenlehner et al. (2024).

3.3 CMFGEN (C methods)

The stellar atmosphere code CMFGEN is a spherical, non-LTE code that was developed to model stars with strong winds (Hillier 1990; Hillier & Miller 1998) that can also be optically thick in the continuum. Its radiative transfer is performed completely in the CMF with a detailed treatment for line-blanketing using a flexible superlevel approach (Hillier & Miller 1998). While there is a time-dependent branch for the simulation of supernovae spectra (e.g., Dessart & Hillier 2010), the scheme we employ in this work assumes stationary outflows.

Unlike in FASTWIND, CMFGEN uses the notation R_* for the inner boundary in its input file. In general, R_* will correspond to τ_Ross ≫ 2/3. The temperature T_eff = T_2/3 and log g (also referring to R_2/3) are only a relevant input when iterating the density structure for the hydrostatic domain, which is necessary for objects such as the O stars studied in this work, where the spectrum is not completely formed in the wind.

3.3.1 C1 –χ² analysis

The C1 method is a grid-based approach with subsequent refining via additional model sets. The underlying grids of CMF-GEN models are computed for stars in the Magellanic clouds, which will be presented in Marcolino et al. (2024). The models of these grids are computed for scaled solar abundances (CNO and beyond). A microturbulence velocity of 10 km s⁻¹ is adopted for the line profiles in the CMF radiative transfer and in the spectrum calculation throughout the initial grids. No clumping is assumed in the initial grids (f_V = 1). Using a χ² criterion on a series of optical-only lines (cf. Appendix B.2), the best-fitting grid model is determined.

For the two LMC stars, the spectral fit is subsequently improved by calculating further models that can go beyond the grid assumptions for wind clumping, microturbulence, abundances, and so on. In this process, another χ² analysis is performed on an extended line set including CNO lines. Finally, the wind parameters were adjusted to reproduce the wind (UV and Hα) spectral features. More details on these steps are given in Appendix B and the whole methodology is presented in more depth in Martins et al. (2024, paper V in the series). Uncertainties are derived with the help of the obtained χ² values from the different models.

3.3.2 C2–Individual fit

The C2 method employs the traditional method of calculating a series of individual atmosphere models to obtain a reasonable manual (“by eye”) fit to the observed spectrum. The initial assignment of T_eff, L, and log g is based on the spectral types of the sample stars and the O-star calibration by Martins et al. (2005). No specific parameter restrictions were made apart from testing only single-ß velocity laws and assuming optically thin wind clumping.

Only the two LMC stars are analyzed in the framework of C2. Initial abundance estimations for them employ Geneva evolutionary tracks from Eggenberger et al. (2021), but further refinements were made when necessary. The uncertainties are estimated by comparing the best-fitting solution to models with varied parameters.

3.4 PoWR (P methods)

The fundamental physical approach to the non-LTE stellar atmosphere modeling is similar between PoWR and CMFGEN, while the development and implementation of these codes are completely independent. While similarly designed for hot stars with significant winds, including those with optically thick continuum, the numerical approaches differ considerably, for example regarding the implementation of iron-line blanketing (Gräfener et al. 2002) or the determination of the temperature stratification (Hamann & Gräfener 2003). PoWR provides the opportunity to calculate models completely from scratch, but the most common way to start the spectral analysis is to select a model from an older study or from a previously calculated grid (e.g., Hainich et al. 2019).

3.4.1 Specific notations and hydrostatic domain treatment

Due to its roots in the analysis of Wolf-Rayet stars, POWR does not use T_eff = T_2/3 as an input parameter, but instead uses the inner boundary radius, termed R_*. This radius is associated to a maximum Rosseland continuum optical depth, which is typically set to τ_Ross,cont = 20. As this is also the case for the present models, we denote the corresponding radius as R₂₀ to avoid any confusion from the different usages of the label R_*. The associated effective temperature – typically termed T_* in papers employing PoWR or CMFGEN – is labeled T₂₀ here. The quantities T_eff and R_2/3 are output parameters from converged models. The input surface gravity log g is also specified at R₂₀, but we instead provide the corresponding value for g(R_2/3) to ease the direct comparison.

In the hydrostatic regime, PoWR integrates the hydrostatic equation using directly the radiative force Γ_rad = a_rad/g calculated in the comoving frame to obtain the density and velocity stratification. When starting a new model, the initial integration is either performed with a mean Γ_rad or a depth-dependent description is taken from an old model.

3.4.2 P1/P2 – Individual fits with tailored models

The P1 and P2 methods consist of a series of individual model calculations to obtain a reproduction of the observed spectrum that is deemed sufficient upon visual inspection. The initial models are selected from publicly available OB grids (Hainich et al. 2019) with parameters as close as possible to the assumed stellar parameters. After constraining the rotational velocity (and macroturbulence), tailored models with adjusted stellar and wind parameters are calculated until the synthetic spectrum sufficiently reproduces the observed spectrum. The models further vary the depth-dependent optically thin clumping and the P2 models further include additional X-rays.

P1 and P2 differ slightly in their detailed assumptions and fixed inputs (e.g., in microturbulence entering the hydrostatic equation and the usage of the wind velocity law; see Appendix B for details). P1 is only applied to the SMC star AzV 377, while P2 is limited to the two LMC stars. The error margins quoted for P1 and P2 are determined by varying the individual stellar and wind parameters of the final synthetic model spectrum and include model parameters that still mimic the observed spectrum.

4 Comparison of the spectral analyses

For each of our three sample stars, the results from all methods applied to the specific object are listed in a separate table. For the O5-dwarf AzV 377, the results are listed in Table C.1. The values for the peculiar O7(n)(f)-star Sk -69° 50 are provided in Table C.2, and the results for the O9-supergiant Sk -66° 171 are given in Table C.3. A brief overview of the maximum spread in the derived fundamental parameters is shown in Table 3. The last column further indicates a “typical” spread in each parameter arising from our sample. This value is calculated as a rounded average of the values for the different targets and approximately reflects the systematic uncertainty arising from the different analysis codes and methods.

Fig. 1 Comparison between the different methods in the Balmer lines and He λ14541 region for AzV 377.

4.1 Spectral line reproduction

In Figs. 1–3, we present panels of the first four Balmer lines and He II 4541 Å for each target, showing the observation (in black) and the different synthetic spectra resulting from the different methods. Most of the lines are reproduced well by all of the different methods, with good agreement between them. In general, the grid-based fits provide a slightly poorer reproduction of the precise line shapes compared to the tailored and Kiwi-GA approaches, which is expected due to the finite spacing in the parameters.

For the more tailored approaches, the different reproduction of the line profiles does not reflect the ability of a certain code, as evident from the examples where the same underlying code yields different profile shapes. Instead, the panels illustrate the different choices made in the fitting process. This is especially evident when comparing the F2 and F3 results, which use the same code but take a different amount of data into account. Considering for example AzV 377, the profile fits from the Kiwi-GA (F2), which only considers the optical spectrum, are quite similar to manually derived ones from the P1 method. When the UV spectra are also taken into account, the algorithm needs to make a compromise between both sets of data, and the overall fit of the shapes gets slightly worse.

In general, the reproduction of the Ha profiles for the two LMC stars is not satisfactory. The variable nature of the Hα profile in OB supergiants is well known (e.g., Ebbets 1982; Markova et al. 2005; Prinja et al. 2006) and the profiles tend to change even when the other diagnostics do not. Introducing different velocity laws or sophisticated clumping prescriptions with radial dependencies and/or optically thick clumps can improve the spectral fits (e.g., Oskinova et al. 2007; Bouret et al. 2012; Šurlan et al. 2013; Bernini-Peron et al. 2023; Rübke et al. 2023), but in particular pronounced P Cygni profiles in Hα that differ in their velocity diagnostics from the UV profiles are challenging to reproduce. When examining and comparing the results for the different wavelength regimes, one further has to take into account that the available spectra for each object from FUSE, HST, and X-shooter were not taken simultaneously. Consequently, given the intrinsic wind variability, imperfections in the spectral reproductions between the wind-affected diagnostics from different regimes are to be expected. The aim of this study is not to obtain a detailed reproduction of the Hα profile as it does not significantly impact the total set of derived parameters, and would require an extensive additional modeling effort. The main UV diagnostic lines are featured in Figs. 4, 5, and 6, where the observations are compared to the synthetic spectral lines from the different methods. Unlike in the optical regime, the identification of the continuum level in the UV is cumbersome due to the forest of iron lines creating what is sometimes referred to as a “pseudo continuum”. Several methods therefore work with flux-calibrated data in this regime, with the normalization performed afterwards via the employed model. Therefore, discrepancies between observations and model near some of the diagnostic lines are not uncommon. Moreover, interstellar absorption affects the observation. While general reddening is accounted for in the models depicted in Figs. 4, 5, and 6, not all methods have applied line reddening for interstellar Lyα and Lyß absorption on their synthetic spectra. In particular, Lyα absorption can become broad enough to affect the blue edge of N V 1238/1242 Å, limiting its v_∞ diagnostic in some cases. Moreover, no correction for the considerable spectral imprint of the H₂ Lyman and Werner band lines below 1107 Å has been performed, which needs to be taken into account when inspecting the panels for O VI 1032/1038 Å and S IV 1063/1073 Å.

A complete overview of the spectral fitting results is illustrated in Figs. C.1, C.2, and C.3.

Fig. 2 Comparison between the different methods in the Balmer lines and He II λ4541 region for Sk-66° 171.

Fig. 3 Comparison between the different methods in the Balmer lines and He II λ4541 region for Sk-69° 50.

4.2 Prominent spectral discrepancies

Figures 7, 8, and 9 show panels of those lines, which show considerable disagreement between the models in the blue optical region, which is usually the range with the best diagnostics for the stellar parameters of O stars. The depicted lines are He I 4471 Å, the N III multiplet between 4510 and 4524 Å, the N III complex around 4635 Å and 4645 Å, as well as He II 4686 Å. All of these lines are sensitive to the wind onset region, which is one of the most uncertain regimes in stellar atmosphere calculations. In addition to its inherent physical uncertainties (e.g., Cantiello et al. 2009; Sundqvist et al. 2011; Grassitelli et al. 2016; Schultz et al. 2023), there can also be numerical artifacts stemming from connecting the (quasi-)hydrostatic domain to the wind domain described by the ß-law. Consequently, lines that are formed in this onset region are subject to these inherent uncertainties and can for example be affected in their appearance when making minor changes to parameters related to the connection criteria, such as the assumed line broadening in the radiative transfer (u_Dop) or the ß-value. Consequently, He lines should not be trusted blindly as a diagnostic for O stars between 30 and 35 kK, which essentially covers both of our LMC targets in this work. The observed emission of some of the N III lines in Figs. 7, 8, and 9 is also very sensitive to these connection settings. Moreover, there is an overlap of two resonance lines from N III and O III around 374 Å in the extreme ultraviolet (EUV) that affects models at least in the temperature domain between ~33 and 35 kK (Rivero González et al. 2011). Minor details of the modeling approach affect the resulting optical lines, including the wavelengths in the atomic data, the broadening assumed in the radiative transfer, as well as the treatment of line overlaps (see Rivero González et al. 2011; Puls et al. 2020, for a more in-depth discussion). Consequently, there are notable issues in the reproduction of N III with different results for the different stars and no clear preference for any of the methods. The situation gets generally better for the hotter SMC O5 dwarf with the remaining differences in the He lines mainly arising from fixed-grid approaches.

For late O supergiants, He II recombination can either set in or be avoided even when parameters such as T_eff or $\dot{M}$ undergo only minor changes. This effect is likely able to explain some the observed deficiencies for the He II 4686 Å line, in particular for methods such as C2 or P2, which provide a good reproduction of the UV lines. UV and optical diagnostics can in practice favor slightly different temperatures, which cannot be resolved within a given atmosphere code version and thus demand a compromise in the fit, regardless of the specific method used. A similar problem exists for the N III 4512 Å line. While such discrepancies between different wavelength regimes can arise due to the nonsimultaneous observations, they can also reflect limitations in the current treatment of 1D model atmospheres, for example with respect to the assumptions of a single wind velocity and ionization structure.

In the UV range, there is generally a good consensus between the methods taking this regime into account. Larger discrepancies mainly occur for high-ionization lines such as O VI 1032/1038 Å in the case of AzV 377 or N v 1238/1242 for the LMC stars. When calculating atmosphere models with the necessary T_eff derived from the remaining diagnostics, the population of the levels corresponding to these lines is not sufficiently large to reproduce the observed strength. To remedy this shortcoming, the model codes have the ability to include additional X-rays⁶ in the wind, but not all methods make use of this (see Appendix B for the handling of the individual methods). The P V 1118/1128 Å doublet is further known to be sensitive to optically thick clumping (e.g., Sundqvist et al. 2011; Šurlan et al. 2013), which is not taken into account in most of the methods, except F3.

Fig. 4 Comparison between the main UV profiles for AzV 377. The panels from (a) to (h) depict respectively the profiles of O VI 1032/1038 Å, S IV 1063/1073Å (and He II 1085 Å), PV 1118/1128 Å, C III 1176 Å, NV 1238/1242 Å, Si IV 1394/1403 Å, C IV 1548/1551 Å, and N IV 1718 Å. Interstellar Lya and Lyß absorption affects some of the diagnostics; most notably the wing of NV 1238/1242 Å in case of larger terminal velocities (e.g., seen in AzV 377).

Fig. 5 Comparison between the main UV profiles for Sk-66° 171. The spectral windows are the same as in Fig. 4, following the same order.

Fig. 6 Comparison between the main UV profiles for Sk-69° 50. The spectral windows are the same as in Fig 4, following the same order.

4.3 Spectral energy distribution

In Fig. 10, we show the reproduction of the observed spectral energy distribution (SED) from flux-calibrated UV spectra plus optical and IR photometry with the different models for all three targets. Apart from the methods that do not take the UV spectra into account, all methods yield an acceptable reproduction of the SED, illustrating the phenomenon that all O stars essentially appear “blue”, meaning that the flux we see in the optical and beyond simply maps what would be the Raleigh-leans of a blackbody (of ≈ 0.8 … 0.9 T_eff). Nevertheless, hot stars can also deviate significantly from a general black body shape. This is most obvious in the UV, where the large number of iron lines (“iron forest”) forms a “false continuum” and significantly alters the emerging shape of the flux distribution. Extending with their wavelengths far into the usually unobservable EUV, the large number of transitions from iron (and other elements with complex electron configurations) lead to a “blanketing” effect that alters not only the ionization and temperature structure of a hot star (e.g., Dreizler & Werner 1993; Hillier & Miller 1998; Gräfener et al. 2002; Lanz & Hubeny 2003) but also the spectral shape, in that the continuum emission is enhanced at longer wavelengths (e.g., Hummer 1982; Abbott & Hummer 1985). For stars with stronger winds, the slope of the flux decline is further altered by additional free-free emission contributing to the continuum with the relative contributions getting larger for longer wavelengths. Consequently, the temperature determination for hot stars cannot be achieved with photometry (see also Hummer et al. 1988, who in particular discuss the effect of blanketing), but requires a detailed spectroscopic analysis with lines of different ionization stages acting as crucial temperature indicators (cf. the method descriptions in Sects. 3.2, 3.3, 3.4).

To investigate the reproduction of the SED in more detail, we show zoom-ins for the three targets around the photometric magnitudes and the flux-calibrated UV spectra in Figs. 11, 12, and 13. For most of the filters, there is excellent agreement between the model spectra and the photometry. Notable shifts occur in particular for the F1, F2, and F3 methods, which do not take the flux-calibrated UV spectra directly into account, but either use normalized parts of the UV spectrum (F3) or do not consider the UV part of the spectrum. To derive the luminosity, these methods use anchor magnitudes and a reddening law as described in Appendix B.3. This seems to lead to a slight overestimation of the UV flux for the hottest target in the sample. Nevertheless, while the shift in Fig. 11 appears quite dramatic, this is mostly a result of employing an extinction law that is not adjusted for the UV regime, and the difference in the derived luminosity is less than 0.1 dex compared to the other methods, which is a typical error margin (see also Table3).

Fig. 7 Comparison between the lines with more prominent disagreement on the optical for AzV 377.

Fig. 8 Comparison between the lines with more prominent disagreement on the optical for Sk -69° 50.

Fig. 9 Comparison between the lines with more prominent disagreement on the optical for Sk-66° 171.

Fig. 10 Modeled and observed SED of the targets. The flux points correspond to the magnitudes listed in Table 1, following the same order (from bluer to redder). The observed flux-calibrated UV spectra correspond to those acquired by ULLYSES.

Fig. 11 Zoom-in comparison of the SEDs from the different models for AzV 377 around the applied photometry (small upper panels, crosses mark photometric measurements) and the flux-calibrated UV spectra (big lower panel).

Fig. 12 Same as Fig. 11, but for Sk-66° 171.

Fig. 13 Same as Fig. 1, but for Sk-69° 50.

4.4 Abundances

Unless one purely relies on a fixed grid of models, the determination of the elemental abundances usually requires additional rounds of iteration among the necessary model calculations. In OB stars, the temperature in many cases cannot be sufficiently constrained without taking metal lines of different ionization stages into account. Consequently, abundance effects can overlap with temperature effects and in several (though not all) methods the finer tuning of the abundances is only performed once the main stellar parameters are robustly determined. Moreover, some of the involved grids only have a fixed set of abundances. Given that we do not aim to further “tune” the derived values after comparing our initially obtained values, we do not expect our abundances to be as robust as they usually would be in studies focusing on particular stars. Still, we can identify general trends and discrepancies between the analysis methods.

For both of the LMC O supergiants, almost all methods yield a He enrichment. Almost unanimously, all methods predict a He mass fraction of ~0.35 for the O7(n)(f)p target, while the scatter is larger for the O9 Ia star with values reaching from almost zero enrichment (X_he = 0.26) up to X_he = 0.39. For the SMC O5 SMC dwarf, there is a similar scatter, interestingly now with different approaches yielding the higher enrichment of up to X_He = 0.40. Unless hydrogen is strongly depleted, the imprints of He enrichment are more subtle, making the determination more cumbersome than that of other abundances such as CNO.

All methods that determined CNO abundances found strong nitrogen enrichment for all of the studied targets. Converted to mass fractions, the LMC baseline abundances (cf. Vink et al. 2023) are X_c = 9.06 • 10⁻⁴, X_n = 1.11 • 10⁻⁴, and X_o = 2.96 • 10⁻³, yielding a combined CNO abundance of 3.98 • 10⁻³. For the SMC, values are X_c = 2.34 • 10⁻⁴, X_n = 0.47 • 10⁻⁴, and X_o = 1.33 • 10⁻³, yielding a total CNO mass fraction of 1.61 • 10⁻³. With nitrogen mass fractions between 7.4 • 10⁻⁴ (F2) and 1.9 • 10⁻³ (C1), enrichment factors between 15 and 40 are found for the SMC dwarf AzV 377. The situation is less clear for carbon, while the depletion of oxygen is clearly confirmed. The total CNO abundance found for AzV 377 scatters between 0.9 and 2.0 times the baseline value, preventing any more robust conclusions as to whether our sample star is actually slightly more metal rich than is presumed to be typical for the SMC.

A similar scatter around the total CNO baseline abundances is found for the two LMC targets, with Sk-69° 50 yielding slightly higher factors (0.8 … 2.4) than Sk-66° 171 (0.78 … 1.96). Nevertheless, the two targets are quite different in their nitrogen enrichment, which is found to be much higher for the O7(n)(f)p target, where most methods yield enrichment factors of ~40 (except C2) compared to more moderate factors of ~10 for the O9 supergiants. Clearly, Sk-69° 50 seems to be the most evolved target in our sample, as all methods find carbon to be depleted, while the situation is less clear for the other two sample stars.

4.5 Ionizing fluxes

All of our sample stars show a considerable flux beyond the hydrogen ionization edge. As also evident from the tabulated results (Tables C.1, C.2, and C.3), the ionizing fluxes depend strongly on the temperature and luminosity of the stars with higher temperatures and luminosity yielding higher fluxes. The fluxes beyond the He II ionization edge are more complicated, as they also strongly depend on the wind density. For denser winds, even very hot stars can yield essentially no He II ionizing flux. This is the case for the two supergiants in the sample, for which the models formally yield photon fluxes of up to ~ 10⁴² s⁻¹. These values are orders of magnitude below considerable contributors such as early O dwarfs (e.g., Smith et al. 2002; Martins & Palacios 2021), hot, thin-wind Wolf-Rayet stars (e.g., Crowther & Hadfield 2006; Sander et al. 2023), or luminous envelope-stripped stars below the WR regime (e.g., Götberg et al. 2023; Ramachandran et al. 2023). We note that the absolute numbers of the reported magnitude must be taken with care as these low values can be subject to numerical uncertainties, for example if models are optically thick up to the outer boundary at some of the corresponding wavelengths. The studied O5 dwarf has a He II-ionizing photon flux of ~10⁴³… 10⁴⁴ s⁻¹, which is actually in line with expectations for the derived HRD position (Martins & Palacios 2021).

Fig. 14 HRD with the obtained positions for the two LMC stars Sk -69° 50 and Sk -66° 171 – indicated by the respective labels. For comparison, tracks from Brott et al. (2011, up to Minit = 50 M⊙) and Köhler et al. (2015, from 60 M⊙) are shown.

5 Discussion

5.1 HRD position and evolutionary status

In Figs. 14 and 15, we provide an overview of all the obtained positions of our three sample stars with the different methods in the Hertzsprung-Russell diagram (HRD). For all targets, there is a noticeable trend that hotter solutions tend to come with higher luminosities. This can be understood when considering the SED fit (cf. Sect. 4.3). If one aims to fit the same photometric SED with a hotter model atmosphere, this enforces a slightly higher reddening and thus a higher luminosity⁷. In particular, we see the F2 and F3 methods yielding the highest luminosities (and P1/P2 the lowest), in line with our findings for the SED fits (cf. Sect. 4.3). The obtained stellar parameters are therefore not independent and changes in one parameter can propagate into other quantities. Given the high number of input parameters into stellar atmosphere models and the nontrivial effects of their variation, only the most obvious ones can be calculated in the form of a rigorous error propagation. In all other cases, the only feasible option is to assume a larger error than obtained from statistical considerations. As for cool-star atmospheres, systematic errors are usually not considered at all. These can arise for example because of uncertain atomic data, method-inherent approximations, or code-specific numerical treatments.

Fig. 15 HRD with the obtained positions for the SMC star AzV 377. For comparison, tracks from Brott et al. (2011, up to Minit = 60 M⊙) and Köhler et al. (2015, from 80 M⊙) are shown.

5.2 Mass discrepancy

From comparing the derived positions in the HRD, one can derive “evolutionary masses” M_evol, assuming that a given set of tracks (and the interpolation between them) sufficiently describe the history of the stars. Using the tracks from Brott et al. (2011) and Köhler et al. (2015), also shown in the HRDs (Figs. 14, 15), we derived M_evol for each of our three sample stars and compare them to the spectroscopically derived mass M_spec in Fig. 16, employing the same χ² approach as in Bernini-Peron et al. (2023). The age-dependent M_evol is generally lower than the zero-age main sequence mass M_init because of the decreasing mass along each evolutionary track. In all cases, the ratio between M_evol and M_spec is ≥1.0, meaning that the masses inferred from the evolutionary tracks are similar to or higher than the ones from spectroscopy. This so-called “mass discrepancy” is a long-standing issue in the analysis of hot stars (e.g., Herrero et al. 1992; Markova et al. 2018). Interestingly, detailed spectral analyses for detached pre-interaction binaries (e.g., Mahy et al. 2017, Mahy et al. 2020) obtained spectroscopic masses that match evolutionary estimates, raising questions about whether or not the evolutionary status presumed in the tracks applies to all of our analyzed sample stars.

With the determination of the radius R_2/3 from L and T_eff, the different luminosities for example also affect the derived spectroscopic masses. The highest ratio occurs for the F4 analysis of Sk -66° 171 and is a consequence of the HRD “outlier” position, where F4 yielded a much lower luminosity than the other methods. Disregarding this point, the ratios for AzV 377 and Sk -66° 171 are relatively moderate, with values ranging between only 1.0 and ~1.6. In particular, we note that for the late O supergiant Sk-66° 171 the more tailored methods (F3, C2, P2) yield the best matches. In contrast, the other LMC star, Sk-69° 50, has a systematic shift and never shows a ratio below 1.5. It is therefore likely that this star is not sufficiently described by the evolutionary tracks. Moreover, the origin of the class of Onfp stars of which Sk-69° 50 is a member has been subject to speculation, including the suggestion that these objects are products of stellar mergers (Walborn et al. 2010). For the SMC O5 dwarf AzV 377, the evolutionary situation is less clear, with the methods scattering between good agreement and notable discrepancy.

To remove the mass discrepancies, the derived log g values would have to be larger by 0.1 to 0.3 dex. Even when neglecting Sk -69° 50 due to its probably more evolved status, the necessary increase would have to be up to 0.18 dex to account for a mass discrepancy factor of 1.5. One ingredient that could increase log g is the inclusion of a turbulence term in the hydrostatic equation (cf. Eq. (4)), but so far only one of the codes applied in this work (PoWR) does that. Markova et al. (2018) studied the mass discrepancy of Galactic O-type stars using both CMFGEN and FASTWIND, finding comparable discrepancy values with the two codes. Hence, while it is too early to derive a clear tendency and there are prominent exceptions such as the C2 method for AzV 377 yielding only a small discrepancy (cf. Fig. 16), the inclusion of a turbulence term in the hydrostatic equation and its resulting increase in log g could mark an important step to minimize the occurrence of mass discrepancies. A more focused study on the inclusion of different microturbulent velocity values in the hydrostatic equation – further motivated by recent 2D simulation results from Debnath et al. (2024) – would be necessary to better judge this effect.

Fig. 16 Ratio of the determined spectroscopic masses to evolutionary masses (based on Brott et al. 2011 or Köhler et al. 2015 respectively) for the different methods and sample stars.

Fig. 17 Comparison of the derived (F1, F2: assumed) terminal wind velocity v∞ versus the derived effective temperature. The dashed lines denote the LMC (gray) and SMC (salmon) relations from Hawcroft et al. (2024).

5.3 Wind parameters

5.3.1 Terminal velocities

With the direct availability of the terminal wind velocity from the UV spectra, all of the methods making use of these data obtain very similar values for v_∞. Notably, all of our targets show terminal velocities between ~1800 and 2000 km s⁻¹. This is not expected from their spectral types, as evident also from Fig. 17, where we plot the derived values for v_∞ as a function of T_eff and compare them to the trends derived in Hawcroft et al. (2024, paper III of the XShootU series). The methods assuming a terminal velocity (F1, F2, see Appendix B.4) consequently overestimate v_∞, for the SMC O5 dwarf AzV 377, while the two LMC stars are closer to the expected relation. The low terminal velocity for the SMC star is surprising given its high temperature and – as we discuss below – does not coincide with a higher mass-loss rate.

The values for the two LMC stars align well with the v_∞(T_eff) trend reported by Hawcroft et al. (2024). In particular, the value for Sk -69° 50 matches perfectly, while the value for the late O supergiant Sk -66° 171 is a slightly higher than the value from the trend formula. To compare our findings with the predictions from the mCAK theory for radiation-driven winds of OB stars (Castor et al. 1975; Pauldrach et al. 1986), we also plot the ratio of v_∞ to the effective escape velocity, $$v_{\text{esc},\Gamma }:=\sqrt{\frac{2GM}{R}\left(1-{\Gamma }_{\text{e}}\right)},$$(13)

in Fig. 18. Usually, Eq. (13) is evaluated at R = R_2/3, with M being the stellar mass, Γ_e ∝ L/M denoting the classic Eddington parameter taking only electron scattering opacity into account, and G the gravitational constant. The determination of v_esc,Γ is subject to a variety of error propagations, resulting from the uncertainties in for example the determination of log g and the luminosity L. Consequently, the derived v_∝/v_esc,Γ ratios in Fig. 18 show considerable error bars, with the smallest bars actually resulting from incomplete error estimations (e.g., for C2). Considering the large error bars, one could argue that the obtained values are not in conflict with the presumed ratio of 2.65 times the escape ratio found by Kudritzki & Puls (2000), who slightly updated the value from the factor 2.6 found by Lamers et al. (1995). However, there is a systematic trend for the two LMC stars towards higher ratios, which is qualitatively in line with the mass discrepancy trend found, namely the tendency to determine lower spectroscopic masses for these stars than what would be inferred from evolutionary tracks. For the SMC dwarf, the opposite trend is seen, with the two methods determining v_∞, yielding ratios of ~2.1. This actually aligns nicely with the expected metallicity dependence of v_∞ ∝ Z^0.2 for v_∞ from Vink & Sander (2021) and Hawcroft et al. (2021). Assuming 0.5 Z_⊙ for the LMC and 0.14 Z_⊙ for the SMC, the ratio of 2.6 is expected to reduce to 2.02, which is well in line with the findings for AzV 377.

Fig. 18 Ratio between the terminal wind velocity v∞ and the effective escape velocity vesc,Γ as a function of Teff. The dashed gray line shows the empirical results from Kudritzki & Puis (2000).

5.3.2 Mass-loss rates

The mass-loss rates determined by the different methods are shown in Fig. 19, where we also plot three predictions from the literature, namely Vink et al. (2001), Krtička & Kubát (2018), and Björklund et al. (2023). For F1, the reported value should be considered as an upper limit due to the use of unclumped wind models⁸. To account for the fact that the formulae by Vink et al. (2001) and Björklund et al. (2023) have additional dependencies in addition to those on luminosity and metallicity, we are shading areas that account for the maximum and minimum values of the remaining parameters. With the exception of the model F2, which does not account for the UV spectrum, all methods determine mass-loss rates for the SMC dwarf AzV 377 that fall even below the Björklund et al. (2023) predictions. This seems to be in line with earlier findings of dwarfs in the SMC by Bouret et al. (2003), Ramachandran et al. (2019), and Rickard et al. (2022).

For the LMC targets, the situation is different. The mass-loss rate derived for the late-O supergiant Sk -66° 171 agrees with the Vink et al. (2001) predictions, but also with the Björklund et al. (2023) recipe as this formula turns upwards for higher luminositites. For Sk-69° 50, there are two groups of solutions resulting from the assumption of either low and moderate clumping (f_cl ≤ 5) coinciding with larger mass-loss rates, or strongly clumped solutions (f_cl ≥ 10) and correspondingly lower values for $\dot{M}$. Depending on assumptions or results for clumping parameters and stratification, the derived mass-loss rates are either about 0.3 dex lower than predicted by Vink et al. (2001) or even slightly higher than the prediction. The high-clumping solutions for Sk-69° 50 further align with the Björklund et al. (2023) predictions. In its simple $\dot{M}$(L) form, that is, without any temperature dependency, the Krtička & Kubát (2018) formula is not able to reproduce the derived values and stays between the LMC and SMC solutions. Within this very limited sample, it is hard to draw any robust conclusions about the mass-loss rates, but the results underline the complexity of the situation with our different sample stars spanning from an SMC star with a relatively weak wind to an O7 target that probably exceeds the wind expectations for the LMC.

In addition to the raw mass-loss rate, it is worth also considering the modified wind momentum rate, $$D_{\text{mom}}=\dot{M}{v}_{\infty }\sqrt{R/R_{\odot }},$$(14)

which is expected to be proportional to some (positive) power of the stellar luminosity L/L_⊙ (e.g., Kudritzki & Puls 2000). In Fig. 20, we show the modified wind momentum rates D_mom, taking into account the dumping-adjusted mass-loss rates $\dot{M}\sqrt{f_{\text{cl}}}$, rather than their raw values. By doing so, we lift the split between the two groups seen for Sk-69° 50 in Fig. 19 and both LMC targets now get very similar values in D_mom from most of the methods, with the remaining spread being mainly due to the differences in determining log L. The results for the SMC dwarf, in contrast, remain slightly more spread, but we have to consider that only two methods (F3 and P1) performed a full UV+optical study for this target and the two results agree within their error bars. For comparison, we also show results from Bouret et al. (2013) for the SMC as well as Mokiem et al. (2007) and Ramachandran et al. (2018b) for the LMC. The data compiled in Mokiem et al. (2007) also contain analysis results from Crowther et al. (2002) and Massey et al. (2005). In general, our derived quantities fall within the obtained literature values.

Finally, we also plot the transformed mass-loss rate, $${\dot{M}}_{\text{t}}=\dot{M}\sqrt{f_{\text{cl}}}\frac{1000 \text{km} {\text{s}}^{-1}}{v_{\infty }}{\left(\frac{{10}^6{L}_{\odot }}{L}\right)}^{3/4},$$(15)

defined by Gräfener & Vink (2013) as a function of L/M in Fig. 21. The quantity describes the mass-loss rate the star might have if it had 10⁶ L_⊙ and an unclumped wind with v_∞ = 1000 km s⁻¹. Similar to the D_mom plot (Fig. 20), the SMC dwarf ends up in a very different regime than the two LMC targets, which cluster more than in the D_mom plane. While the absolute values for $\dot{M}$ are relatively close to the regime obtained in Wolf-Rayet studies (Sander & Vink 2020; Sander et al. 2023), neither of our two LMC targets would qualify as an Of/WN star as this would require Hβ to show a P Cygni profile (Crowther & Walborn 2011), which is not observed. The two LMC stars are also too cool to yield notable He II ionizing flux (cf. Sect. 4.5), reflecting that the characteristic values for classical WR stars (−4.5) are not transferable to this parameter regime.

Fig. 19 Mass-loss rate $\dot{M}$ versus stellar luminosity L for our sample stars analyzed with the different methods. For comparison, we also plot the mass-loss recipes from Vink et al. (2001, solid shading), Krtička & Kubát (2018, dashed lines), and Björklund et al. (2021, hatched shading). The SMC comparison data are drawn in light red.

Fig. 20 Modified wind momentum rate versus stellar luminosity L for our sample stars analyzed with the different methods.

Fig. 21 Transformed mass-loss rates versus the ratio between luminosity L and spectroscopic mass M for our sample stars analyzed with the different methods.

6 Conclusions and perspectives for forthcoming papers

In this work, we present an analysis of three O stars from the ULLYSES and XShootU sample, namely the SMC O5V((f)) dwarf AzV 377 as well as the LMC O stars Sk-69° 50 (O7(n)(f)p) and Sk-66° 171 (O9Ia). We analyze these targets using a variety of different methods, applying different model atmosphere codes (FASTWIND, CMFGEN, PoWR), ranging from grid-based approaches to tailored spectral fits. Some methods are only applied to some of the targets and some only take the optical spectra into account, thereby skipping a detailed determination of the wind parameters only accessible from the UV. This study was performed as a “blind test”, meaning that each method was performed without prior knowledge of any outcomes from the other methods. The study is not intended to be a benchmarking of any particular atmosphere code or method, but our aim is to provide an overview of the “natural” spread of results obtained with the different existing approaches in the field. (No fine tuning of the results was performed after comparing the resulting parameters.) Moreover, our detailed descriptions of the individual methods serve as an introduction of the different techniques applied within the “XShooting ULLYSES” collaboration.

Overall, the different applied methods show a reasonable amount of agreement for the obtained parameters. Nevertheless, a spread of up to 3 kK is obtained for the effective temperatures across all three targets, with the GA-based method tending to yield slightly higher values. The differences in log g are on the order of 0.1 dex for the SMC dwarf and up to 0.2 dex for the LMC O stars. One ingredient to minimize the log g discrepancies could be the inclusion of a microturbulent velocity ξ in the hydrostatic treatment, which is so far only possible in one of the three atmosphere codes employed in this work (PoWR). The inclusion of ξ > 0 demands higher log g values and thereby also has the potential to reduce current discrepancies between spectroscopic and evolutionary masses (hinted, e.g., by the good agreement for Sk-66° 171). However, we also obtain some spread in this “mass discrepancy” between methods that do not include turbulence in the hydrostatic solution, and therefore more in-depth studies on this topic are required. Differences in the adopted E (B – V) values tend to be up to 0.1 dex and can affect the derived luminosities. The spread in reddening is mainly caused by the genetic algorithm methods, which tend to yield higher reddening values and consequently also higher luminosities. The remaining methods, regardless of the applied atmosphere code, differ by less than 0.05 dex in E(B – V).

The inspection of our three sample targets illustrates the wide range of regimes found among the O stars in the ULLYSES sample. The SMC dwarf AzV 377 exhibits a very low mass-loss rate, below current theoretical predictions, while the two LMC targets show comparably strong winds aligning with or even exceeding current recipes. The wind and the abundance analysis as well as the derived mass discrepancy suggest that the O7(n)(f)p star Sk-69° 50 is a particularly evolved object that might no longer be in the stage of central hydrogen burning. The two other targets are also nitrogen-enriched, but the methods differ in their conclusion regarding a potential mass discrepancy.

Overall, the different methods scatter slightly more than expected, which we attribute in particular to the “blind test” scenario. Several of the heterogeneous choices among the various methods in the current study (e.g., the selection of the clumping law or differences in the reddening approach) could be avoided in a more coordinated effort. Hence, the scatter obtained in this work should reflect the amount of scatter to be expected when “blindly” combining data from different literature sources. For coordinated efforts within the collaboration, we recommend harmonize assumptions and selections as much as possible in order to minimize uncertainties in future studies. Our derived value spreads also provide a good indicator for future comparisons of data from heterogeneous sources.

Acknowledgements

This manuscript has been inspired by discussions during the workshop “ULLYSES – new horizons in massive star spectroscopy”, hosted and supported by the Lorentz Center, Leiden, Netherlands. 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. AACS, MBP, VR, RRL, and CJKL are supported by the Deutsche Forschungsgemeinschaft (DFG – German Research Foundation) in the form of an Emmy Noether Research Group - Project-ID 445674056 (SA4064/1-1, PI Sander). AACS and RRL acknowledge funding by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) – Project-ID 138713538 – SFB 881 (“The Milky Way System”, subproject P04). AACS, MBP, VR, and ECS are further supported by funding from the Federal Ministry of Education and Research (BMBF) and the Baden-Württemberg Ministry of Science as part of the Excellence Strategy of the German Federal and State Governments. ECS acknowledges financial support by the Federal Ministry for Economic Affairs and Climate Action (BMWK) via the German Aerospace Center (Deutsches Zentrum für Luft- und Raumfahrt, DLR) grant 50 OR 2306 (PI: Ramachandran/Sander). CJKL gratefully acknowledges support from the International Max Planck Research School for Astronomy and Cosmic Physics at the University of Heidelberg in the form of an IMPRS PhD fellowship. PAC and JMB are supported by the Science and Technology Facilities Council research grant ST/V000853/1 (PI. V. Dhillon). DP acknowledges financial support by the Deutsches Zentrum für Luft und Raumfahrt (DLR) grant FKZ 50OR2005. RK acknowledges financial support via the Heisenberg Research Grant funded by the German Research Foundation (DFG) under grant no. KU 2849/9. JSV gratefully acknowledges support from STFC via grant ST/V000233/1. OM acknowledges support from the project RVO:67985815 of the Astronomical Institute of the Czech Academy of Sciences. A.u.-D. acknowledges NASA ATP grant number 80NSSC22K0628 and support by NASA through Chandra Award number TM4-25001A issued by the Chandra X-ray Observatory 27 Center, which is operated by the Smithsonian Astrophysical Observatory for and on behalf of NASA under contract NAS8-03060. ACGM acknowledges the support from the Polish National Science Centre grant Maestro (2018/30/A/ST9/00050) and from the Max Planck Society through a “Partner Group” grant. F.N. acknowledges support by grants PID2019-105552RB-C41 and PID2022-137779OB-C41 funded by MCIN/AEI/10.13039/501100011033 by “ERDF A way of making Europe”. SRB acknowledges support by NextGeneration EU/PRTR and MIU (UNI/551/2021) through grant Margarita Salas-ULL. This work used the Dutch national e-infrastructure with the support of the SURF Cooperative using grant no. EINF-3257.

Appendix A Detailed physical treatments within the different expanding atmosphere codes

A.1 Velocity and density stratification

A.1.1 Fastwind

All Fastwind models in this study employ a ß-velocity law as mentioned in Sect. 3.1, which is smoothly connected to the hydrostatic regime at v ≈ 0.1 a_s(T_eff), with a_s denoting the (isothermal) sound speed. The ß-law is implemented in the form of Eq. (1), with ß and v_∞ given as input parameters. The quantity R₀ is usually labeled b in their literature and fixed by the connection demand. To account for the radiative acceleration in the solution of the hydrostatic equation, in the very first iterations FASTWIND approximates the flux-weighted opacity by the Rosseland opacity, which is further approximated by a Kramer’s-like formula, $${\kappa }_{\text{Ross}}(r)\approx {\sigma }_{\text{e}}{n}_{\text{e}}(r)\left[1+k_{\text{c}}\rho (r)T{(r)}^{-x}\right]$$(A.1)

(Santolaya-Rey et al. 1997), for the Rosseland opacity κ_Ross with n_e denoting the electron number density and σ_e the Thomson cross section. The parameters k_c and x are obtained from fits to the actually computed κ_Ross. These approximations are relaxed in the course of the non-LTE iterations. After an initial phase, the Rosseland opacity is calculated from its actual definition and the radiative acceleration is updated using the current flux-weighted opacities.

A.1.2 CMFGEN

For the initial setup of the (quasi-)hydrostatic layers, either a TLUSTY or an old CMFGEN model can be used. To jointly account for the supersonic and hydrostatic parts, CMFGEN implements a velocity law of the form $$v(r)=\frac{v_0+\left(v_{\infty }-v_0\right){\left(1-\frac{R_{*}}{r}\right)}^{\beta }}{1+\frac{v_0}{v_{\text{core }}}\exp \left(\frac{R_{*}-r}{h_{\text{eff }}}\right)},$$(A.2)

as introduced in Hillier et al. (2003) for this single-ß form. The “effective scale height” h_eff can be updated during the calculations with the current radiative acceleration to closely fulfill the hydrostatic equation of motion (e.g., Martins et al. 2012b). Other approaches analytically calculate v(r) without using the ß-law by means of the W-Lambert procedure (Müller & Vink 2008; Gormaz-Matamala et al. 2021), but it requires a significant effort to determine the extra parameters.

A.1.3 PoWR

The PoWR models employed in this work use the default PoWR branch with a pre-described ß-law in the outer part instead of PoWR^HD where the entire velocity field v(r) can be obtained from solving the hydrodynamic equation of motion (e.g., Gräfener & Hamann 2005; Sander et al. 2017). The ß-law implementation in PoWR reads $$v(r)=p_1{\left(1+\frac{R_{20}}{r+p_2}\right)}^{\beta }.$$(A.3)

The two parameters (p₁, p₂) are determined at the beginning of the model and after updates of the density stratification. The latter happens either when the hydrostatic equation is not sufficiently fulfilled or when the optical depth at the inner boundary deviates notably from the pre-specified value. In the default setting, the two parameters p₁ and p₂ are adjusted such that v(R_max) = v_∞ (i.e., p₁ ~ v_∞) and one obtains a smooth connection of the velocity gradient dv/dr when the analytic Eq. (A.3) is merged with the numerically solved v(r) for the hydrostatic part (cf. Sander et al. 2015). Alternatively, the connection can be forced at a specific fraction of the sound speed, which usually yields a different value for the parameter p₂. For large ß-values, the default gradient connection can lead to an unfavorable situation as the gradients for larger ß-laws could enforce a connection too far beneath the sonic point, which can in some cases impact the derived absorption line profiles. In such cases, the forced sonic connection is preferred. In the current work, ß-values around 1 are employed in the PoWR models and thus the selection of the connection criterion is of no importance. However, a good fit of the LMC targets required the use of a 2ß-law, where Eq. (A.3) is extended to $$v(r)=p_1{\left(1+\frac{R_{20}}{r+p_2}\right)}^{{\beta }_1}+p_3{\left(1+\frac{R_{20}}{r+p_4}\right)}^{{\beta }_2}.$$(A.4)

Instead of a single ß-parameter, now two values ß₁ and ß₂ are required with ß₂ > ß₁. Moreover, a weighting factor q_2ß for the impact to v_∞ from the second term needs to be provided, such that p₁ ≈ v_∞ (1 – q₂ß) and p₃ ≈ v_∞q_2ß. The calculation of p₂ is similar to the single-ß case and p₄ = R₂₀ − R_con follows from the connection radius R_con determined as described above, as the determination of R_con uses only the gradient from the ß₁ term.

A.2 Wind clumping

A.2.1 Fastwind

To account for wind clumping, Fastwind can use two different approaches. In the standard “microclumping” approximation, the wind is composed only of clumps with a void “inter-clump medium”. The clumps are further optically thin and a “clumping factor” f_cl = ⟨ρ²⟩/⟨ρ⟩² needs to be defined. Under the above conditions, f_cl = D, where D describes the density contrast between the clumps and the mean density. Hence, for a wind with mean density ⟨ρ⟩, the clumps as such have a density D⟨ρ⟩. The volume filling factor f_V can then simply be calculated as $f_{\text{V}}=f_{\text{cl}}^{-1}$. For f_cl = 1 , a smooth (“unclumped”) wind situation is recovered. Alternatively, a more extensive effe tive opa ity formalism an be employed to statisti ally a ount also for opti ally thi k clumps (Sundqvist & Puls 2018). The formalism introduces three additional parameters, f_ic, f_vel, and the “porosity length” h. For the latter, the default scaling in Fastwind is h(r) = h_∞v(r)/v_∞ with h_∞ = R_2/3 (Sundqvist & Puls 2018; Brands et al. 2022). This scaling is employed in all models employing the detailed clumping formalism in this work. With values of f_ic > 0, a non-void inter-clump medium can be introduced. Similar to D, the quantity f_ic = ρ_ic/⟨ρ⟩ describes the ratio of the density to the mean density, but now for the inter-clump medium. With f_ic > 0, the volume filling factor reads $$f_{\text{V}}=\frac{{\left(1-f_{\text{ic}}\right)}^2}{f_{\text{cl}}-2f_{\text{ic}}+f_{\text{ic}}^2}.$$(A.5)

The mean density is then $$\langle \rho \rangle =f_{\text{V}}{\rho }_{\text{cl}}+\left(1-f_{\text{V}}\right){\rho }_{\text{ic}}.$$(A.6)

The parameter f_vel describes a velocity filling factor that can account for the so-called “velocity porosity” (Owocki 2008).

The value of f_vel enters in the calculation of the line optical depth for the clumps (Sundqvist & Puls 2018), which is proportional to the Sobolev optical depth of the mean wind velocity field divided by the factor f_vor := f_vel/(1 − f_vel). For f_vel → 0, the line optical depth of the clumps approaches infinity, and the effective opacity is determined by the inter-clump medium alone.

The FASTWlND-default clumping stratification is set by the two parameters v_cl,start and v_cl,max. Below v_cl,start, a smooth atmosphere is assumed. Between v_cl,start and v_cl,max, the clumping factor is increased linearly until the specified value of f_cl is reached at v_cl,max. Beyond this point, the constant f_cl is employed. The formalism from Sundqvist & Puls (2018) is an effort to describe two wind components without introducing significant computational overhead. Therefore, such models cannot account for a difference in the level populations (and ionization stages) between clumps and inter-clump medium at a given radius as this would require the calculation of an additional set of non-LTE population numbers and come with a significant computational cost.

A.2.2 CMFGEN

Clumping is treated in the “microclumping” approximation, assuming a void inter-clump medium. In all CMFGEN models presented in this work, the volume filling factor f_V has a velocity-dependent behavior with $$f_{\text{V}}(r)=f_{\text{V},\infty }+\left(1-f_{\text{V},\infty }\right)e^{-v(r)/v_{\text{cl}}},$$(A.7)

as introduced in Hillier et al. (2003). The quantity f_V,∞ denotes the value at r — ∞. For a void inter-clump medium and optically thin clumping, we obtain D = f_cl and $D_{\infty }=f_{\text{V},\infty }^{-1}$ (see Sect. 3.1.2). The free parameter v_cl is a characteristic velocity varied in a tailored analysis to describe the onset of the clumping. Due to the exponential increase in Eq. (A.7), this does not imply a completely unclumped wind for v < v_cl, contrary to the meaning of v_cl in FASTWIND.

A.2.3 PoWR

The “microclumping” approximation is the default setting in PoWR (Hamann & Koesterke 1998), assuming a clumped wind with an overdensity Dρ compared to a smooth wind with the density ρ. Similar to CMFGEN, the inter-clump medium is assumed to be void, meaning that the volume filling factor is f_V = D⁻¹. The clumping factor D = f_cl can be radially dependent with a multitude of choices for the prescription of D(r). For all models in this work, the radius-based prescription $$D(r)=\{\begin{array}{ll}{D}_{\text{in }}\hfill & \text{ for }r<R_1\hfill \\ \frac{D_{\text{in }}}{2}[1+\cos (x)]+\frac{D_{\infty }}{2}[1-\cos (x)]\hfill & \text{ for }{R}_1<r<R_2\hfill \\ D_{\infty }\hfill & \text{ for }r>R_2\hfill \end{array}$$(A.8)

is chosen, with $$x=\frac{\pi }{R_2-R_1}\left(r-R_1\right)$$(A.9)

and D_∞ denoting the clumping factor reached at the outer boundary. The inner value is set to D_in = 1. Between the radii R₁ and R₂ ≥ R₁ the clumping factor increases from D_in to D_∞ in the shape of a cosine. In all models presented in this work, R₁ is set to the radius of the sonic point. The outer value differs for each target to improve the spectral fit. The values for R₂ are listed as R_D in Tables C.1, C.2, and C.3. PoWR also has options to account for optically thick clumping or detailed rotational broadening of wind lines in the calculation of the emergent spectrum (Oskinova et al. 2007; Shenar et al. 2014), but we do not make use of either feature in the present paper.

Fig. A.1 Comparison of the different radial behaviors of the clumping factor fcl = D, assuming optically thin clumping and input values used for the analysis of Sk-66° 171. To focus on the clumping stratification, all curves assume the same maximum clumping factor (fcl = 20) and the same underlying velocity field (v∞ = 1800 km s−1 ).

With the clumping parametrization used in the PoWR models applied in this work, the maximum clumping factor D_∞ is reached much further out than in FASTWIND and CMFGEN models employing clumped winds. Moreover, the increase from the sonic point starts rather slow if described by Eq. (A.8) due to the differences between radius and velocity space. Both effects are illustrated in Fig. A.1 and have to be considered when interpreting and comparing the resulting values. The PoWR example shown in Fig. A.1 uses R_D = 3 R₂₀, corresponding to a velocity of ~ 900 km s⁻¹.

Appendix B Detailed method descriptions

In this section we provide a detailed description of the different approaches, sorted by the different physical aspects. While some individual results are mentioned, a complete overview of all derived values is provided in Tables C.1 to C.3 with a more broad discussion of the results and their implications given in Sects. 4 and 5.

B.1 Rotation and macroturbulence

F1. To determine the projected rotational and macroturbulent velocities, iacob-broad is employed. Different lines are used for the different stars as metal lines depend on the spectral type and they are faint at low metallicities. N IV λ 5200, 5204, 6380 are used for AzV 377; O III λ 5592 and C IV λ 5812 for Sk-69° 50; and Si IV λ 4089, N III λ 4510, 4514, O III 5592 and C IV λ 5801, 5812 for Sk-66° 171. The obtained projected rotational velocities are, respectively, v sin i= 74±7, 196+20 and 97±10km s⁻¹, where we adopt a 10% uncertainty when less than four lines are available. v_mac is not well constrained, and we adopt estimated values of 60 km s⁻¹ for AzV 377 and Sk -69° 50 and 30km s⁻¹ for Sk-66° 171.

F2/F3. In the optical-only analyses (F2) the projected rotational velocity v sin i is determined by convolving the spectra within the Kiwi-GA run (cf. Appendix B.4). In the optical + UV analyses (F3) the rotational velocity is not fitted, but instead fixed at the best fitting value as found from the optical-only analyses (F2). No explicit macroturbulence is assumed (v_mac = 0km s⁻¹).

F4. The rotational velocities are obtained within the spectral fitting process. For this, the synthetic spectra from the calculated model grids are pre-convolved with a set of projected rotational velocities spanning v sin i = 0, 20, 50, 100, 150, 200, 250, 300, and 400 km s⁻¹ (see also Bestenlehner et al. 2024). For all targets, a macro-turbulent velocity of v_mac = 20km s⁻¹ is adopted.

C1. We determine v sin i with the Fourier transform analysis from iacob-broad. For this, we used He I 4713, O III λ 5592 and C IV λ5801, and computed an overall mean for v_sin i. We estimate that this approach yields a typical error of ± 10 km s⁻¹ on v sin i. Broadening caused by macroturbulence proved uncertain to constrain with this method, and it was derived from a comparison of synthetic spectra with observations. A reference spectrum is chosen from the C1 initial model grid where T_eff and log g are close to the expected values for the star’s spectral type and luminosity class. The spectrum is convolved with the instrumental resolution (assuming a Gaussian profile for the instruments) and the radial velocity is obtained from direct comparison with the observed spectrum. The error on v_rad with such a simple approach is ± 10 km s⁻¹. The synthetic spectra are then convolved with the broadening parameters (v sin i and v_mac, assuming a radial-tangential profile for the latter), plus instrumental broadening, before comparing them to observations. The shape of line wings gives access to macroturbulence in case v sin i is not too large. For Sk -69° 50 and Sk -66° 171, we adopted v_mac = 80 km s⁻¹ and 20 km s⁻¹, respectively, while for AzV 377, we simply assumed v_mac = 0 km s⁻¹.

C2. The projected rotational velocity v sin i is determined by convolving the spectra during the fitting process. No explicit macroturbulence is assumed (v_mac = 0km s⁻¹).

P1. We employ the iacob-broad to measure the projected rotation velocity v_rot sin i, assuming no macroturbulence. As metal lines are, contrary to helium and hydrogen lines, not significantly affected by pressure broadening, we choose O IV at 3397 Å and of N IV at 3404 Å, 3412 Å, and 3414 Å which are all clearly visible in the optical spectrum of the O5 V((f)) star. The mean value obtained from the different lines is v_rot sin i = 60 km s⁻¹.

P2. We used a combined Fourier transform (FT) and goodness-of-fit (GOF) analysis employing the iacob-broad tool. We applied this method to several absorption lines such as Si IV λ 4089, He I / 4713, O III λ 5592, C IV λ 5801, 5812, and He I λ 5876. The overall mean v sin i and macro-turbulent (v_mac) velocities are given in Tables C.3 and C.2. Subsequently, these values, along with instrumental broadening, were used to convolve the model spectra to match the observations.

B.2 Temperature, surface gravity, and abundances

Depending on the method, the reddening is either determined before, after, or in parallel to other parameters.

F1. The IACOB-GBAT package (Simón-Díaz et al. 2011) is used to determine the stellar parameters from the optical H and He lines, without re-normalizing the normalized spectrum provided in the XSHOOTU data release (Sana et al. 2024). The package uses a grid of FASTWIND models to determine the effective temperature T_eff, the surface gravity log g, the He/H number ratio y_He, the wind strength parameter Q_ws, the ß exponent (see Eq. 1), and the microturbulence velocity ξ. The underlying grid for the LMC comprises of nearly 10⁵ model atmospheres at metallicity Z = 0.008 calculated with the FASTWIND version 10.5 using HTCondor⁹ . The grid covers the parameter ranges listed in Tab. 2. For the SMC, a dedicated grid centered around the parameters of AzV 377 at Z = 0.004 was used. IACOB-GBAT allows to automatically find the best-fitting model within the grid by looking for the minimum χ² obtained by adding the individual ${\chi }_i^2$ of each individual line, in such a way that lines with a larger uncertainty or photon noise have a smaller weight in the final ;χ². A full description of the procedure is given in Holgado etal. (2018, Appendix A)¹⁰.

F2/F3. The basic stellar parameters are determined together with the wind parameters, see Appendix B.4. To compare FASTWIND models to the data, the spectra need to be normalized. The UV spectra were normalized using a CMFGEN model from the grids of Bestenlehner et al. (2014) and Marcolino et al. (2024), with the method as described in Brands et al. (2022). For the fitting process of Sk -69° 50 and Sk -66° 171, the following optical lines, sorted by ion and wavelength, are considered: Hϵ, Hδ , Hγ, Hß, Hα, He I λ4026, He I λ4387, He I λ4471, He I λ4922, He I λ5875, He I λ7065, He II λ4200, He II λ4541, He II λ4686, He II λ5411, He II λ6406, He II λ6527, He II λ6683, C III λ4650, C III λ5695, C IV λ5801, N III λ4640, N IV λ3480, N IV λ4058, N IV λ6380, N V λ4603, N V λ4620, and O III λ5592. For the combined optical+UV fit, we additionally include: P V λλ1118-1128, C IV λ1169, C III λ1176, N V λ1240, O IV λ1340, and Si IV λλ1394-1402, and C IV λ1550 for both Sk -69° 50 and Sk -66° 171. For Sk -66° 171, we further include N IV λ1718 and N III λ1751. For AzV 377 the same lines are used as for Sk-66° 171, excluding He II λ6406, N IV λ4058, N IV λ6380, and the N V lines, but including the He I λ4143, He I λ4713, He I λ5015, He II λ1640, N III λ4515, and Si III λ1113.

F4. The spectroscopic analysis pipeline minimizes the χ² and returns a 4D (T_eff, log g, M, and y_He) probability distribution function. The 3 variations of CNO abundances provide only a rough estimate as they are not treated as independent parameters. The methodology is described in Bestenlehner (2022) and Bestenlehner et al. (2024). Uncertainties of the stellar parameters have been derived by defining 1 σ confidence intervals. The resulting uncertainties are typically larger compared to the other methods presented in this work due to correlations between all four stellar parameters.

Without weighting specific spectral lines, the ;χ² is dominated by Balmer lines, which are the strongest and broadest lines in early-type stars. Aiming to be able to eventually analyze the entire ULLYSES sample from mid-B dwarfs to early-O super-giants, we chose an approach that is not tailored to a specific parameter regime. In the wavelength range of the Balmer lines we remove every second wavelength point, thereby decreasing the wavelength resolution by a factor of two. In regions containing metal lines we add points to artificially increase the wavelength resolution by a factor of two. Helium lines are left untouched. In addition, we mask interstellar lines and bands. For mid- and late-O stars, the temperature determination is dominated by the ionization balance of He I and II. For early-O stars and for B stars, the nitrogen and silicon lines respectively contribute enough to the total χ² to obtain reliable temperatures. The surface gravity is derived on the basis of the Balmer line wings. The He/H number ratio y_He is determined from the H and He lines, except for Sk-69° 50 where grid issues prevented a determination of this quantity.

C1. T_eff and log g are constrained using a χ² analysis on the initial model grid, relying on the optical spectra only. For AzV 377, we re-normalized the spectrum by eye to optimize the use of the entire wavelength range, in particular close to the Balmer jump. For the LMC stars, we applied no re-normalization at this step. The χ² calculations is limited to a list of lines which includes all the hydrogen and helium lines as quoted above for models F2/F3, but we did not use the CNO-lines listed therein. We first determine T_eff and log g from using the grids of models. We stress that since He I lines are predominant in this list, we gave more weight to the He II in order to avoid a systematic bias of the results towards low T_eff. This weight was set to 2 in this work.

In a second step after the luminosity determination (see below), y_He and CNO abundances are determined from the analysis of UV and optical (normalized) spectra. Starting from a model with the determined T_eff, log g, we computed new sets of models varying the abundances and microturbulence ξ. We then apply the same type of χ² analysis as before, but this time for a set of He, C, N, or O lines. Essentially, we computed new atmosphere models for typically five to seven different abundances and used values of 5, 10, 15 and 20 km s⁻¹ for ξ (applied in the line absorption profiles in the CMF radiative transfer as well as in calculation of the synthetic spectra). We refer to Martins et al. (2024) for more details on the methodology as well as a discussion of the dispersion of the results on the abundances. This approach does not propagate uncertainties in the fundamental parameters into the errors on the abundance measurements (see also Martins et al. 2012a, Martins et al. 2015).

C2. Utilizing an earlier set of calculations previously used in Maryeva et al. (2014); Gvaramadze et al. (2018, 2019) we ran a series of custom CMFGEN models for the two LMC stars. The initial assignment of T_eff, L, and log g is based on the spectral types of the sample stars and the O-star calibration by Martins et al. (2005). We then proceed with a more detailed parameter determination. To measure T_eff, we used the Si III-Si IV, N II-N III-N IV, C III-C IV, and He I-He II lines. In the initial models for estimating T_eff and L (see below), we assume y_He = 0.11 for Sk-69° 50 and y_He = 0.14 for Sk-66° 171. We further assume increased abundances of nitrogen, and decreased carbon and oxygen abundances in these initial models, while for all other elements (Al, Si, P, S, Fe), subsolar abundances scaled to the metallicity of LMC were taken. After the determination of T_eff and L, we compared the locations of the stars in the H-R diagram to Geneva evolutionary tracks from Eggenberger et al. (2021) for a better estimate of H, He, C, N, and O. Finally, these abundances are then slightly adjusted to better match the observed spectra. The resulting models show a considerably deficiency in the C III 4647-50-51 feature, which could not be resolved without affecting the fit quality of other lines. While it is possible to decrease the intensity of this triplet by decreasing the carbon abundances, this would also decrease the otherwise sufficient intensities of the C III 5696 and C IV 5801,5812 lines.

P1. To estimate the effective temperature and the surface gravity of the star, we use the ratio between the He I and He II lines. The He I lines play a crucial role in setting the lower limit of the temperature. In our fitting procedure, we excluded He I lines with lower level 1s2p ¹P as these are known to be unreliable (Najarro et al. 2006). To determine the surface gravity of the star, the wings of Hß, Hγ, and Hδ are used, as they are almost insensitive to the wind. The best fit is achieved with an effective temperature of T_eff = 44.6 ± 2.0 kK and a surface gravity of log(g/cms⁻²) = 4.0 ± 0.1. Notably, temperature and surface gravity are linked to each other and obtaining the best fit is an iterative process. For the microturbulent velocity in the hydrostatic integration, a value of 10 km s⁻¹ is chosen, similar to what is used as a minimum in the formal integral. Compared to a model that does not include turbulence in the hydrostatic equation, we expect a higher log g of about 0.04 according to Eq.(4).

In the optical and UV spectra of AzV 377, one can see several metal lines associated with CNO elements, yielding valuable information on their surface abundances. For carbon, we use the C IV λ1169 and the C III λ1175 lines in the UV and the C IV λλ5801,5812 doublet in the optical. We find that these lines are best matched when lowering the carbon abundance to X_C = 15 × 10⁻⁵ ± 5 × 10⁻⁵. For nitrogen, we see strong absorption lines of N IV at 3404 Å, 3412 Å, and 3414 Å as well as weak emission of N IV at 4058 Å. These lines can only be matched when using an increased nitrogen abundance of X_N = 80 × 10⁻⁵ ± 15 × 10⁻⁵. To determine the oxygen abundance, the O IV λλ1339, 1343 doublet and the O III λλ1409,1410,1411,1412 multiplet in the UV and the O IV λλ3381,3386,3390,3397,3410 multiplet in the optical are used, yielding $X_{\text{O}}={50}_{-10}^{+20}\times {10}^{-5}$.

P2. The stellar temperature is primarily constrained by the helium ionization balance. Since He I singlet lines are more susceptible to model details (Najarro et al. 2006), we give them less weight during the fitting process. In addition, we used the Si III-Si IV, C III-C IV, and N III-N IV line ratios to determine the temperature. After getting a constraint on the temperature, we measured the surface gravity using the pressure-broadened wings of the Balmer lines. The main diagnostic lines are Hδ, Hγ, and Hß, since they are less affected by wind emission. For Sk-66° 171, the combined optical and UV spectra are best reproduced using a model with T_eff = 29.3 kK and log g = 3.08. The N III emission lines in the optical mark a notable exception, where the temperature is too low to replicate them. On the other hand, higher temperature models reproducing the N III emission lines would spoil the overall spectral fit, especially the strength of He I and He II lines. Moreover, the C IV λ1169 and C III λ1176 lines in the UV are temperature-sensitive, and we combine them with other optical diagnostics to obtain the best fit. The combined UV and optical spectra of Sk -69° 50 are found to be best reproduced by T_eff = 33.4 kK and log g = 3.42. The microturbulent velocity of ξ = 14km/s⁻¹ included in the hydrostatic equation leads to an increased log g of about 0.1 dex following Eq. (4). The microturbulence entering the hydrostatic equation is slightly lower than our eventually preferred value for the minimum microturbulence in the formal integral (18 and 20km s⁻¹).

We adjusted the CNO abundances in the models to reproduce the observed strength of their respective absorption lines. For Sk-69° 50, we increased the nitrogen mass fraction in the model by a factor of 60 times the baseline LMC values. The best-fit model reproduces N IV and N III absorption lines in the optical and UV ranges while under-predicting the N III emission. To match the strength of the N III λ4097 absorptio-n line in Sk-66° 171, we increased the nitrogen mass fraction in the model by a factor of ten. The He line strength for both stars is found to be better predicted by models with a modestly higher He enrichment. For the remaining elements we adopted typical LMC abundance values derived from OB stars (Trundle et al. 2007) if available, or were otherwise adopted as half-solar abundances.

B.3 Luminosity and reddening

F1. We derive the stellar radius and the resulting luminosity by comparing the model magnitudes from the best fit with the B, V, J, H and K observed magnitudes listed in Tab. 1. To this end, we also need to determine the extinction. We have used the extinction law by Maíz Apellàniz et al. (2014) and characterized it by the R₅₄₉₅ index. Ideally, all photometric bands will give the same radius for the correct extinction. In practice, we take the extinction where the dispersion for the different radii values is lowest. We note that this dispersion is the major contribution to the error budget for the stellar radius. We get a best value of R5495 = 3.8 and 4.0 for Sk-69° 50 and Sk-66° 171, whereas for AzV 377 the best value is at R₅₄₉₅ = 2.5. This points to a different extinction nature for both cases, although the differences with the canonical value R₅₄₉₅ = 3.1 are small. For AzV 377, the K-band was not included in the luminosity determination due to an inaccuracy in the K-magnitude distributed among the team for this star which was only corrected after the F1 calculations were already completed.

F2/F3. To derive the luminosity, Kiwi-GA requires an anchor magnitude, for which we use the K_s-band motivated by its generally low extinction. We derive the absolute K_s-magnitude by fitting the Fitzpatrick (1999) extinction law through the U, B, V, J, H, and K_s magnitudes listed in Vink et al. (2023), adopting the CMFGEN model employed for the normalization as the intrinsic model. With the obtained $A_{K_{\text{s}}}$ value and the adopted LMC distance given in Sect. 2, we derive $M_{K_{\text{s}}}=-5.92\text{\hspace{0.17em}mag}$ (Sk -66° 171), $M_{K_{\text{s}}}=-4.89\text{\hspace{0.17em}mag}$ (Sk-69° 50), and $M_{K_{\text{s}}}=-3.67\text{\hspace{0.17em}mag}$ (AzV 377). For F2 and F3, the flux-calibrated UV spectra do not enter the luminosity determinations.

F4. The photometric data from Table 1 is utilized to determine the reddening and the luminosity L. We fit the model SED with the limfit routine¹¹ to the optical to near-IR photometry by applying the extinction law by Maíz Apellániz et al. (2014), similar to the procedure performed in Bestenlehner et al. (2022). The uncertainties from the SED fit are used to estimate the uncertainties on the photometric data. Together with the uncertainties on the stellar parameters, these are propagated into the error of L.

C1. Once T_eff, log g have been derived from the analysis of the optical spectrum, the stellar luminosity is constrained by comparing the SED of the corresponding model to the “observational” SED built from the flux-calibrated FUSE + HST spectra in the UV, as well as optical and NIR photometry. For the latter, we employed the photometric data listed in Table 1 and converted the magnitudes to flux, using the SYNPHOT Python packages available on ASTROPY. The bolometric corrections were derived from Martins & Plez (2006) for the previously derived effective temperature. In the determination of the luminosity, the interstellar extinction has to be derived as a side-product. For this, the synthetic spectra were reddened, using the Fitzpatrick et al. (2019) extinction law for the Milky Way foreground reddening, and Gordon et al. (2003) for the internal reddening in the Magellanic Clouds. The distances given in Sect. 2 are used to scale the synthetic flux for the sample stars, enabling proper comparisons to the observations.

C2. At first, the luminosity L is roughly estimated by comparing the spectral energy distribution (SED) of the model spectrum with the photometric measurements over the whole spectral range. Beside the photometric data given in Table 1, we also used data from the XMM-Newton Serendipitous Ultraviolet Source Survey (XMM-SUSS; Page et al. 2012), Spitzer SAGE infrared photometry (Bonanos et al. 2009), and the Wide-field Infrared Survey Explorer (WISE) (Cutri et al. 2021). For a more accurate determination of the luminosities, the synthetic magnitudes of the stars are then calculated in the U, B, and V filters and compared with observations. In order to obtain the magnitudes for the model spectra, we first rescale the fluxes to the distance of the LMC. The resulting fluxes are then corrected for the interstellar extinction with the IDL procedure FM-UNRED. The LMC2 option is selected, employing average LMC reddening parameters from Misselt et al. (1999). Afterwards, the calculated spectra are convolved with the transmission curves of the standard U, B, and V filters, and the corresponding zero points are applied using the PYSYSP Python package (Casagrande & VandenBerg 2014).

P1. The luminosity is constrained by fitting the synthetic model flux to the photometry from Table 1, yielding a luminosity of log(L/L_⊙) = 5.3 ± 0.1. The reddening is modeled as a combination of the color excess arising from the Galactic foreground with E_b–v,Gal = 0.03 mag and the local SMC reddening with E_b–v,SMC = 0.028 mag, using the reddening laws by Seaton (1979) and Bouchet et al. (1985), respectively.

P2. The luminosity and color excess E_B–V are derived by fitting the model SED to the photometry from Table 1 and the flux-calibrated UV spectra from ULLYSES. We provide the total extinction including the contribution from the Galactic foreground (assuming E_B–V = 0.04 mag), adopting the reddening law from Seaton (1979) for the Milky Way part, and the reddening law described in Howarth (1983) with R_V = 3.2 for the LMC contribution.

B.4 Wind parameters, including clumping and additional X-rays

Some of the methods separate the determination of the stellar parameters from the wind parameters, while others vary both at the same time.

F1. The optical spectrum for the studied targets does not allow to derive the terminal wind velocities. We thus rely on Eq. (5) from Hawcroft et al. (2024) for AzV 377 – which consequently has a large error – and the terminal velocity values by Brands et al. (in prep., part of the XShootU series) for Sk-69° 50 and Sk-66° 171. To determine the mass-loss rates, Eq. 12 is evaluated. For the dwarf star AzV 377 only an upper limit for Q_ws and thus Ṁ can be derived, while the two LMC stars with stronger winds enabled to constrain Q_ws and thus Ṁ (with the additional information abut v_∞). All models are unclumped, and no additional X-rays were added to any of the models.

F2/F3. Once the input for a Kiwi-GA run is prepared, the normalized spectrum, the luminosity anchor, and the line selection are fed to the algorithm, and the best fit including uncertainty values are derived. We do this twice for each star: once with only the optical line selection, and once with the combined optical and UV line selection. For the optical-only run, we have to make assumptions about the wind clumping due to the lack of more diagnostics. We adopt fixed “macroclumping” parameters, with clumping factor f_cl = 10, inter-clump density contrast f_ic = 0.1, velocity-porosity f_ic = 0.5, and onset velocity of clumping v_cl,start = 0.05_∞ (for details about these clumping parameters, see Appendix A.2.1 and Sundqvist & Puls 2018). The terminal velocity is also fixed in F2, to a value estimated from the position of the blue edge of C IV λλ1548-52. We further fix the wind acceleration parameter ß = 1.0, and the micro turbulence velocity to 15 km s⁻¹. The parameters we fit are: T_eff, log g, v_rot sin i, He, C, N and O abundance, and the mass-loss rate Ṁ. For the optical + UV run, we adopt the He abundance and the v_rot sin i to the best fit values from the optical-only run, and then fit the following free parameters: T_eff, log g, C, N, and O abundance, Ṁ, ß, terminal velocity v_∞, wind turbulence velocity v_windturb, f_cl, inter-clump density contrast f_ic, velocity filling factor f_vel, and clumping onset velocity v_cl. We assume an LMC metallicity of Z = 0.5Z_⊙. Furthermore, we include X-rays following the prescription of Carneiro et al. (2016). The X-rays are described by a radius-dependent shock temperature, $$T_{\text{S}}\left(r\right)=\frac{3}{16}\frac{\mu m_{\text{H}}}{k_{\text{B}}}u{\left(r\right)}^2,$$(B.1)

with a jump velocity, $$u(r)=u_{\infty }{\left[\frac{v(r)}{v_{\infty }}\right]}^{{\gamma }_{\text{x}}},$$(B.2)

that has two free parameters u_∞ and γ_x. We fix γ_x = 0.75 (Brands et al. 2022) and assume u_∞ = 0.3v_∞ for the maximum jump velocity u_∞. The necessary X-ray volume filling fraction is chosen such that the X-ray output luminosity is approximately 10⁻⁷ times the stellar luminosity (for details, see Brands et al., in prep., part of the XShootU series).

F4. As we apply this method only to the optical spectra, no value for v_∞ can be measured. As mentioned for F1, Q_ws can be estimated from the emission line profiles of Hα, and, in the case of mid-early O stars, also from He II λ4686, but v_∞ needs to be known to derive the mass-loss rate Ṁ. For F4, no assumptions for u₀ were made and thus no mass-loss rates were constrained. For AzV 377, only an upper limit for Q_ws can be determined (cf. Table C.1). For the two LMC targets, we find best fits using models with Q_ws = −12.3. All models are unclumped (f_cl = 1) and no additional X-rays were added to any of the models.

C1. The initial grid models do neither include clumping nor X-rays, and their wind mass-loss rates are computed using the Vink et al. (2001) formula for the corresponding metallicity, adopted to be 0.2 Z_⊙ and 0.5 Z_⊙ for SMC and LMC, respectively. Moreover, v_∞ = 3.0 v_esc,eff is assumed to be consistent with Martins & Palacios (2021) and references therein. The wind acceleration parameter is further fixed to ß = 1.0 For the LMC stars, we then refined the wind parameters by calculating dedicated models beyond the initial assumptions of the LMC grid. The wind terminal velocities were determined from the UV PCygni profiles (C IV resonance profiles are saturated for both Sk-69° 50, and Sk-66° 171). All the models computed later on adopted the measured values of v_∞. The mass-loss rate, the ß exponent of the velocity law, and the volume filling factor for the wind clumping were constrained from the UV and optical spectra. For the models with clumping, we adopted v_cl = 30 km s⁻¹ (see, e.g., Bouret et al. 2005). We started with log L_X/L_bol = -7 but further tuned this quantity to improve the strength of N IV λλ1238–1242. For the SMC star AzV 377, no wind parameters were determined.

C2. The terminal velocity v_∞ is estimated using spectral lines with P Cygni profiles in the UV. Intensities of emission lines were used for refining the mass-loss rate Ṁ. For the clumping parameters, we adopted f_V,∞ =0.2 and 0.3 as well as v_cl=20 km s⁻¹ and 10 km s⁻¹ for Sk-69° 50 and Sk-66° 171, respectively. We further included additional X-rays in the models for both stars to reproduce the N v λ1239,1243 doublet. In a last step, the microturbulent velocity ξ was also adjusted to better match the spectra.

P1. For AzV 377 no optical wind diagnostics are available. In the UV, N V λλ1239,1243 doublet has a strong pronounced P Cygni profile, while the C IV λλ1548,1551 doublet only shows a weak absorption trough. After measuring the terminal velocity from these lines, we find that this odd morphology can only be explained with a mass-loss rate of log (Ṁ[M_⊙ yr⁻¹]) = −7.8(1). The density contrast D is assumed to be 10, which corresponds to the clumping factor fcl in this case as we assume a void inter-clump medium. No clumping is assumed below the sonic point. No additional X-rays are included in the models as AzV 377 has a high enough temperature to intrinsically show N V.

P2. After fixing T_* and log g_* (see above), we adjust the wind parameters and recalculate the models accordingly. We measure the terminal velocities (v_∞) from the blue edge of the absorption trough of the C IV λλ1548−1551 P-Cygni line. The ß parameter of the velocity law is varied such that the synthetic spectrum can reproduce the profile shapes of the UV resonance lines and the optical wind emission lines, such as Ha and He II λ4686. To get the best fit, we had to assume a double-ß law. For Sk-69° 50 we use ß₁ = 0.8 and ß = 1.5 with q_2ß = 0.80. In the case of Sk-66° 171, higher ß values (ß₁ = 1.1 and ß = 1.8 with q_2ß = 0.85) provide a better fit to the observed spectra. By consistently fitting the wind lines in the UV and optical, we infer the mass-loss rate and the clumping parameters. The primary diagnostic lines used are the UV resonance doublets such as C IV λλ1548–1551, Si IV λλ1393–1403, and N V λλ1238–1242 along with Hα emission. Since the models for both stars are found to be too cool to reproduce the N V λλ1238–1242 line profile, we incorporated an additional X-ray field in the models (Baum et al. 1992). For Sk -69° 50, we used T_x = 1 MK and an X-ray onset radius of R_x,min = 1.1 R_*, while for Sk-66° 171, we used T_X = 3 MK and R_X,min = 1.2 R_*. The X-ray filling factor was set to 0.1 for both stars.

Appendix C Best-fit result overview and atomic data coverage

C.1 Overview of best-fit values

In Tables C.1, C.2, and C.3, we provide an overview of the best-fitting parameters determined with each of the methods. Abundances are given in different formats for easier comparison. Depending on the method, not all parameters were determined. Fixed or assumed values are denoted in italics.

C.2 Wide-range spectral comparison

In addition to the discussions about individual lines in Sects. 4.1 and 4.2, we provide comparison plots in a large wavelength range in Figs. C.1, C.2, and C.3. Each of these figures shows the normalized observations compared to all the model spectra. For clarity, we show only one normalization of the observations, while some methods performed individual (re-)normalization of the data in some wavelength regimes. This is especially apparent in the UV range, where normalization is tricky due to the forest of iron lines and commonly done with the help of a continuum from an atmosphere model itself. As the resulting normalization is then model-dependent in such cases, this can lead to an apparent disagreement between the normalized observation shown in our figures and in the continuum level of methods using a (slightly) different normalization.

C.3 Ion coverage

Table C.4 provides an overview of all the ions taken into account by each of the methods. The handling of the atomic data varies considerably between the different model atmosphere codes. Below, we give information for each of the codes (denoted by their associated methods) to better interpret the values given Table C.4.

F1/F2/F3/F4. The number of levels and lines provided in Table C.4 for the methods F1 to F4 refer to the explicit and background elements in FASTWIND (see Sect. 3.2) as used in the specific models. Model F1 has been calculated with only H and He as explicit elements. The models for F2 and F3 treat H, He, C, N, O, Si, P as explicit elements, and the models for F4 do this for H, He, C, N, O, Si.

Most of the denoted levels are packed with respect to fine-structure, and are depacked when calculating the radiative transfer (for details, see Puls et al. 2005, 2020). In addition to the displayed levels and lines, transitions from other elements are considered as well, covering almost all atoms from H to Zn (except those with a negligible abundance). Moreover, the complete line list used for calculating the SEDs comprises not only the quoted explicit lines (when depacked), but also the multitude of additional lines with the ground-state or a metastable level as the lower level and a high-lying upper-level, where the source function is approximated by a two-level atom approach.

C1/C2. In CMFGEN, there is a distinction between the levels taken explicitly into account for the non-LTE solution of the statistical equations and the levels in the radiative transfer.

To keep the number of statistical equations manageable, levels are groups into superlevels with only each superlevel entering the rate equations. In the radiative transfer however, the levels are instead accounted for explicitly (i.e., no opacity sampling or redistribution). The grouping into superlevels can be changed by the user. Consequently, Table C.4 lists two columns for the super-levels and the “total levels” which refer to the number of levels accounted for in the radiative transfer. A more in-depth description of the superlevel approach is provided in Hillier & Miller (1998).

P1/P2. POWR differentiates between “normal” and iron-group elements. For normal elements, recommended sets of atomic data are available where the fine-structure of the higher levels is usually packed into single levels. The user then usually just requests which elements and ions should be taken into account. Each of these levels is then treated in full non-LTE, i.e., explicitly enters the rate equations and their transition are taken into account for the radiative transfer (unless they are radiatively forbidden). As this can lead to very large sets of equations, the user further has the option to limit the total number of levels for a given ion.

For iron and other elements of the “iron group” (Sc to Ni), a different approach is used in order to ensure that full blanketing effect is taken into account while also the forest of UV lines can be reproduced sufficiently. Grouped by energy bands, and in more recent years also by even and odd parity, the large list of iron levels is reduced to a set of superlevels. The allowed transitions between the superlevels are grouped into “superlines” for which a wavelength-dependent cross-section is calculated. To account for different Doppler velocities (and thus profile broadenings) in the radiative transfer, different pre-calculated sets are available. In the actual model calculation, only the population numbers for the superlevels are calculated, taking into account their different nature compared to the “normal” levels. In the radiative transfer, the pre-calculated cross-sections ensure that all transitions are taken into account at their correct wavelengths. The POWR superlevel concept (without parity splitting) is described in more detail in Gräfener et al. (2002).

Fig. C.1 Comparison of the different model spectra (colored lines, see legend) with the normalized observations of AzV 377 (thick black line) across a wide range of UV and optical wavelengths.

Fig. C.2 Comparison of the different model spectra (colored lines, see legend) with the normalized observations of Sk –69° 50 (thick black line) across a wide range of UV and optical wavelengths.

Fig. C.3 Comparison of the different model spectra (colored lines, see legend) with the normalized observations of Sk -66° 171 (thick black line) across a wide range of UV and optical wavelengths.

References

All Tables

All Figures

	Fig. 1 Comparison between the different methods in the Balmer lines and He λ14541 region for AzV 377.
In the text

	Fig. 2 Comparison between the different methods in the Balmer lines and He II λ4541 region for Sk-66° 171.
In the text

	Fig. 3 Comparison between the different methods in the Balmer lines and He II λ4541 region for Sk-69° 50.
In the text

Fig. 4 Comparison between the main UV profiles for AzV 377. The panels from (a) to (h) depict respectively the profiles of O VI 1032/1038 Å, S IV 1063/1073Å (and He II 1085 Å), PV 1118/1128 Å, C III 1176 Å, NV 1238/1242 Å, Si IV 1394/1403 Å, C IV 1548/1551 Å, and N IV 1718 Å. Interstellar Lya and Lyß absorption affects some of the diagnostics; most notably the wing of NV 1238/1242 Å in case of larger terminal velocities (e.g., seen in AzV 377).

In the text

	Fig. 5 Comparison between the main UV profiles for Sk-66° 171. The spectral windows are the same as in Fig. 4, following the same order.
In the text

	Fig. 6 Comparison between the main UV profiles for Sk-69° 50. The spectral windows are the same as in Fig 4, following the same order.
In the text

	Fig. 7 Comparison between the lines with more prominent disagreement on the optical for AzV 377.
In the text

	Fig. 8 Comparison between the lines with more prominent disagreement on the optical for Sk -69° 50.
In the text

	Fig. 9 Comparison between the lines with more prominent disagreement on the optical for Sk-66° 171.
In the text

	Fig. 10 Modeled and observed SED of the targets. The flux points correspond to the magnitudes listed in Table 1, following the same order (from bluer to redder). The observed flux-calibrated UV spectra correspond to those acquired by ULLYSES.
In the text

	Fig. 11 Zoom-in comparison of the SEDs from the different models for AzV 377 around the applied photometry (small upper panels, crosses mark photometric measurements) and the flux-calibrated UV spectra (big lower panel).
In the text

	Fig. 12 Same as Fig. 11, but for Sk-66° 171.
In the text

	Fig. 13 Same as Fig. 1, but for Sk-69° 50.
In the text

	Fig. 14 HRD with the obtained positions for the two LMC stars Sk -69° 50 and Sk -66° 171 – indicated by the respective labels. For comparison, tracks from Brott et al. (2011, up to Minit = 50 M⊙) and Köhler et al. (2015, from 60 M⊙) are shown.
In the text

	Fig. 15 HRD with the obtained positions for the SMC star AzV 377. For comparison, tracks from Brott et al. (2011, up to Minit = 60 M⊙) and Köhler et al. (2015, from 80 M⊙) are shown.
In the text

	Fig. 16 Ratio of the determined spectroscopic masses to evolutionary masses (based on Brott et al. 2011 or Köhler et al. 2015 respectively) for the different methods and sample stars.
In the text

	Fig. 17 Comparison of the derived (F1, F2: assumed) terminal wind velocity v∞ versus the derived effective temperature. The dashed lines denote the LMC (gray) and SMC (salmon) relations from Hawcroft et al. (2024).
In the text

	Fig. 18 Ratio between the terminal wind velocity v∞ and the effective escape velocity vesc,Γ as a function of Teff. The dashed gray line shows the empirical results from Kudritzki & Puis (2000).
In the text

	Fig. 19 Mass-loss rate $\dot{M}$ versus stellar luminosity L for our sample stars analyzed with the different methods. For comparison, we also plot the mass-loss recipes from Vink et al. (2001, solid shading), Krtička & Kubát (2018, dashed lines), and Björklund et al. (2021, hatched shading). The SMC comparison data are drawn in light red.
In the text

	Fig. 20 Modified wind momentum rate versus stellar luminosity L for our sample stars analyzed with the different methods.
In the text

	Fig. 21 Transformed mass-loss rates versus the ratio between luminosity L and spectroscopic mass M for our sample stars analyzed with the different methods.
In the text

	Fig. A.1 Comparison of the different radial behaviors of the clumping factor fcl = D, assuming optically thin clumping and input values used for the analysis of Sk-66° 171. To focus on the clumping stratification, all curves assume the same maximum clumping factor (fcl = 20) and the same underlying velocity field (v∞ = 1800 km s−1 ).
In the text

	Fig. C.1 Comparison of the different model spectra (colored lines, see legend) with the normalized observations of AzV 377 (thick black line) across a wide range of UV and optical wavelengths.
In the text

	Fig. C.2 Comparison of the different model spectra (colored lines, see legend) with the normalized observations of Sk –69° 50 (thick black line) across a wide range of UV and optical wavelengths.
In the text

	Fig. C.3 Comparison of the different model spectra (colored lines, see legend) with the normalized observations of Sk -66° 171 (thick black line) across a wide range of UV and optical wavelengths.
In the text