© The Authors 2024

1 Introduction

Massive stars (M > 8 M_⊙) are born as hot objects (earlier than spectral type B2 V; e.g., Harmanec 1988; Hillier 2020) and are far outnumbered by low-mass stars. However, they deeply impact their surroundings and shape the galaxies’ dynamical and chemical evolution, both due to their usually violent deaths and to their powerful winds, energetic radiation largely emitted in the ultraviolet (UV), and the metal-enriched material they depose in their vicinity. Thus, understanding how these objects behave throughout cosmic history is paramount to comprehending many aspects of the Universe since the birth of the first stars.

Among the hot stars, B-type stars, especially supergiants, contribute significantly to that feedback in stellar populations younger than 50 Myr (de Mello et al. 2000). Additionally, most B stars observed in external galaxies are B supergiants (BSGs)¹, due to their usually high luminosity (e.g., Kudritzki et al. 2012). Recently, this has reached even cosmological dimensions, as lensed BSG candidates were inferred from spectra taken of high-redshift objects (e.g., z ~ 4.8 and z ~ 6.2) using JWST (see, e.g., Welch et al. 2022; Furtak et al. 2024). BSGs can even be progenitors of supernovae (e.g., SN 1987A, Walborn et al. 1987), resulting in compact objects such as neutron stars (NSs) and black holes (BHs). Their properties and evolution are therefore a major puzzle piece toward understanding high-mass stellar evolution and its endpoints, which are crucial in order to properly interpret gravitational wave events produced by BH and NS mergers (e.g., Abbott et al. 2016, 2017).

B supergiants play a crucial role in our understanding of line-driven wind physics and mark an important stage among the puzzling diverse evolutionary paths of massive stars. BSGs lie at the cool end of the hot stars regime, with effective temperatures ranging from ~30 kK down to ~10 kK. Within this parameter range, important and drastic changes in the wind properties (e.g., terminal wind velocities, clumping, and X-rays) appear to happen, especially at spectral types around B1, or temperatures of ~22 kK (e.g., Lamers et al. 1995; Driessen et al. 2019; Berghoefer et al. 1997; Petrov et al. 2014). Traditionally, these changes are attributed to the “bi-stability phenomenon” (e.g., Pauldrach & Puls 1990; Lamers et al. 1995; Vink et al. 1999), whose understanding is paramount to constraining the feedback and the evolution of these stars properly.

From a more observational perspective, such changes include a relatively sharp variation in terminal wind speed (e.g., Lamers et al. 1995; Markova & Puls 2008) and in the X-ray luminosity (Berghoefer et al. 1997) within this temperature range. Such drastic changes align with theoretical findings that massive stars might also experience variations in their wind structure (e.g., Driessen et al. 2019), ionization of elements (Vink et al. 1999; Petrov et al. 2014, 2016), and a steep jump in mass-loss rates (Pauldrach & Puls 1990; Vink et al. 1999, 2000). The last is particularly relevant for the evolution of high-mass stars, as different mass-loss rates can significantly change a star’s evolutionary path (see, e.g., Vink et al. 2010; Renzo et al. 2017). However recent theoretical works challenge the existence of a steep jump in the mass-loss rates in that region, either predicting a mild increase (Krtička et al. 2021) or even predicting a continuous decline (Björklund et al. 2023).

At solar metallicity, recent empirical studies point toward an increase in mass-loss rates in the bi-stability region (e.g., Bernini-Peron et al. 2023), albeit less pronounced than predicted by Vink et al. (1999). However, it is unclear whether and how the bi-stability phenomenon would manifest at lower metallicities. For instance, Evans et al. (2004b) shows that terminal wind velocities (v_∞) of Small Magellanic Cloud (SMC) BSGs are similar to those from Galactic BSGs, on the hot side and and on the cool side of the jump region, which resonates with theoretical predictions finding weak (Leitherer et al. 1992) to almost no (Vink & Sander 2021) dependence of v_∞(Z), in particular for cooler BSGs. When looking at mass-loss rates ($\dot{M}$), most studies in the literature invoke a power-law scaling of $\dot{M}$ with (Z/Z_⊙) with the most commonly employed scaling being $\dot{M}$ ∝ Z^0.70 (Vink et al. 2001). However, the actual $\dot{M}$(Z) dependence could vary significantly depending on the temperature, luminosity, and proximity to the Eddington limit (see, e.g., Marcolino et al. 2022; Krtička et al. 2024, for recent studies).

For this work we analyzed a representative sample of 18 SMC BSGs across the bi-stability jump region at ~22 kK. To determine their stellar and wind properties, we used data from the ULLYSES initiative and the X-Shooting ULYSSES (XShootU) collaboration plus previous X-Shooter data from the ESO archive. We focused particularly on deriving the stars’ terminal wind velocities, mass-loss rates, clumping factors, and X-ray properties. Additionally, with the stellar properties at hand, we discuss the evolutionary status of the BSGs. Since the stars are in the SMC, the distance is well-constrained which is advantageous for obtaining luminosities and masses.

For our analysis we used the non-LTE comoving-frame stellar atmosphere code CMFGEN (Hillier & Miller 1998) to obtain the stellar and wind properties. A similar approach was previously taken by Evans et al. (2004a), albeit for a smaller sample of objects and only up to early-type BSGs. In addition to extending the sample, in particular to the crucial regime of later-type BSGs, we also take into account X-rays in this regime for the first time. The outcomes of our study provide insights into the general behavior of BSGs at low metallicity (e.g., relations between wind and photospheric properties). Additionally, it will reduce the necessary parameter space for follow-up studies that will address the full ULLYSES BSGs sample. Specifically, we provide constraints on the relative X-ray luminosity (L_X/L), clumping, and evolutionary properties.

In Sect. 2 we discuss the selection of our targets, the XShootU data acquisition, the selection of additional archival data, and interesting features we noticed in the spectra. We explain our analysis methods in Sect. 3 before discussing the results in Sect. 4 (stellar properties), Sect. 5 (X-rays), and Sect. 6 (wind properties). The conclusions are drawn in Sect. 7.

2 Observational data

2.1 Sample selection

In the recent final ULLYSES data release DR7² there are 30 SMC stars classified as BSGs. For this work, we performed a pilot study of the different spectral subtypes. Thus, we selected, whenever possible, two stars per spectral subtype as representatives of the full sample. In the target selection process, the literature spectral classifications from Paper I (Vink et al. 2023) were used, some of which have been revised (cf. Table 1) following classification of new VLT/Xshooter observations following the SMC B supergiant scheme of Lennon (1997) with luminosity classes assigned using H-gamma morphologies from Galactic B templates of Negueruela et al. (2024). In total, we cover 18 of the 30 ULLYSES targets. A quick inspection of the optical and UV spectra indicated that some of our targets have similar spectral types and a similar spectral appearance in the optical while showing extreme differences in the UV. This hints at important differences in wind properties despite similar photospheric properties, underlining the importance of picking more than one target per subtype where possible. We discuss our findings on this matter in Sect. 4.

In this work, we focus only on apparently single stars. Therefore, we exclude objects such as the known high-mass X-ray binary AzV 490 (e.g., Dickey et al. 2023) or Be/B[e] stars (namely, AzV 261 and AzV 16). While this does not rule out the presence of so far undetected companions, we checked the existing, partially multi-epoch, spectral data for signatures of multiplicity (e.g., radial velocity variations, steps in the UV absorption troughs, additional broad components in absorption lines) and found no clear evidence for them in any of the targets. The full list of our selected targets is shown in Table 1.

Fig. 1 Positions of BSGs in the SMC from ULLYSES. The filled diamond symbols indicate the BSGs analyzed in this study, whereas the empty diamonds indicate the remaining BSGs of the total ULLYSES sample.

2.2 Sample spectra

In Fig. 1, we plot the positions of the ULLYSES BSGs in the SMC. We assume that all the stars are to be located at a distance of 62.45±2.61 kpc (Graczyk et al. 2020) and correct the spectra for the radial velocities of each star in our analysis.

For the spectral analysis, we use HST UV spectra from the ULLYSES dataset and optical spectra from the XShootU collaboration (Sana et al. 2024) – see, e.g., Fig. 2. The UV spectra were obtained with FUSE, from 950 Å to 1190 Å, and COS plus STIS instruments for the rest of the UV range. FUSE operated using the LWRS grating offering a resolving power R = Δλ/λ ~ 10 000. COS operated in medium resolution modes, yielding R ~ 13 000. STIS E140M and E230M modes reach R ≳ 45 000 and R ≳ 30 000 respectively. The X-Shooter spectra were taken using the UVB and VIS arms, which cover respectively the range of 3100–5600 Å with R ~ 8000 and 5600 Å to 10 240 Å with R ~ 11 000. In Table A.2, we list for each of the stars we analyzed the instruments that were used together with their respective gratings and the spectra acquisition dates for the UV and optical regions, which for most of BSGs were not obtained simultaneously.

For targets that were not included in the XShootU data release (DR3), we retrieved the spectra from the ESO-Archive, also taken with the X-Shooter spectrograph, and normalized them by fitting the continuum manually. Some targets in our sample show peculiar line profiles in the optical and/or UV region. Such behavior may highlight potential binaries and/or atypical phenomena happening in the BSGs photospheres. To inspect eventual variability in spectral lines of such BSGs, we further used additional archival spectra retrieved from the ESO archive and MAST archive, which were subsequently normalized using the same technique. Information on these spectra can also be found in Table A.2. We discuss the results of our brief variability investigation in Sect. 4.

Fig. 2 UV and optical spectra used to derive stellar parameters of the sample stars. This figure illustrates the variation in some spectral features with the spectral type. From left to right, the columns show spectral regions around Si IV λ1400, C IV λ1550, O+N+Hδ+He I+Si lines, and Hα+C II.

2.3 Photometry

The UBV magnitudes were retrieved mainly from sources listed in (Vink et al. 2023), namely: Massey (2002), Ardeberg & Maurice (1977), Ardeberg (1980), and Azzopardi et al. (1975). Additionally, we included Gaia DR3 (Gaia Collaboration 2022) magnitudes (Gbp, G, and Grp). The near-infrared JHK magnitudes from 2MASS, also retrieved from Vink et al. (2023), are sourced from Cutri et al. (2003) and Cioni et al. (2011). The only exception is AzV 78, whose JHK magnitudes were collected from Bonanos et al. (2010). The far-infrared Spitzer and WISE magnitudes were collected from Bonanos et al. (2010) and the AllWISE catalog (Cutri et al. 2013) respectively. Each magnitude value can be found in the Appendix F available at Zenodo³.

3 Atmospheric modeling

For the detailed analysis of our sample stars, we use CMFGEN (Hillier & Miller 1998; Hillier et al. 2003), an open-source, non-LTE stellar atmosphere code that solves simultaneously and consistently the radiative transfer and the ionization/excitation population equilibrium equations. The program requires an input velocity structure and atomic data, which encodes the information about the transition levels and energies. Our atmosphere models included H, He, C, N, O, Ne, Al, Mg, Si, S, P, Ca, Cr, Ma, Fe and Ni, adding up to a total of more than 5 × 10⁵ individual atomic transitions – see Hillier et al. (2003) for details on the atomic data sources. Our initial models were calculated with a metallicity of 0.2 Z_⊙ using Asplund et al. (2009) for the elements listed in Vink et al. (2023). There is evidence for non-solar relative ratios of individual elements for the SMC, especially carbon (C), nitrogen (N), and oxygen (O). However, as BSGs in general present alterations in their CNO surface compositions (e.g., Trundle & Lennon 2005; Crowther et al. 2006), and we individually determine the abundances of each of the CNO elements, their initial baseline values were unimportant.

The procedure to obtain the stellar properties is similar to Bernini-Peron et al. (2023), albeit we explore more parameters in this work. Below, we summarize the procedure to obtain the stellar properties, highlighting the differences to Bernini-Peron et al. (2023). Except for the stellar luminosity, which is determined by the overall spectral energy distribution (SED), all the properties were obtained via fitting specific diagnostic lines, whose strength and shape respond well to certain stellar properties.

Given the high dimensionality, the non-linearity, the model computation times, and the intricate cross-dependence between the parameters in the atmosphere models, a fully automatic approach exploring a high-dimensional grid to obtain the stellar properties and their uncertainties is computationally unfeasible. Instead, we pick a suitable model from an existing set of models and iteratively adjust the stellar parameters to find the best compromise to reproduce the various diagnostics, as common in the analysis of hot stars with computation-intensive atmosphere codes such as CMFGEN and PoWR (e.g., Crowther et al. 2006; Martins et al. 2015; Ramachandran et al. 2019). Knowing the reactions of the lines to changes in the parameters from a sequence of model calculations, the quality of the spectral fitting is inferred from visually inspecting the different solutions (“by-eye fit”). This is a standard procedure in quantitative spectral analysis that has been successfully applied since the availability of complex non-LTE atmosphere models (e.g., Hillier & Miller 1999; Trundle & Lennon 2005; Martins et al. 2015; Ramachandran et al. 2019). As shown in the discussion of the different approaches and analysis methods in XShootU Paper IV (Sander et al. 2024), this type of analysis yields in general coherent and similar results compared to more automatized approaches.

To establish the uncertainties of fundamental stellar and wind parameters (e.g., T_eff, log g, abundances, $\dot{M}$, v_∞), we test variations of the individual parameters deviating from the best fit we could achieve. Given the inherent systematic uncertainties within the atmosphere modelling and the scatter between different atmosphere codes (cf. Sander et al. 2024), we refrain from doing this individually for each star, but perform our error margin study for selected individual hot and cool BSGs, adopting the obtained error margins to the rest of our sample. For properties based on more approximate or ad-hoc descriptions (e.g., clumping, X-rays, β), we do not provide uncertainties in this study. While it is formally possible to explore this parameter space, e.g., in the context of a Genetic Algorithm coupled with a fast atmosphere calculation code (see, e.g., Brands et al. 2022, for applications with FASTWIND), codes like CMFGEN, which compute the full spectrum in the comoving frame in much more detail are very time consuming and computationally expensive. Therefore, in this analysis framework, uncertainties for each parameter for each star are much more challenging to be quantified. Usually one need to rely on the analysis of a few objects/models (e.g., Bernini-Peron et al. 2023) and/or usually can only be better constrained in dedicated studies for specific properties (e.g., Rübke et al. 2023).

3.1 Atomic data

For the detailed NLTE modeling of the atmosphere, we followed the same approach as described in Bernini-Peron et al. (2023), who analyzed Galactic BSGs cooler than 20 kK. However, as we also model hot BSGs in this work, we additionally consider higher ionization stages. The ion sets included in hotter and cooler models are slightly different, as certain ions are not significantly populated at some temperature ranges. We list our considered elements and ions in Table B.1.

The atomic data we used are those present in the default, publicly available, installation of CMFGEN (v.05.05.2017). Most of the sources used by the code can be found in Hillier et al. (2003), though.

3.2 Determination of the photospheric parameters

In the following, we describe the process to obtain the fundamental stellar (i.e., non-wind) properties, which are summarized and discussed later in Sect. 4.

Luminosity and reddening. To obtain the stellar luminosity L and the reddening E(B − V), we fitted the SED from the UV to the infrared using the flux points derived from the different magnitudes (listed in Sect. 2) and flux-calibrated spectra from ULLYSES and the XShooter counterpart when available. We adopt the extinction law by Howarth (1983) with R_V = 2.7 and we consider a fixed galactic foreground extinction with E(B − V) = 0.03 and R_V = 3.1 using Cardelli et al. (1989) law. The generally low extinction toward the SMC limits the diagnostic values of the absorption bump at 2175 Å. For inferring the E(B − V), we thus focus on reproducing the overall spectral energy distribution given by the shape of the flux-calibrated UV and optical spectra as well as infrared photometry. We estimate a conservative uncertainty of 0.03 for our derived E(B − V) values. As the reddening also enters the determination of the luminosity L, we adopt an error of 15% for our derived values of L.

Rotation and macroturbulence. The rotation v sin i and macroturbulence v_mac were obtained by using the IACOB Broad tool (Simón-Díaz & Herrero 2014). We used He and metal lines, namely He I λ4713 and Si III λ4552 for stars earlier than B2 and Mg II λ4481 for later types (see Fig. C.1 and Fig. C.2 for an example). The tool also provides the v sin i computed from the Fourier transform of the line profile (hereafter v_fft), which we use as a sanity check and to compute uncertainties (see Appendix C). As the metal lines are less affected by Stark and other additional broadening mechanisms (see, e.g., Simón-Díaz & Herrero 2007, 2014), we give preference to their resulting v sin i, v_mac, and v_fft and report them in Table 2. To quantify the uncertainty, we also obtain the v sin i, v_mac, and v_fft via He I λ4713 and compute the standard deviation between the He and metal lines for each the quantities (Δv sin i, Δv_mac, and Δv_fft). For all stars, we consider a minimum uncertainty of 25% in v sin i and of 40% for v_mac – the average relative standard deviation of these quantities in the sample. In case of the values being smaller than Δv_fft, we select the latter as the final uncertainty. This conservative approach avoids mathematically small errors and does not artificially underestimate the errors of targets whose uncertainties are larger. In the spectral analysis, the derived rotation is implemented by convolving an elliptical kernel, while for macroturbulence a Gaussian kernel convolution is applied (cf. Gray 2005). We further evaluated the overall spectral fitting, and if necessary made adjustments to the initially obtained values of rotation and macroturbulence. This happened in particular for AzV 187, where v sin i is close to 11 km s⁻¹, which is near the limit imposed by the Nyquist frequency for the star’s spectral resolution – see Simón-Díaz & Herrero (2007). Thus, the determinations were deemed as not reliable – see discussion in Appendix C.

Effective temperature and surface gravity. The effective temperature T_eff is obtained via evaluating the ionization balance of Si and He by fitting multiple ionization stages, when available. For the “hot” BSGs (B0 to B1.5) Si IV/Si III and He II/He I are used, namely Si IV λ4089, 4116 and Si III λλ4556,69,76 for silicon, and (ii) He I λ4471, He I λ4387 and He II λ4552 for helium. He II λ4686 was used as an auxiliary diagnostic as it can be affected by wind properties (Martins et al. 2015). For the “cool” (B2 to B5) and “cold” BSGs (later than B5) we aim to fit Si III/Si II and the He I lines (as the He II lines are non-existent), namely Si II λλ4128,32 and Si III λλ4556,69,76 for silicon, and especially He I λ4471, He I λ4387 and He I λ4713. Additionally, other lines such as Mg II λ4481 and the Balmer lines are taken into account. In some cases, T_eff went through further revisions during the determination of He abundance. As Vink et al. (2023) discusses, there is evidence that the Si and Mg content in SMC could be individually lower than the 0.2 Z_⊙ scaling −0.14 and 0.17, respectively. However, this does not introduce major uncertainties in the temperature, as it affects Si lines more or less equally. Testing models with varying temperatures only, we find in general an accuracy of 1 kK. To account for influence of other parameters, we consider a standard error of 1.5 kK, which is within the range of similar studies. We specify otherwise in the case of targets which we encountered more difficulties finding a satisfactory spectral fitting. The log g values in the photosphere are obtained by fitting the wings of the Balmer lines, given their high sensitivity to pressure in the stellar photosphere (see, e.g., Gray 2005). Additionally, some lines such as Si III λ4550–65–79 and He I λ4471 are also affected and were used as a sanity check. Likewise for the temperature, we estimate an uncertainty of 0.1 dex for log g.

Spectroscopic masses and radii. These quantities were indirectly determined from log g, T_eff and L. The uncertainty for mass M and R are, respectively, 35% and 15%. These values were obtained considering our errors in the primary quantities to follow normal distributions (see Bernini-Peron et al. 2023). Because of the small relative errors in the adopted distance (e.g., ≲0.01 dex, Graczyk et al. 2020), the luminosity is very well constrained in the case of SMC stars. Hence, log g is the dominant source of the mass error.

Microturbulence. The microturbulence near the photosphere ξ_phot, reflecting small-scale perturbations (see, e.g., Moens et al. 2022; Debnath et al. 2024), is estimated mainly by analyzing the strength of Si III lines as these are sensitive to the velocity field (e.g., Catanzaro & Leone 2008). The intensity of other lines such as He, H (mostly in the core of the line), and other metals are also affected and thus used as sanity checks. For the Si lines, the most prominent in the optical and the UV are considered, namely, Si III λλ4553, 69, 76 and lines around 1300 Å. ξ_phot affects temperature and surface gravity diagnostics, T_eff and log g were re-adjusted if necessary. Especially for the earlier-type stars, the central depth of the Balmer lines is also affected by the density of the inner wind, thus introducing a sensitivity to the mass-loss rate, clumping, and the velocity gradient – see Sect. 3.3. The microturbulence is set to increase linearly with the velocity across the wind, following: $$\xi (r)={\xi }_{\text{phot }}+\frac{({\xi }_{max}-{\xi }_{\text{phot }})v(r)}{v_{\mathrm{\infty }}},$$(1)

where ξ_phot is the turbulence near the photosphere and ξ_max is the turbulence when the wind reaches v_∞. We obtain an accuracy of ~3 km s⁻¹. However, in cases where the actual Si abundance is lower than the adopted baseline, a higher microturbulence would be inferred to reproduce the same line strength. To include this effect in our uncertainty estimates, we assume a higher margin of about 5 km s⁻¹.

He and CNO surface abundance. As presumably evolved objects, BSGs are expected to present surface alterations in H, He, and CNO. In this study, to obtain the helium mass-fraction Y, we adjusted the He abundance aiming for a better fit of the many He I (and He II for earlier types) lines along the optical range. From that procedure, we estimate an error of a factor of 2. As discussed in Hillier et al. (2003) or Crowther et al. (2006), a precise determination of the He abundance in supergiants is particularly challenging and we thus give this rather large margin. Enhanced He abundances should thus be rather seen as a qualitative sign of enrichment rather than a precise number yielding hard constraints. To determine the CNO abundances, we employ the following sets of lines: For carbon, we fit C II λ4267 and C II λλ6578, 82 in late BSGs and C III λ4069 and C III λλ4648–50 in early BSGs. For nitrogen, the N II series of lines at 4600 Å (in late BSGs), N III 44097, N II λ4447, N III λλ4510–20 (early BSGs) were used. N III λ4097 is heavily affected by the temperature and, therefore was not considered as a main diagnostic. In the case of oxygen, the optical region between 4000 and 500 Å is crowded with oxygen lines. From these, we give most emphasis on O II λ4069–92, O II λ4590–96, 4661 and O III λ4367. The determinations of precise abundances for individual hot massive stars can be quite challenging, as different diagnostic lines may point to different abundances (Martins et al. 2024, XShootU Paper V) and fundamental parameters (T_eff, log g, and ξ) affect these lines significantly. Given that uncertainties in T_eff, log g, and ξ affect CNO lines, we consider an error of 0.3 dex also for the CNO abundances. In Fig. E.2 we show an example of how CNO diagnostic lines are affected by different abundances.

3.3 Determination of wind parameters

After determining the photospheric parameters, we focused on the analysis of the UV and recombination lines, which encode most of the wind properties. The resulting properties are listed further in Table 3.

Wind terminal velocity and microturbulence. The wind terminal velocity is obtained by fitting the width of the UV P Cygni wind lines, especially the absorption trough. For stars where the UV P Cygni profiles are very unusually shaped (e.g., AzV 104) or very weak (mostly for later type BSGs) the determination of the terminal velocity is not obvious. In those cases, we adopt the values that yield a better overall fit of the line. The adopted microturbulence at the terminal wind speed (ξ_max), which is computed in the formal integral to producing the observer’s frame spectrum, can influence the derived values of v_∞. To keep uniformity with previous works that used CMFGEN, we took a value of ~10% of the respective v_∞. For targets with unusual or very weak UV P Cygni profiles, we explored a wide range of values, reaching, in some cases, values of ξ_max ~ v_∞.

Mass-loss rates. To determine the mass-loss rates, we aim for a simultaneous fit of Hα and the wind P Cygni profiles in the UV. For some targets, fitting the profile with CMFGEN is not possible, in particular, due to a peculiar Hα shape. In these cases, we perform a different approach described in Appendix D. Additionally, Hα is severely affected by clumping parameters whereas the UV wind lines are in general very affected by the X-ray emission. From exploratory model calculations, we infer a standard uncertainty of 0.15 dex. However, when considering the effects of clumping – see below – we take a cautious approach and estimate a combined error of 0.40 dex.

Velocity profile. As CMFGEN does not compute the velocity stratification v(r) consistently with the radiative acceleration, the velocity structure needs to be provided. This is described by a modified β-law, describing the wind, connected to a quasihydrostatic structure that describes the inner atmosphere (subsonic regime) which is updated based on the estimated line acceleration. Wind lines, in particular optically thick lines are affected by different values of the parameter β (Lefever et al. 2023, e.g.,). In the BSG regime, Hα is notoriously sensitive to the density profile. However, other Balmer lines and P Cygni profiles are also affected, though with much less sensitivity than Hα. Thus we estimated β by aiming to find a better fit to Hα and used the UV features as auxiliary diagnostics. As clumping also affects significantly Hα, we also needed in many cases to revise our β values. In this study, we set the connection velocity v_con, reflecting the transition between the photospheric/subsonic and wind/supersonic regime, to be 10 km s⁻¹ by default, which is around 70% of the sound speed in the photosphere of BSGs. The value is lowered if necessary to ensure a monotonic velocity field (which is mandatory in co-moving frame atmosphere calculations).

Clumping parameters. In CMFGEN clumping is pictured as optically thin overdensities, whose diameters are smaller than the photons’ mean free path, and whose surroundings (interclump regions) are void. Mathematically, it is by standard implemented the following profile in radius stratification: $$f(r)=f_{\mathrm{\infty }}+(1-f_{\mathrm{\infty }})e^{-v/v_{\mathrm{c}\mathrm{l}}},$$

describing a scenario where clumping reaches its maximum in the outer wind. The two clumping parameters are f_∞ and v_cl, which, respectively, represent the volume filling factor (i.e., how clumped is the wind at r → ∞) and the characteristic onset velocity (describing where clumping starts to be relevant in the wind). These are obtained by aiming at a consistent fit of Hα (and Hβ in some cases) and the UV lines in general. However, for some targets, due to their unusual Hα profiles, we explored other approaches – cf. Appendix D. The Balmer lines as recombination lines are proportional to ρ²(r), while the UV P Cygni profiles as scattering lines are proportional to ρ(r) only. Therefore, as essentially perturbations to the density, clumping will impact much more the latter. However, as e.g. Bouret et al. (2012) has demonstrated, clumping also has some impact in the UV by altering the ionization stage of different ions via a radiative field which encounters a different distribution of matter. For the O and early-B targets, N IV and P V are particularly good diagnostics. As Bouret et al. (2005) discusses, these ions increase their population in the wind due to the presence of clumps (i.e., enhanced density). Therefore we aimed for a reasonable fitting of these lines when obtaining clumping. For later-type stars, clear diagnostic lines specific to this parameter are weak or absent, so we can only aim for a good overall UV spectral line fitting. Moreover, in later BSGs profiles that are sensitive to wind stratification are also affected by X-rays (Bernini-Peron et al. 2023).

X-ray parameters. The X-ray emission in the wind is by default parametrized in CMFGEN as a Bremsstrahlung component plus tables of emissivities generated by the APEC code (Smith et al. 2001). Three free parameters define the emissivity profile, namely T_X, which represents the typical temperature reached by the X-ray production mechanism; v_X, the characteristic velocity of the production of X-rays that reflects the point in the wind where the X-ray emissivity profile would reach ~30% of its possible maximum; and f_X, the X-ray volume filling factor. As this X-ray treatment is a simplification, T_X and v_X should not be interpreted as actual physical properties of the wind and their precise values do not bear larger physical significance. With the basic X-ray parameters, the code integrates the emissivity profile and gives the produced as well as the observable X-ray luminosity (hv > 0.1 keV). The latter is usually given in terms of the stellar luminosity as log(L_X/L). The X-ray luminosity (L_X) is constrained by fitting the UV lines. For stars earlier than B1.5, the main diagnostic used is the famous N V λ1238 line. This ion can only exist due to the presence of this additional ionization source. For cooler stars, C IV λ1550 becomes the main diagnostic, as analogously, C IV would not be present in winds of stars cooler than ~20 kK without the presence of X-rays. Conversely, the C II λ1335 profile without the presence of X-rays tends to be overpredicted, as such lower temperatures are more fertile to the existence of this less ionized stage. The “shock temperature parameter” T_X is set to be 1.0 MK for models with v_∞ ≥ 900 km s⁻¹ and 0.5 MK otherwise, motivated by the results of Bernini-Peron et al. (2023). Finally, the characteristic velocity parameter v_X is initially set to ~80% of v_∞ and changed if necessary to improve the fitting. For some targets, this value had to be severely changed to improve the fitting. In Sect. 5 we discuss these aspects in more detail.

4 Stellar properties

In the process of analyzing the spectra of BSGs, prior to the modeling phase, we noticed that some targets initially classified with similar spectral types had very similar optical spectra – which foreshadow similar photospheric properties, but looked very different on the UV. Conversely, some other targets present a similar UV spectrum between them, but very different optical spectra. This is illustrated in Fig. 3, where we can see that AzV 210 and AzV 472 have an almost identical optical spectrum (except for Hα) and similar luminosities, but very different UV lines. This highlights that their winds must be considerably different.

Contrary to that, AzV 22 and Sk 179 display very similar UV features but quite different optical lines. The differences between the latter two could be related to their difference in luminosity (or proximity to the Eddington Limit), since AzV 22 is much more luminous and has a lower log g, noticeable by its much narrower Balmer lines. By covering all of these cases in our analysis sample, we aim to get a representative overview of the different kinds of BSGs in the SMC.

Some stars within our sample display peculiar signatures in their UV, optical, or IR spectra – namely AzV 235, AzV 104, Sk 191, and AzV 362, which could not be properly modeled by CMFGEN. We examined the spectra from different epochs and found no clear signatures of multiplicity. Yet, we noticed variability in their Balmer profiles (except for AzV 104), which is especially strong for AzV 362. For Sk 191, a mild variability can also be seen in He I λ4471, which we interpret as a wind feature. Specific comments on individual targets are further given in Appendix D.

In Fig. 4, we provide an overview of the spectral fits in selected diagnostic ranges for the whole sample to illustrate the quality of the results and the trends across the BSG subtype range. The figure illustrates also some common issues in the quantitative spectroscopy of BSGs, such as the under-prediction of N III λ4097 in the earlier-type stars (e.g., Crowther et al. 2006; Searle et al. 2008). Full comparison plots between models and observations as well as the SED fits can be found in the Appendix F, available at the platform Zenodo⁴. The resulting photospheric properties we derived are listed in Table 2.

Fig. 3 Comparison between observed spectra of AzV 210 (B1.5 Ia) vs. AzV 472 (B1.5 Ia), and AzV 22 (B3 Ia) vs. Sk 179 (B3 II). For the comparison between the B1.5 stars, their optical and photospheric spectra are very similar to each other, but their wind lines at the UV and Hα are considerably different, despite the similar luminosity and temperature. In the juxtaposition between B3 and B8, the profiles in their UV are similar, while their optical spectra differ, especially due to broadening and different surface gravities. In each comparison, we list the old and the updated spectral types (Crowther, priv. comm.).

4.1 Spectral-Type calibration

As we analyzed a representative sample of SMC BSGs in this work, we can check the correspondence between the effective temperatures and spectral types. In Fig. 5, we show the comparison between T_eff and spectral type alongside the calibrations from Schootemeijer et al. (2021) for SMC supergiants, as well as with calibrations for Galactic BSGs from the studies by de Burgos et al. (2023) and Searle et al. (2008), who analyzed Galactic BSGs using FASTWIND and CMFGEN, respectively. We additionally plot results from Evans et al. (2004a), Trundle et al. (2004), and Trundle & Lennon (2005) for SMC BSGs as well.

In general, there is a good agreement between the spectral analysis results from our work and the literature when compared to the relation of Schootemeijer et al. (2021). Yet, up to a spectral type of B2, the T_eff of our sample BSGs fall below the calibration by Schootemeijer et al. (2021), which was based on previous literature results including those marked in fuchsia in Fig. 5. The divergence is largest between spectral types B0.5 to B1.5. Within this range, we also tend to find lower temperatures than the Galactic calibration by de Burgos et al. (2023). On the other hand, our results for this spectroscopic range agree better with the relation from Searle et al. (2008), even though falling slightly lower still. For spectral types later than B2, our temperatures appear to agree well with the Schootemeijer et al. (2021) relation, but diverge from Searle et al. (2008).

Fig. 4 Spectral fits of exemplary subsets of diagnostic lines for the sample BSGs. The black curves are the X-Shooter data for each star, whereas the colored curves are their respective CMFGEN models.

4.2 Evolutionary status of the BSGs

When comparing our derived luminosities and T_eff to previous spectroscopic analyses on SMC BSGs (e.g., Evans et al. 2004a; Trundle et al. 2004; Trundle & Lennon 2005), we notice no relevant systematic differences (see Fig. 6). Likewise, we also do not find systematic disagreements between our investigation and the study of Schootemeijer et al. (2021), who used Gaia distances to infer the stellar luminosities. Even the largest discrepancies in the parameters are within the expected frame of uncertainty when comparing different analysis methods (Sander et al. 2024, see).

With the luminosities and temperatures of our 18 BSGs determined, we can locate the stars in the Hertzsprung Russell diagram (HRD) and discuss their evolutionary context. In Fig 7, we plot the HR diagram of the BSGs together with the evolutionary tracks’ curves by Brott et al. (2011, hereafter B11) and Georgy et al. (2013, hereafter G13). Both sets of evolution models have similar rotation, but their treatment of convection is considerably different. The B11 models determine convection using the Ledoux criterion (Brott et al. 2011), while the G13 models use the Schwarzschild criterion (Ekström et al. 2012; Georgy et al. 2013). We find our sample stars to be displaced relative to the zero age main sequence, which is consistent with previous results in the literature and the expectation that BSGs are already evolved objects (e.g., Trundle et al. 2004; Evans et al. 2004c). However, this does not automatically imply that BSGs cannot be core-hydrogen-burning. As Higgins & Vink (2019) discuss in detail, evolutionary models with high convective mixing can produce tracks whose main sequence covers the BSG temperature and luminosity regime – see also Martins & Palacios (2013).

Notably, the B11 tracks can reach the BSGs earlier than B1 in the HRD. As discussed in Brott et al. (2011), the overshooting adopted in B11 tracks is calibrated for stars with initial masses M_init ~ 16 M_⊙ to match the TAMS inferred from the VFTS sample (Hunter et al. 2008). This aligns with more recent investigations that argue more massive stars may have progressively higher mixing (e.g., Castro et al. 2014; Martinet et al. 2021). As a consequence, the main sequence might reach down to T_eff ~ 20 kK – corresponding approximately to the B 1 spectral type (McEvoy et al. 2015).

Contrary to that, the G13 tracks interpret the HRD positions of the early BSGs as objects beyond core-H burning. The “terminal age main sequence (TAMS) hooks” of the G13 tracks happens at much hotter temperatures than those in the B11 tracks. This is rooted in the much lower convective mixing in the G13 high-mass tracks compared to B11. The G13 high-mass tracks are calibrated to match the TAMS width for intermediatemass stars (Georgy et al. 2013). According to the G13 models, more massive/luminous stars spend much less time as early-type BSGs (i.e., between ~30 and ~20 kK) than as later-type BSGs.

Below T_eff ≲ 19 kK, the redward evolution predicted by the G13 models slows down and the stars spend a decent amount of time still in the hot star regime, where central He-burning eventually sets in. This is different for the B11 models, which evolve considerably fast in the Hertzprung gap, as indicated by the scatter symbol distribution along both track sets which mark intervals of ~50 kyr. This discrepancy is due to the different convection criteria between the two evolution model grids. The Schwarzschild criterion in the G13 models produces a larger intermediate convective zone (ICZ) at the onset of hydrogen shell burning (see also Sibony et al. 2023; Josiek et al. 2024). The additional energy required to establish a larger ICZ is not available for the expansion of the star and thus the G13 models spend more time crossing the Hertzsprung gap than the B11 models.

In contrast to the Galactic BSGs discussed in Bernini-Peron et al. (2023), neither the B11 nor the G13 tracks predict a blueward evolution from the red supergiant (RSG) regime. Blue loops only happen for more massive G13 models and at temperatures more compatible with yellow supergiants/hypergiants and luminous blue variables (LBVs). For this regime, we do not have corresponding objects in our sample. Consequently, assuming the evolution framework of B11 and G13, none of the studied SMC BSGs would be post-RSG objects. However, certain mixing settings (e.g., Schootemeijer et al. 2019) or enhanced mass-loss in the RSG or LBV regime (e.g., Schootemeijer et al. 2024) could produce tracks evolving blueward from the RSG regime at SMC metallicity. Therefore, we cannot rule out that some of our targets are post-RSG objects. This would be supported by our finding of considerably lower masses compared to evolutionary predictions, which we further discuss in Sect. 4.4.

Fig. 5 Comparison between effective temperatures and spectral types. The red diamond symbols represent our spectroscopic determinations, whereas the hexagons are literature values from Evans et al. (2004a), Trundle et al. (2004), and Trundle & Lennon (2005). The plotted curves are the calibrations from Schootemeijer et al. (2021) for the SMC stars (solid blue curve), Searle et al. (2008, dotted green curve), and de Burgos et al. (2023, dashed dark green curve) for Galactic BSGs.

Fig. 6 Comparison between the derived effective temperatures and luminosities of our sample SMC BSGs with results from the literature, namely: Evans et al. (2004a, fuchsia hexagons), Trundle et al. (2004, pink hexagons), Trundle & Lennon (2005, orange hexagons), and Schootemeijer et al. (2021, black squares).

Fig. 7 HR diagram of the sample stars (HRD). Tracks of Georgy et al. (2013) (G13, thick orange lines) and Brott et al. (2011) (B11, thin gray lines) with initial rotation of veq/vcrit ~ 0.4 (i.e., ~300 km s−1). The small triangular and circular points on the tracks mark intervals of about 50 kyr for G 13 and B 11 respectively, and the red diagonal lines represent the iso-radii lines (1, 10, and 100, R⊙ from left to right).

Fig. 8 HR diagram showing helium enrichment (Y/Yinit, where Yinit is the baseline value of 0.27) in the B11 and G13 SMC tracks compared to our sample BSGs (diamonds). The tracks are the same as in Fig. 7.

4.3 Surface chemical abundances

In Table 2, we also list the surface helium mass fraction Y. Several of our objects show signs of He enrichment, in particular the BSGs above 19 kK, corresponding approximately to spectral types earlier than B1.5. This is qualitatively in line with the recent finding by de Burgos et al. (2024b), who also found hardly any He-enriched objects in the Galactic BSGs below this temperature. Although, it is important to note that for T_eff ≲ 20 kK, He abundance determination rely solely on He I. In our sample, the BSGs without He enrichment are in the clear minority in the hotter regime, which does not seem to be the case for the Galactic sample, albeit the work by de Burgos et al. (2024b) mixes actual supergiants with objects from other luminosity classes. When considering the He abundances of both evolution models, we do not see the observed amount of He enrichment in the early stars (see Fig. 8), which hints at remaining issues. One possibility is that the models do not assume sufficient mixing. Schootemeijer et al. (2019) showed that for certain conditions of overshooting and semiconvection, it is possible to produce BSGs evolving blueward from the RSG regime, which would be richer in surface He. Still, these models cannot reach the T_eff of our earlier targets. Recent work from Gormaz-Matamala et al. (2024) showed that adopting a lower $\dot{M}$ (obtained via evolutionary models whose mass loss is computed consistently) in the main sequence may produce stellar models that reproduce the He enrichment for the early BSGs. Another possibility, outside of the single-stellar evolution scenario, is that the hotter BSG population might have a more complex origin, like stripping or mergers (Klencki et al. 2022; Menon et al. 2024) – which we discuss further below.

When considering the evolution of the rotation and CNO abundances, we find a clear displacement from the SMC carbon and nitrogen baseline⁵ (dashed lines in Fig. 9). This confirms previous findings and underlines that our targets are evolved objects. For oxygen, we do not observe much variation and find values close to the baseline (see panel E in Fig. 9), implying that the material in the outer layers has not reached full CNO equilibrium. The evolution models by B11 and G13 each have their own, slightly different, initial baseline abundances. Yet, the effect on the resulting tracks is likely minor compared to the larger differences in the physical treatments inherent to the two different evolution model codes. The G13 tracks and the slow-rotating (i.e., with initial rotation velocity $v_{\text{rot }}^{\text{init }}\lesssim 100\mathrm{k}\mathrm{m}{\mathrm{s}}^{-1}$) B11 models reproduce the data better than the fast-rotating B11 models. Similar to our comparison in the HR diagram in Fig. 7, it is evident that the G13 evolutionary models predict a fast evolution through the early-BSG stage (panel A in Fig 9).

Notably, our derived v sin i are lower than previous determinations (Evans et al. 2004a; Trundle et al. 2004; Trundle & Lennon 2005), whose values are higher than the rotation predicted from G13 for the late-O and BSG temperature range. This can be understood when considering that earlier work conflated v sin i and v_mac in their total line broadening, which can overestimate the rotation and thus explain their sistematically higher values of v sin i.

Our obtained v_mac is generally high compared to the derived v sin i, which is consistent with the results obtained by Simón-Díaz et al. (2017) for Galactic OB stars. Moreover, as e.g., discussed in Sundqvist et al. (2013), disentangling rotation and macroturbulence is considerably more challenging for stars with low rotation (v sin i ≲ 50 km s⁻¹). Additionally, as discussed in Simón-Díaz & Herrero (2014), the photospheric microturbulence also hampers an accurate determination of v sin i below 40 km s⁻¹. Hence, the derived v_mac could even be underestimated if v sin i is overestimated for our slowest rotating targets. A lower rotation for the late BSGs would increase the agreement with the predictions from G13, whereas the disagreement would be aggravated for the earlier ones.

We can only obtain a reasonable compromise between the stellar rotation, T_eff, and partially the involved evolution timescales, when we consider B11 tracks with lower initial rotation, preferably with $v_{\text{rot }}^{\text{init }}=60\mathrm{k}\mathrm{m}{\mathrm{s}}^{-1}$. However, for the chemical surface abundances, the B11 tracks with an initial rotation of 120 km s⁻¹ would be preferred. This occurs especially for nitrogen, which also is very similar to the values predicted by G13. Only for oxygen, the B11 models of both 60 and 120 km s⁻¹ give a better representation of our results than the G13 models. In general, none of our sample stars exceeds a v sin i of 100 km s⁻¹. Our object with the highest (projected) rotational velocity (85 km s⁻¹) is Sk 179, a B6 supergiant that defies the general trend seen in Fig. 9. However, even this value would not classify the star as a fast rotator. We do not see any fast rotating star in our BSG sample, not even above 21 kK. While fast-rotating BSGs are also not common at Galactic metallicity, there is an observed group of them, in particular above the presumed bi-stability jump (e.g., de Burgos et al. 2024a).

In Bernini-Peron et al. (2023), we discussed stellar mergers as a possible solution for the presence of cooler BSGs (B2 to B5) in the Milky Way, given that evolutionary models predict a fast evolution in this regime. In the galactic context, mergers would produce core-He-burning objects, spending a considerable time of their post-main-sequence evolution in the observed parameter range. Given the high binary fraction among massive stars (e.g., Mason et al. 2009; Sana et al. 2012, 2013), such a scenario would not be unlikely. Specifically, for apparently single stars, the merger incidence considering coalescence events at different moments of the evolution can amount up to 20% (de Mink et al. 2014). To get a basic idea of the enrichment of a merger scenario, we compare the N/C vs. N/O ratios derived for our SMC BSGs with those from Large Magellanic Cloud (LMC) BSG merger evolution models from Menon et al. (2024) in Fig. 10. The figure also shows data points for known LMC and SMC stars with (partially) stripped envelopes are plotted for comparison. Despite agreement within the large error bars, there seems to be a systematic shift between our BSGs and the merger models. The merger models further overpredict the He enrichment, notably in contrast to the underprediction by the G13 and B11 (single-star evolution) models. In terms of the predicted mass range, the LMC merger models predict similar values as obtained in our spectroscopic analysis. Yet, it is hard to draw any firm conclusions as none of these diagnostics can only be explained by a merger scenario. Moreover, tests for the SMC are necessary as the metallicity affects how the stars ought to interact (Klencki et al. 2022), and thus at which stage they may merge or not.

Considering the stripped OB stars and our sample BSGs, we find a good agreement in their He enrichment, except for 2dFS-163 (O8 Ib). We further find that part of the stripped sample shows much higher N/C and N/O ratios. Yet, we see that two objects align well with the BSGs. One is an early O-giant primary and the other has a primary star with a BSG spectral type (B1.5 Ia, Ramachandran et al. 2024), but is located in the LMC. The masses of our sample BSGs are higher than the stripped stars found by Ramachandran et al. (2024), but overall they agree with those from the more massive stars analyzed by Pauli et al. (2022) and Rickard & Pauli (2023), even though the latter are O giants.

Fig. 9 Comparison of effective temperature, rotation, and CNO enrichment–depletion between G13 (orange tracks with triangle symbols, Georgy et al. 2013), B11 (gray, blue, and pink tracks, Brott et al. 2011), and our spectroscopic determinations (white diamonds). In this figure, we include tracks with initial masses from 12 to 50 M⊙ (to 40 M⊙ for the G13 set). The tracks’ equatorial rotation velocities are scaled by π/4 to account for the stars’ inclinations (e.g., Hunter et al. 2009). As in Fig. 7, the points along the tracks mark intervals of about 50 kyr. Baseline CNO values of B11 models are respectively 7.37, 6.50, and 7.98, while for G13 those are 7.56, 6.95, and 7.83. The dashed horizontal lines mark the CNO average baselines listed in Vink et al. (2023, Table 2) determined as the average values from previous studies. Panel A and B show the evolution of the rotation against Teff while panels C, D, and E show respectively the evolution of the surface C, N, and O vs. the rotation.

Fig. 10 N/C vs. N/O diagram. The dashed lines show the CN and CNO equilibrium curves for the SMC (in fuchsia) and LMC (in black) using the respective baseline abundances from XShootU Paper I (Vink et al. 2023) and from Hunter et al. (2009), which was the same used by the LMC merger models of Menon et al. (2024, circles). The star symbols with green contours represent the ratios determined for SMC stripped stars from Pauli et al. (2022, six-pointed star), Rickard & Pauli (2023, seven-pointed star), and XShootU Paper VIII (Ramachandran et al. 2024, five-pointed star). The diamonds show the ratios of our sample BSGs and the fuchsia hexagons represent ratios determined by Evans et al. (2004a); Trundle et al. (2004); Trundle & Lennon (2005). Those do not provide He abundance determination, whose enrichment relative to the baseline is represented by the color scale applied to the symbols.

Fig. 11 Comparison between evolutionary and spectroscopic masses. The color bar represents the effective temperature of the targets with lighter colors reflecting lower Teff. The fuchsia hexagons correspond to the data of Trundle et al. (2004) and Trundle & Lennon (2005) and the dotted blue lines connect the same stars. The darker horizontal orange and gray bands represent, respectively, the intervals of 1σ obtained for G13 and B11. Likewise, the lighter colored bands are for 2σ.

4.4 Spectroscopic and evolutionary masses

In Fig. 11, we show a comparison between the spectroscopic and inferred evolutionary masses for each BSG in our studied sample. The evolutionary masses were determined by minimizing the χ² between the spectroscopically derived stellar parameters and the available tracks of G13 and B11 each considering T_eff, log g, and L (Eq. (7) in Schneider et al. 2014). This method, despite not considering interpolated tracks between the available set, is more systematic than what is commonly done by comparing visually in the HRD and/or log g-T_eff diagram (see, e.g., Bouret et al. 2021, for O stars).

Globally, there is a clear trend of higher evolutionary masses relative to the spectroscopic mass, especially below M_spec ≈ 30 M_⊙ (G13) and M_spec ≈ 45 M_⊙ (B11), respectively. As shown in the figure, the B11 comparison trend aligns very well with the discrepancy obtained by Trundle et al. (2004) and Trundle & Lennon (2005), especially the comparison with the B11 tracks. Juxtaposed to the Trundle et al. (2004) results, the mass discrepancy is slightly lower for the comparison with G13 tracks. The more and less intense shaded regions represent deviations from unity of 1σ and 2σ respectively.

When considering the mean results, there is a decent agreement between the spectroscopic and evolutionary masses for both track sets with G13 yielding a ratio within 1σ, whereas the B11-derived mean mass ratio only agrees within 2σ. The amount of scatter is similar to what Bouret et al. (2021) found for SMC O stars, but unlike them, our results show the aforementioned systematic shift. Similar to previous works which analyzed BSGs both in the Milky Way and in the SMC, we further find that for higher masses (or luminosities) the ratio between M_evol and M_spec tends to swap and the evolutionary models predict slightly lower masses (see, e.g., Trundle & Lennon 2005; Searle et al. 2008). This did not occur for the SMC O stars studies by Bouret et al. (2021), who attributed the discrepancy to a lack of photospheric turbulence in their CMFGEN O-dwarf models. This could lead to an underestimation of log g – as also discussed in Paper IV. Given that BSG atmospheres are more extended than OB-dwarfs (which implies higher scale heights), they are expected to have higher turbulent velocities, up to ~20 km s⁻¹.

Fig. 12 Comparison between models of AzV 488 (B0.5 Ia) with and without X-rays. The green dashed curves show the model without X-rays and the green full curves are the model with X-rays included. For this model, the log(LX/L) −7.4 with a TX parameter of 1.0 MK. All the other properties are identical between each model.

5 Limits of the X-ray luminosities

The UV spectra of BSGs reveal the presence of highly ionized ions (e.g., N V λ1240 and O VI λ1038), which are incompatible with their effective temperatures (≲30 kK) and the inferred wind stratification (see, e.g., Puebla et al. 2016). This “superionization” can only be modelled when including a hot component in the form of an X-ray field into the atmosphere calculations, providing therefore an indirect evidence for the presence of X-rays in the wind of BSGs.

In Fig. 12, we show a comparison between the models for AzV 488 (B0.5 Ia) with and without X-rays where we illustrate their impact on diagnostic wind lines. The biggest impact of the X-rays is seen in N V (as discussed before) and C IV. In the case of C IV, the presence of X-rays is necessary to reach the observed profile saturation – even though this comes at the cost of an overprediction of the emission part. The Si lines are barely affected, as well as N IV λ1720, which shows that in this regime these lines are good mass-loss rate and clumping diagnostics. Regarding the O VI and P V lines, we noticed no significant changes after the inclusion of X-rays as these lines look essentially photospheric.

For Galactic BSGs later than B2, Bernini-Peron et al. (2023) showed that the inclusion of X-rays is paramount to reproduce most of the superionization ubiquitously observed in the UV of such stars (e.g., Cassinelli & Olson 1979). They also showed that the X-ray luminosity in general falls slightly below the “standard relation” of log L_X/L = −7 (Sana et al. 2006) observed in many early-B and O supergiants. Unlike Galactic BSGs, the SMC targets show much less developed UV P Cygni profiles, an expected consequence of their low metallicity. However, one can still notice a clear superionization due to the presence of (i) C IV λ1548–50 and N V λ1238–42 as P Cygni profiles (the latter especially at the earlier type spectra), and reciprocally, (ii) C II λ1335 and Al III λ1855 as a weaker profile. In comparison to the Galactic targets though, the L_X/L ratio of the SMC BSGs seems to be similar, suggesting that the X-ray luminosity in these objects could be independent of metallicity.

In Bernini-Peron et al. (2023), it was further shown that even the very low log(L_X/L)~−12 with temperatures of ~0.1 MK produced successful simultaneous fits of C II and C IV lines of cool Galactic BSGs. However, the presence, even though weak, of N V in the observations, which could not be produced by these models, showed that these values are likely not real. For the SMC stars, on the other hand, this ion is not observed in targets later than B2 (T_eff < 19 kK), which prevents us from discarding those extremely low L_X/L ratios. Still, our modeling efforts apply the more likely value of 0.5 MK for the cooler BSGs in this work, which provides a more realistic upper limit for L_X/L. This is also motivated by the fact that N V is seen in the denser wind of the only B-hypergiant in our sample (AzV 78, B1 Ia+).

In Fig. 13, we compare models with and without X-rays for AzV 187 (B3 Ia), a star with T_eff = 16.5 kK. It is evident that the inclusion of X-rays improves the fit of all the lines, with the possible exception of Si IV λ1400. The absence of this extra ionization source would yield a distinct population of C II in the outer part of the wind, which is not observed in the C II λ1335 line, which appears essentially photospheric. Reciprocally, the models predict a photospheric profile for C IV λ1550, in clear disagreement with the clear P Cygni profile in the observations.

Similar to C II, Al III is predicted too strong in the model without X-rays. By including it, we managed to improve the fit, although the wind feature in the doublet is still too strong. An even higher ionization could improve the fit of this doublet, but would not align with the over-prediction by our model of Si IV λ1400, which, despite matching the width of the absorption trough better than without X-rays, predicts too much emission and a saturation. This problem in simultaneously modeling Si IV and Al III is similar to what occurs for the Galactic BSGs in Bernini-Peron et al. (2023). As they discuss, optically thick clumps, which have strong evidence to be present in BSGs even at low metallicity (see, e.g., Prinja & Massa 2010; Parsons et al. 2024), may improve this slightly, but more in-depth studies to understand how different optically thick clumping formulations affect spectral lines in BSGs are necessary.

Even the coolest star in our sample, AzV 343 (B8 Iab, one of the latest BSGs in ULLYSES SMC sample), for which we determined an effective temperature of 12.5 kK, shows spectroscopic features which can only be reproduced when including X-rays. This suggests that even much cooler BSGs with very low terminal velocities harbor somehow an extra source of ionization in their atmospheres. However, explaining these in the paradigm of shock-generated X-ray emission is challenging since winds as slow as v_∞ ≲ 300 km s⁻¹ should not produce the necessary high-velocity dispersions to generate X-rays. Specifically, from the relation T_X ~ (Δv/300 km s⁻¹)² · 10⁶ K, presented by Cohen et al. (2014), the velocity dispersion of the shocks Δv would need to be of the same order of v_∞ (or larger for some of the later-type BSGs) to be able to produce X-ray temperatures of T_X ~ 0.5 kK.

Considering the full sample, we plot the relative X-ray luminosity as a function of the stellar temperature in Fig. 14, color-coding also the terminal velocities. Our results indicate that BSGs with v_∞ < 900 km s⁻¹ have log(L_X/L) ≲ −8. Except for AzV 104, this regime also coincides with T_eff ≲ 22 kK. This aligns with the findings reported by Berghoefer et al. (1997) (but see Nazé 2009), that in spectral types later than B1, BSGs fall in general below the X-ray detectability limit.

On the other hand, there is also the possibility that our result might be biased by our initial setup where models with v_∞ > 900 are assigned with T_X = 1.0 MK (after Bernini-Peron et al. 2023). In Fig. 15, we show a comparison between models for AzV 488 with the same X-ray luminosity but employing different T_X. It is noticeable that the model with T_X = 0.5 MK produces a much stronger N V λ1240, the main diagnostic for L_X/L for these stars. Only when we reduce the X-ray luminosity (dashed line) by a factor of five – thus becoming similar to what we find for cool BSGs – we can obtain a compatible profile with the higher T_X. This is also a manifestation of the parameter degeneracy of the current X-ray implementation. Under these conditions, our results point to limits of −8.0 ≳ L_X/L ≳ −7.0 for early BSGs (< B1.5, not considering AzV 104 lowest value) and −9.0 ≳ L_X/L ≳ −7.5 for later BSGs. This relation however might be rooted more in the wind terminal velocity (which tends to be proportional to T_eff, see Hawcroft et al. 2024; Parsons et al. 2024) than in the temperature itself. Even AzV 104, which is an early BSG but has a very low v_∞, has a low L_X/L regardless of whether we employ a T_X = 0.5 or 1.0 MK.

The derived X-ray luminosites of all stars in our sample are well below log(L_X[erg s⁻¹]) < 32.5 which correspond to X-ray fluxes < 5 × 10⁻¹⁶ erg s⁻¹ cm⁻¹. This is far below current point source detection limits in the SMC. For instance, “The X-ray point-source catalogue” obtained during the XMM-Newton survey of the SMC has a flux limit for point sources 10⁻¹⁴ erg s⁻¹ cm⁻¹ (0.2–4.5 keV band) (Sturm et al. 2013). Similarly, even deep Chandra observations of the SMC (Laycock et al. 2010; Oskinova et al. 2013) were not sensitive enough to detect BSGs in the SMC. Therefore, it is not surprising that none of our sample stars are detected in X-rays so far. In the region of 30 Dor in the LMC, which is closer to us, X-ray properties of OB stars appear generally consistent with Milky Way counter-parts (Crowther et al. 2022). Therefore, it is possible that the X-ray emission in OB winds does not change significantly with metallicity.

We conclude this section by highlighting that since direct X-ray measurements for SMC BSGs are currently unfeasible, the only way to obtain information on the X-ray content of such stars is via UV spectroscopy. Within the current analysis framework, big systematic studies focused on specific targets (analogous to e.g., Puebla et al. 2016) could improve the constraints on L_X/L by investigating the effects of each parameter independently. Additionally, more sophisticated formulations need to be tested (e.g., two-component winds, multicomponent plasmas, different emissivity laws, see Zsargó et al. 2008). However, as the theoretical understanding of how X-rays are generated in such stars is reaching its limit, detailed simulations of the wind launching in these stars are required to fundamentally constrain any such approximative treatments.

Fig. 13 Comparison between models of AzV 187 (B3 Ia) and AzV 343 (B8 Iab) with and without X-rays. The dashed curves show the model without X-rays and the full curves are the model with X-rays included.

Fig. 14 X-ray luminosity compared to effective temperature. The thick diamonds with numbers represent the SMC targets AzV (or Sk for 179 and 191), while the circles with letters, in alphabetical order, represent the Galactic targets ε Ori (B0 Ia, Puebla et al. 2016), κ Ori (B0.5 Ia, Huenemoerder et al. 2011; Haucke et al. 2018), 9 Cep (B2 Ib), 55 Cyg (B2.5 Ia), o2 CMa (B3 Ia), and 64 Oph (B5 Ib/II, Bernini-Peron et al. 2023, as for the latter three). The colors of the circle and broad diamond symbols represent the wind terminal velocity. The purple thin diamonds represent the alternative values considering a TX = 0.5 MK for the corresponding early target connected by gray dashed lines. Likewise, the fuchsia crosses represent the alternative value for AzV 104 if considering a TX = 1.0 MK.

Fig. 15 Comparison between models with different relative X-ray luminosity and shock temperature parameters. The thin gray line is the normalized UV spectrum and the colored lines are the CMFGEN models. The blue filled line is the model with TX = 0.5 MK and log(LX/L) = −7.4. The purple dashed line has the same temperature, but with a lower log(LX/L) = −8.0. The red dot-dashed line has the same log(LX/L) = −7.4, but a higher TX = 1.0 MK. The different LX are obtained by varying the X-ray emission filling factor. All the additional X-ray parameters are kept the same for all the models.

6 Wind properties

The wind properties of our models are compiled in Table 3.

6.1 Wind terminal and turbulence velocities

In Fig. 16, we compare our inferred terminal velocities with previous literature results. For the small sample overlapping with Evans et al. (2004a), who also derived their v_∞-values by fitting CMFGEN models to the UV spectral lines, we find an excellent agreement. Hawcroft et al. (2024) employed the SEI method (Hamann 1981; Lamers et al. 1987) to obtain homogeneous values for v_∞ and ξ_max of most of the ULLYSES OB sample. In general, there is no striking discrepancy between their values and our analyses for the terminal velocities, although a fraction of our values is systematically lower, even when accounting for the error bars. We additionally compare our v_∞ determinations with those from the recent study by Parsons et al. (2024), who analyzed BSG winds using the SEI method to constrain optically thick clumping. We find that except for AzV 104, all of the v_∞ values agree well with Parsons et al. (2024) within the error margins. Three key differences in their work compared to Hawcroft et al. (2024) are that (i) Parsons et al. (2024) use different lines (Si IV and Al III instead of C IV), (ii) do not enforce β = 1, and (iii) consider the individual radial velocities of the stars. In general, Hawcroft et al. (2024) assigned the lowest quality flag to many targets with v_∞ < 1000 km s⁻¹, so that it is not surprising that more detailed follow-up studies that consider lines beyond C IV 1550 Å do find differences in the derived v_∞ values. The comparison of ξ_max instead, shows a large scatter between the studies accompanied by a lower trend of our values as well.

When we consider v_∞ + ξ_max, which roughly translates to the width of the UV P Cygni absorption troughs, we notice that the comparison trend is tighter, also when looking at the temperature, represented by the color gradient of the scatter points in Fig. 16. However, our values are still systematically lower than their determinations.

In the lower right (D) panel of Fig. 16, we check the agreement of the effective escape velocity $$v_{\mathrm{e}\mathrm{s}\mathrm{c},\mathrm{\Gamma }}=\sqrt{\frac{2GM}{R}(1-{\mathrm{\Gamma }}_{\mathrm{e}})},$$(2)

where Γ_e is the Eddington parameter considering only the electron scattering opacity. Again, we obtain lower values than those reported by Hawcroft et al. (2024). Interestingly, one can see a curved trend, where there is more agreement among the extreme early- and late-BSGs and less in the middle. In their study, v_esc,Γ was obtained via Bayesian analysis with B11 evolutionary tracks (BONNSAI, Schneider et al. 2014), where the masses are the corresponding evolutionary masses. As we consider spectroscopic masses for computing v_esc,Γ, the systematic differences may be a consequence of the mass discrepancy discussed in Sect. 4.4.

Temperature and wind terminal velocity relation. With T_eff and v_∞ obtained from the spectral fits, we compare our results to the v_∞(T_eff)-relation derived by Hawcroft et al. (2024) in Fig. 17. Considering our total BSG sample (red diamonds), we plot the derived T_eff and v_∞ values with the corresponding trend being reflected as a thick red line. Interestingly, the fit hardly changes when using the values for SMC O and B supergiants from Hawcroft et al. (2024, blue points and blue-thick line). When comparing our BSG and O+BSG SMC relations to the “total SMC” relation from Hawcroft et al. (2024), which takes into account all ULLYSES SMC OB stars in their sample (thin dot-dashed-blue line), it is evident that the slope of the latter is notably steeper. However, for our BSG dataset, the difference in the slope between our BSGs’ relation and Hawcroft et al. (2024) total SMC OB stars almost vanishes when excluding the peculiar slow-wind B0.5 supergiant AzV 104 (cf. Appendix D), as shown by the thick-dotted red line. What remains is a systematic shift toward higher v_∞ for the same T_eff. Additionally, the SMC O+BSG sample of Hawcroft et al. (2024) also contains hot, “slow-wind” stars effectively reducing their obtained slope. If these are excluded, the obtained slope of the O+BSG relation (filled-blue line) would be much closer to the slope of their total SMC OB relation (dash-dotted-blue line).

In the regime of cooler BSGs below ~19 kK, it is hard to identify a clear trend between v_∞ and T_eff. Still, it appears that a majority of the stars follow the general slope, but the amount of outliers, in particular below the curve, gets larger. In principle, the weaker UV profiles make it more difficult to derive accurate values for v_∞. Yet, the narrow profiles unequivocally show similar low velocities despite differences in T_eff.

One could imagine that such stars, which have a comparably slow wind for their temperature, might have an increased massloss rate, but as we will see in Sect. 6.3, this is not the case and the mass-loss rates are usually even lower. Hence, the low terminal velocities of the “outliers” in the v_∞(T_eff)-diagram cannot be attributed to objects lying closer to the Eddington limit such as hypergiants or LBVs (e.g., Sander et al. 2014, 2018; Vink 2018), but their winds seem to be suppressed instead for yet unknown reasons.

Bi-stability jump in velocity. The bi-stability jump has also been associated with a change in the ratio between the terminal velocity v_∞ and the effective escape velocity v_esc,Γ along the temperature domain (Lamers et al. 1995). To check whether our SMC BSG sample shows this, we plot the terminal velocity normalized by the escape velocity ($V_{\mathrm{e}}^{\mathrm{\infty }}:=v_{\mathrm{\infty }}/v_{\mathrm{e}\mathrm{s}\mathrm{c},\mathrm{\Gamma }}$) as a function of T_eff in Fig. 18. We find that the SMC BSGs seem to behave very similar to the Galactic BSGs. The existence of a drop in $V_{\mathrm{e}}^{\mathrm{\infty }}$ at T_eff ≲ 25 kK is a well-established phenomenon, but whether this is a sharp “jump” or a more smooth drop is not clear so far (cf. Markova & Puls 2008; de Burgos et al. 2024a). In Fig. 18, we show the relations from Kudritzki & Puls (2000) compiling previous literature (black dashed lines) for Galactic stars in comparison to the results inferred from our sample. As evident, the average values (thick red and blue lines) at each side of the jump are very close to the Galactic ones. However, the individual data points suggest that the “jump” in $V_{\mathrm{e}}^{\mathrm{\infty }}$ happens at a slightly lower T_eff ~ 19 kK compared to ~21 kK in the Milky Way.

If we scale the Galactic relation from Kudritzki & Puls (2000) for $V_{\mathrm{e}}^{\mathrm{\infty }}$ by a presumed metallicity-dependence of v_∞(Z), we can test whether the resulting scaled relation align with our dataset for $V_{\mathrm{e}}^{\mathrm{\infty }}$. Thus, we plot relations scaled by Z/Z_⊙^0.1 (Leitherer et al. 1992), Z/Z_⊙^0.2 (Hawcroft et al. 2024), as well as Z/Z_⊙^0.19 (hot side) and Z/Z_⊙^0.003 (cool side) from Vink & Sander (2021). The resulting curves are all in agreement with our data within the 1σ (shaded regions). On the hot side, however, the Hawcroft et al. (2024) and Vink & Sander (2021) scaling already start to deviate from most of the data. Notably, Hawcroft et al. (2024) discuss that their $V_{\mathrm{e}}^{\mathrm{\infty }}$ are mostly in agreement with the scaled Galactic values for T_eff > 21 kK, but this is likely due to their larger values for v_esc,Γ.

Overall, our empirical findings argue for a relatively weak scaling of v_∞(Z), even on the hot side of the bi-stability jump. This is further supported by our findings for v_∞(T_eff) discussed above. Considering our v_∞(T_eff)-relation without the “outlier” AV 104 (cf. Fig. 17), we more or less obtain the slope of Hawcroft et al. (2024), but shifted upward by ~400 km s⁻¹. Notably, this shifted relation essentially aligns with the LMC relation found by Hawcroft et al. (2024), implying that there is not really a metallicity difference in the terminal velocities of “normal” OB stars. Instead, the lower metallicity environment of the SMC seems to bear more objects that defy the v_∞(T_eff) and $V_{\mathrm{e}}^{\mathrm{\infty }}\left(T_{\text{eff }}\right)$ relations (“outliers”).

Fig. 16 Panel A: Comparison between derived terminal velocities (v∞) by this study and those from Hawcroft et al. (2024, small blue diamonds) and Parsons et al. (2024, green squares). The fuchsia hexagons are data from Evans et al. (2004a), which also obtained wind properties via optical and UV analysis with CMFGEN. Panel B: Comparison of the wind microturbulent velocities (ξmax) between this work and Hawcroft et al. (2024). Panel C: Comparison between the sum of v∞ and ξmax, which would roughly translate the width of the absorption component of the UV P Cygni profiles. Panel D: Comparison of the Γ-corrected escape velocities (vesc,Γ) between this study and Hawcroft et al. (2024).

Fig. 17 Terminal velocities as a function of effective temperature: The results from this work (red diamonds) are compared to the values from Hawcroft et al. (2024, blue circles). The thick solid red and solid blue lines show the resulting relation when considering only supergiants (BSGs and OBSGs). The thick red dotted line shows the relation for only BSGs, but excludes the peculiar AzV 104. The blue dot-dashed line shows the linear relation derived for all SMC OB stars (of all types and classes) in the Hawcroft et al. (2024) paper also considering the literature results. Conversely, the dark blue dashed line refers to the LMC.

Fig. 18 Escape-normalized terminal velocity vs. effective temperature. The dashed black lines are the relation presented by Kudritzki & Puls (2000) compiling previous data. The dot-dashed line is the relation presented by Krtička et al. (2021). The red and blue thick solid lines are the average value of our determinations and their respective shaded regions are the standard deviation. The thin dashed fuchsia lines are the Kudritzki & Puls (2000) relation scaled by (Z/Z⊙)0.1 metallicity dependence from Leitherer et al. (1992) and the thin dashed green line below it is the metallicity scaling found by Hawcroft et al. (2024). The brown horizontal lines indicate the velocity scaling proposed by Vink & Sander (2021).

6.2 Mass-loss rates and clumping

The derived mass-loss rates of our SMC B supergiant sample are plotted against their luminosities in the upper panel of Fig. 19, where we compare them to literature values for other OB stars in the SMC from Bouret et al. (2013, 2021), Ramachandran et al. (2019), and Rickard et al. (2022). We further added the recent SMC binary analysis results by Pauli et al. (2022, 2023) and Rickard & Pauli (2023). It is evident that the SMC is an environment where objects with a wide range of mass-loss rates can be found. The results from Rickard et al. (2022) of the O star population in NGC 346 unveiled a sample of non-supergiant O stars with very weak mass-loss rates that fell below previous theoretical expectations, raising the question whether a star actually loses a relevant amount of mass via stellar winds during the main sequence in this low-metallicity environment. However, our BSG results as well as studies of O supergiants and other evolved stars underline that wind mass loss is definitely not negligible after the main sequence.

One has to be careful when assigning spectral types to massloss regimes as the implicit assumption about the evolutionary status may not be correct, as for example recently demonstrated by Pauli et al. (2022) for the ULLYSES target AV 476. Still, already the color-coding of the different luminosity classes in Fig. 19 highlights that there is no simple scaling of $\dot{M}$ with the stellar luminosity L for all OB-type stars. Instead our B supergiants as well as the O supergiants analyzed by Bouret et al. (2013, 2021) and Ramachandran et al. (2019) are located above the bulk of the dwarfs in the $\dot{M}$-L-plane. Further detailed studies of more objects will have to see if the current overlap will remain or is largely an effect of missing wind information (due to the lack of UV data) or unresolved multiplicity.

In addition to the direct scaling with the luminosity L, it can also be insightful to consider the “modified wind momentum” (e.g., Kudritzki et al. 1999) $$D_{\text{mom }}=\dot{M}{v}_{\mathrm{\infty }}\sqrt{R/R_{\odot }},$$(3)

which is expected to scale with stellar luminosity L for radiationdriven winds following the (modified) CAK theory (Castor et al. 1975; Pauldrach et al. 1986; Kudritzki et al. 1989). The lower panel of Fig. 19 showing the corresponding quantities reveals that even within the regime of our analyzed BSGs there is a considerable scatter, despite the overall trend toward larger D_mom with higher L. For the highest luminosities, the scatter seems to get smaller, but there is also a smaller sample and a mixture of different types of objects. Notably, also the recently analyzed binary stars seem to follow the general trends. As the studies include both pre- and post-interaction binaries, this indicates that disentangled binary components do not severely differ in their wind properties from isolated stars with similar parameters.

Considering the clumping parameter f_∞, we do not see a clear trend along the BSG domain (see Table 3). For the cool Galactic BSGs, Bernini-Peron et al. (2023) found a tendency toward smoother winds (f_∞ → 1), in line with recent simulations (Driessen et al. 2019). We still find some of the cool BSGs in the SMC preferring a smooth-wind solution, but also have objects that require higher amounts of clumping, in particular AzV 362 and AzV 22, where f_∞ = 0.1 is demanded by their Hα profile – the latter, in line with the significantly high clumping Parsons et al. (2024) also found. However, Parsons et al. (2024) also found evidence of clumping in AzV 234, for which we required a basically smooth wind.

The hot BSGs show generally a higher degree of clumping, but again with notable exceptions. For several targets, there is no strong preference for a particular value of f_∞ (see also the discussion in Bernini-Peron et al. 2023) and we took this into account by assigning a larger error bar to the corresponding mass-loss rates.

As the BSG regime is a stage which most of massive stars ought to pass – and some stars may spend a considerable time in (cf. Sect. 4.2) –, it is important to benchmark theoretical predictions (which are applied to stellar evolution codes) with spectroscopic empirical values. Given that the temperatures of our BSGs span across the bi-stability jump region, we can use their properties to study the behavior of $\dot{M}$ versus T_eff. In the first panel of Fig. 20 we compare both quantities for each star as well as with their respective Γ_e.

In the figure, we notice an overall downward trend, which is related to the progressively lower luminosities at lower T_eff. There is also a correlation between Γ_e and $\dot{M}$, where stars farther from the Eddington limit are showing lower mass-loss rates, notably with a large scatter at ~19 kK. This temperature coincides with the drop in $V_{\mathrm{e}}^{\mathrm{\infty }}$ discussed above.

To better address whether or not there are “jumps” in $\dot{M}$(T_eff), it is worth correcting for the different intrinsic stellar luminosities in this comparison. In the lower panel of Fig. 20 we show the transformed mass-loss rates $\dot{M}$_t, introduced by Gräfener & Vink (2013), where the stars are scaled to the same luminosity and terminal velocity v_∞ – essentially tracing the emission measure of the mass loss. Numerically, this scaling of the mass-loss rates is defined as ${\dot{M}}_{\mathrm{t}}=\dot{M}\cdot f_{\mathrm{\infty }}^{-1/2}\cdot \left({10}^3\mathrm{k}\mathrm{m}{\mathrm{s}}^{-1}/v_{\mathrm{\infty }}\right)\cdot {\left({10}^6{L}_{\odot }/L\right)}^{3/4}$.

From Fig. 20, it is evident that there is only a weak decline of $\dot{M}$_t(T_eff). Moreover, the correlation between Γ_e and $\dot{M}$ is much more diluted. This result hints at a more constant behavior of the mass-loss rate with temperature.

Fig. 19 Mass-loss vs. luminosity (upper panel) and modified wind-momentum vs. luminosity (lower panel) diagrams of SMC OB stars of different sources and luminosity classes. The different symbols mark different studies and different colors/sizes mark different luminosity classes. Error bars are not shown for visual clarity. The literature values are slightly translucent whereas the diamonds represent our sample BSGs. The error bars in the lower right corner indicate the typical associated errors for our sample. The scatter points with downward black triangles denote upper limits.

6.3 Comparison to mass-loss rate recipes

Our derived mass-loss rates can be compared with different theoretical predictions (“mass-loss recipes”) from the literature that apply to the SMC metallicity, namely: Vink et al. (2001, V2001), Krtička et al. (2024, K2024), and Björklund et al. (2023, B2023). There are also recipes by Bestenlehner (2020) and Gormaz-Matamala et al. (2023), but they only apply to the O-star regime. We discuss the three applicable to BSGs in the following paragraphs:

V2001. Vink et al. (1999, 2000, 2001) derived mass-loss rates based on Monte-Carlo simulations for computing the radiative acceleration, ensuring a global hydrodynamical consistency of the wind (i.e., the consistency between radiative acceleration and mechanical acceleration can be violated locally, but the values integrated over the total wind are equal). This series of studies predicted that, due to the recombination of Fe IV into Fe III, the mass-loss rates ought to increase relatively sharply within a short effective temperature range between 27 kK ≳ T_eff ≳ 20 kK. This was associated with the bi-stability jump region and established the paradigm that mass-loss rates are expected to increase when massive stars evolve toward cooler temperatures. Later on, Monte-Carlo models with local hydrodynamical consistency endorsed the presence of the jump (Vink 2018; Vink & Sander 2021). Considering metallicity effects, the derived mass-loss rates are expected to scale with $\dot{M}$ ∝ Z^0.64 for B stars and $\dot{M}$ ∝ Z^0.69 for O stars.

K2024. Krtička et al. (2021, 2024) derived mass-loss rates using the METUJE stellar atmosphere code Krtička & Kubát (2017) to obtain BSG properties. Their models solve the velocity structure consistently with the computed line acceleration, which is performed in the Sobolev approximation, but scaled to account for the difference to a comoving frame calculation. They found that the wind of BSGs on the hot side of the bistability jump should be driven mostly by C, Si, and S while on the cool side Fe III should be dominant. For Galactic stars, they obtained an increase of the mass-loss rate on the cool side of the bi-stability jump. However, their increase is more gradual and less intense than the one in Vink et al. (2000), and happens at lower temperatures (T_eff ~ 19 kK), which seems to align better with recently determined values for Galactic BSGs by Bernini-Peron et al. (2023). The Krtička et al. (2024) findings also align with the results by Petrov et al. (2016), who found a jump at ~20 kK. When considering different metallicities, Krtička et al. (2024) find a scaling of $\dot{M}$ ∝ Z^0.6, which is in agreement with the empirical findings from Marcolino et al. (2022), who compiled literature CMFGEN results for O- and early-BSGs. Additionally, for lower metallicities, the increase in mass-loss rates is less pronounced, so such an increase manifests more as a mild interruption toward lower temperatures in an overall decreasing $\dot{M}$(T_eff)-trend.

B2023. Björklund et al. (2023) uses the comoving-frame version of FASTWIND to calculate a grid of hydrodynamically consistent models. The resulting mass-loss rates are then fitted assuming a multivariate power law to describe the dependence of $\dot{M}$ on the different stellar parameters. The physical and numerical basics of this method is similar to the approach in METUJE (Krtička & Kubát 2010, 2017) and PoWR (Sander et al. 2017, 2020), albeit the three codes differ severely in their detailed assumptions and implementations. However, the recipe derived by Björklund et al. (2021, 2023) differs considerably from the other mass-loss descriptions as Björklund et al. (2021) obtain essentially a monotonic decrease in the relation between T_eff and $\dot{M}$ without any bi-stability jump.

Discussion. A comparison between our derived mass-loss rates and the predictions from the aforementioned recipes is shown in Fig. 21, assuming 0.2 Z_⊙ for the metallicity of the SMC. Previous studies of OB dwarfs in the SMC (e.g., Ramachandran et al. 2019; Rickard et al. 2022) find mass-loss rates far below V2001. Yet, the picture is more complex for supergiants. Among the earliest targets, mostly with spectral types of B0 and B0.5, we find mass-loss rates more compatible with V2001 and K2024 predictions, which in general predict similar $\dot{M}$.

These results align well with those from Evans et al. (2004a) and Trundle et al. (2004), who obtained $\dot{M}$ values slightly higher than Vink et al. predictions for 0.2 Z_⊙. Using the data from Bouret et al. (2021), whose analysis method is similar to ours, we can extend our discussion to also include the regime of O supergiants (OSGs). The trend of the early-type supergiants aligning with the V2001 predictions remains.

Other studies using different codes and methodologies also find a similar behavior for OSGs. In particular, Ramachandran et al. (2019) found that their only OSG, unlike their O giants and dwarfs, agrees with V2001. However, only two of them should be considered seriously as they attribute a low confidence to the other three objects as they are suspected or known binaries. The early O3-supergiant component of SSN 7 in NGC 346 studied by Rickard & Pauli (2023) even has a mass-loss rate above V2001. Unless all these early supergiants are subject to further multiplicity, we suggest there is a dichotomy in the mass loss between dwarfs and supergiants as partially indicated by Fig. 19. Interestingly, the agreements with V2001 for the SMC are contrary to the findings for Galactic OB stars (of all types/classes), where Vink et al. (2000) overestimates $\dot{M}$ by a factor of ~3 when considering the effects of clumping (e.g., Najarro et al. 2011; Krtička & Kubát 2017; Hawcroft et al. 2021). This could imply that OSGs and early-BSGs might have a different $\dot{M}$(Z)-scaling compared to early OB-giants and -dwarfs.

While the root of this remains so far unclear, recent theoretical studies of the winds of WR stars (Sander & Vink 2020; Sander et al. 2023) revealed a more shallow metallicity scaling in the optically thick wind regime and a transition toward a steeper dependence when the winds get optically thin. Possibly, there could also be a different metallicity scaling between OB dwarfs, giants, and supergiants, assuming this coincides with different wind density regimes. Then, the scaling of the supergiant regime would likely be flatter than the scaling for the dwarf regime – where we might see the “weak wind phenomenon” (Martins et al. 2005; Marcolino et al. 2009) already at much higher luminosities in the SMC as recently discussed in Ramachandran et al. (2019). Remarkably, the proximity between empirical mass-loss determinations and the V2001 recipe for early BSGs appears to happen for LMC stars as well (Verhamme et al. 2024), which strengthens this scenario. Yet, it is worth considering the potential origin of the supergiants in the SMC. In the lower metallicity of the SMC, evolution models such as Klencki et al. (2022) predict a considerable number of (partially) stripped stars appearing as hot supergiants. Stripped stars have a higher L/M-ratio and thus show a stronger wind, which would naturally explain why the supergiants are not in line with the trend from dwarfs and giants. However, we do not find clear evidence that our sample BSGs are necessarily stripped.

Moving to cooler temperatures, most of the B1, B1.5, and B2 targets yield values between the predictions by V2001 (cool side) and K2024, with a slight preference for the former. Additionally, in this temperature regime, K2024 is very similar to the $\dot{M}$ values of V2001 (hot side) extended to the cool side (i.e., if we ignore the jump). For stars later than B2, the inferred mass-loss rates agree very well with the Krtička et al. (2024) SMC predictions, thus behaving very similarly to the cool BSGs in the Milky Way (Bernini-Peron et al. 2023).

Considering the overall trend, we notice a general agreement between our findings and the theoretical $\dot{M}$ predictions by K2024 and V2001 for OSGs and early-BSGs. The predictions of K2024 further agree with our derived values for most early and later BSGs. Interestingly, K2024 seems to underpredict $\dot{M}$ in the temperature interval between ~23 and ~19 kK (or between ~B1 and ~B2). Conversely, the predicted increase in $\dot{M}$ by V2001 yields values higher than our empirical data, albeit within the considered error margin for some targets. For the vast majority of stars, B2023 underestimates the empirical $\dot{M}$.

From these results, we do not find clear evidence for a sharp increase in the mass-loss rate when passing the bi-stability jump region predicted by V2001 toward cooler temperatures. This is in line with a recent larger-sample optical analysis of BSGs in the Milky Way by de Burgos et al. (2024a), who also could not find an increase of the mass-loss rate. For the SMC BSGs, we find the mass-loss rates on the cool side of the jump to be lower by more than one order of magnitude compared to V2001, despite the good agreement on the hot side. However, we notice that recipes ignoring the presence of a jump (B2023 and the extrapolated hot-side V2001) underpredict the amount of mass loss by about an order of magnitude. Hence, our results indicate that there is a change in the $\dot{M}$(T_eff)-behavior in the bi-stability region and not just a smooth downward trend.

Fig. 20 Upper panel: Mass-loss rates vs. effective temperature. The scatter points are color-coded by their Γe values. The three points outlined in fuchsia in each plot are the BSGs whose Balmer line profiles were peculiar, and thus have less reliable empirical values for the massloss rates. Lower panel: Transformed mass-loss rates (Gräfener & Vink 2013) vs. effective temperature. The color-coding is the same as in the upper panel.

Fig. 21 Upper panel: different mass-loss prescriptions applied to a generic star with the same properties but varying Teff for illustrative purposes, to facilitate the visualization of the predicted $\dot{M}$(Teff) behaviors. Vink et al. (2001) rates are represented by white circles, Krtička et al. (2024) by white squares, and Björklund et al. (2023) by white triangles. Middle panel: comparison between different mass-loss rate prescriptions and inferred values for the sample stars (black diamonds) and OSGs from Bouret et al. (2021, gray pentagons). The fuchsia-outlined diamonds are the targets where we did not find a satisfactory fitting of Balmer lines. The applied theoretical values have the same symbols as the upper panel. The stars are ordered according to their Teff in the upper axis. The lower axes show the spectral types of each star. Lower panel: similar to the previous panel, but showing the ratio of the predicted according to each recipe ($\dot{M}$th) to our empirical determinations. The symbols follow the same encoding from the previous panels.

7 Conclusions

In this work, we presented a thorough analysis of a representative sample of BSGs in the SMC. Using CMFGEN model atmospheres as well as the ULLYSES/XShootU large and homogeneous dataset plus additional archival data, we could analyze the UV and optical regime of 18 BSGs with unprecedented quality, spanning a temperature range between 12.5 and 27.5 kK and spectral types from B0 to B8. This comprises the widest multiwavelength study on B supergiants in the SMC where photospheric and wind properties were obtained. From the derived photospheric properties we could constrain their evolutionary context, from which we draw the following conclusions:

• The BSGs in general are evolved objects. From comparison with two evolutionary sets of single-star evolution models (G13 and B11), they are either main-sequence objects close to the TAMS (for the early-type stars) or already in the H-shell burning stage for the late-type objects;

• While the G13 tracks reproduce the properties of later-type BSGs considering physical properties and appropriate timescales, they cannot explain the hotter BSGs. On the other hand, the B11 tracks, which have higher overshooting, can reproduce the early BSGs in the HRD, although not the observed rotation and chemical abundances. The post-main-sequence evolution happens quickly in the B11 models, so they cannot explain the existence of late-type BSGs;

• The G13 models also yield a better overall agreement between the spectroscopic and the evolutionary masses compared to B11 and previous M_spec-M_evol comparisons (Trundle et al. 2004; Trundle & Lennon 2005). Additionally, we confirm the previously reported trend of higher evolutionary masses for less massive and less luminous BSGs;

• While none of our objects have clear spectral indication for multiplicity, we cannot rule out multiplicity-based origins for our BSG sample stars. The obtained CNO abundance ratios are roughly in line with existing scenarios. However, the He enrichment is lower than currently predicted in merger scenarios, albeit with specific predictions for merger scenarios in the SMC still missing;

• When comparing the BSG sample with the abundances of (partially) stripped stars, we find partial agreement. However, several stripped stars show stronger enrichment in nitrogen than our sample stars.

Our study marks the first work analyzing a big sample of UV spectra for BSGs in the SMC with a stellar atmosphere code. Moreover, this is the first time both clumping and X-rays have been included when studying UV wind diagnostic lines in this regime. Similarly to the Galactic study by Bernini-Peron et al. (2023), we can put the following constraints on the X-ray emission of BSGs within the SMC:

• We obtained a wide range of L_X/L-ratios for the different stars in our sample. In general, the early-type BSGs have −8.0 ≲ log (L_X/L) ≲ −7.0, while late-type BSGs have −9.0 ≲ log(L_X/L) ≲ −7.5;

• Our derived values are based on the parameterized ad hoc implementation within the wind-shock paradigm. We also find that there are degeneracies in the X-ray parameters with different values leading to a similar effect on the wind lines. However, some degeneracies can be lifted by making reasonable physical assumptions. For instance, it is not expected that cool BSGs, which in many cases have winds slower than 300 km s⁻¹, can produce shocks that are energetic enough to generate T_X ~ 1.0 MK. This in particular sets constraints for the later-type BSGs.

With the wind diagnostics given from the UV spectra and the sample covering a wide temperate range across the presumed bi-stability jump region, we gain the following insights into the wind velocities in this regime:

• We find good agreement between our obtained terminal velocities and those from the previous literature. However, we obtain systematically lower values relative to Hawcroft et al. (2024, XShootU Paper III) while the comparison with Parsons et al. (2024), both employing the SEI method, yields a good agreement;

• We find a clear jump in the escape-normalized terminal velocity at ~19 kK, very similar to what was found for Galactic BSGs at about 21 kK (Lamers et al. 1995);

• For cooler BSGs, we obtained a large spread in the escape-normalized terminal velocities, meaning that none of the proposed Z/Z_⊙-scalings is clearly preferred. For early-type BSGs, a small, possibly even vanishing scaling term seems to reproduce the empirical findings better, while the scaling from Hawcroft et al. (2024) and Vink & Sander (2021) predict velocities that are too low, albeit still within the 1σ error.

Lastly, we can draw the following conclusions on the wind mass loss of BSGs by comparing our derived properties with literature results and theoretical predictions:

• For the early-type BSGs, we find a good overall agreement between the different $\dot{M}$ recipes. However, we tend to find the best overlap with the predictions from Vink et al. (2001), similar to what Evans et al. (2004a) obtained, with some objects even exceeding the predicted values. This result aligns very well with the literature studies for SMC OSGs using different methodologies;

• Given that literature studies for SMC OB giants and dwarfs yield mass-loss rates clearly below the predictions from Vink et al. (2001), we suggest that there may be a dichotomy in the scaling of mass loss with metallicity between supergiant and non-supergiant OB stars;

• In the bi-stability jump region (i.e., between ~23 and ~19 kK), the derived $\dot{M}$ is found to be between Vink et al. (2001) and the other recipes that do not predict a jump. For stars cooler than this T_eff range, the Krtička et al. (2024) predictions that include a small increase align best with the empirical values, similarly to what was recently found for Galactic BSGs (Bernini-Peron et al. 2023; de Burgos et al. 2024a);

• Our findings do not support the smooth, monotonic decrease in $\dot{M}$(T_eff) suggested by Björklund et al. (2023). Instead, our results outline a scenario where $\dot{M}$(T_eff) does not have a clear bump or jump, but rather stays constant or increases very mildly between ~25 and ~16 kK;

• The escape-normalized terminal wind velocities for the SMC BSGs show a rather sharp drop at around ~19 kK. Notably, the temperature of this drop does not directly align with any change in the mass-loss rate behavior.

Acknowledgements

MBP, AACS, and VR are supported by the German Deutsche Forschungsgemeinschaft, DFG in the form of an Emmy Noether Research Group – Project-ID 445674056 (SA4064/1-1, PI Sander). MBP, AACS, and VR 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. JJ acknowledges funding from the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) Project-ID 496854903 (SA4064/2-1, PI Sander) and is a member of the International Max Planck Research School for Astronomy and Cosmic Physics at the University of Heidelberg (IMPRS-HD). DP acknowledges financial support by the Deutsches Zentrum für Luft und Raumfahrt (DLR) grant FKZ 50 OR 2005. JSV gratefully acknowledges support from STFC via grant ST/V000233/1. RK acknowledges financial support via the Heisenberg Research Grant funded by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) under grant no. KU 2849/9, project no. 445783058. We thank Danny Lennon and Philip Dufton for providing additional spectra.

Appendix A Additional observational data

In Table A.2, we list information on the observed spectra. We provide the observation dates for each BSG spectrum. For the UV, this information is retrieved via the ULLYSES search form and for the optical data this information is listed as output from search on ESO-archive portal. For the comments on individual stars provided in appendix D, we further used additional spectra for which analogous information is given in Table A.1.

Appendix B Atomic data

In Table B.1 we list the ions as well as the number of levels, superlevels, and transitions used in our models for the earlier (E) or later (L) spectral types of the BSG sample. The split between these “groups” is set at 21 kK, roughly corresponding to the B1 spectral type.

Appendix C Rotation and macroturbulence

In Table C.1, we list the complete set of broadening velocities for our sample stars obtained with the IACOB Broad tool for the He I λ4713 and Si III λ4552 (or Mg II λ4481 for later BSGs) lines.

For AzV 187, the minimum value derived by the tool is identical to the lower limit that can be obtained. Therefore, we instead use a more realistic value obtained by convolving the observed spectrum. We indicate this in Table C.1 by providing the automatic and the finally adopted value separated by a double-lined arrow (“⇒”). Likewise, the uncertainties Δv sin i and Δv_mac can be severely underestimated if both diagnostics yield a very similar value. In such cases, we increase the adopted uncertainty to the average of the standard deviation of the overall sample. This is indicated in Table C.1 with a single-lined arrow (→) where the values on the left are the derived standard deviations between the He and metal lines, and the new uncertainty is given on the right. However, in some cases, we notice that Δv_fft is larger than the updated uncertainties in v sin i. Hence, we further update Δv sin i in those cases yielding the final values given to the right of the half-arrows (“⇀”). The final values for the broadening values and their respective uncertainties employed in our analysis are written in boldface.

In Fig. C.2 and Fig. C.1, we show the output of the IACOB Broad tool for obtaining the v sin i and v_mac of AzV 18. The first figure depicts the application to Si III λ4552 while the second figure shows the application to He I λ4713. In this case, the values from both lines are in good agreement. However, this is not the case for all targets, as for example shown in Fig. C.3 where we depict the results for the Mg II λ4481 line in AzV 187. For this object, the best-fit v sin i is the minimum value (10.7 km s⁻¹) and the 1σ and 2σ intervals (i.e., the uncertainty) in the goodness-offit plot are very large. In such situations, we disregard the derived v sin i value and instead test different values by convolving the synthetic spectrum. Still, we keep the v_mac provided by IACOB Broad tool for the metal line $v_{\text{mac }}^{\text{Metal }}$. Such procedure is illustrated in Fig. C.4.

Fig. C.1 Output of IACOB Broad tool for AzV 18 applied to Si III λ4552.

Fig. C.2 Output of IACOB Broad tool for AzV 18 applied to He I λ4713.

Fig. C.3 Output of IACOB Broad tool for AzV 187 applied to Mg II λ4481.

Fig. C.4 Model spectra of AzV 187 with different rotational broadening illustrating the fitting process after the unsuccessful application of IACOB Broad tool for this target in particular. The macroturbulence is kept fixed (vmac = 48 km s−1) in all the spectra.

Appendix D Comments on Peculiar Targets

AzV 235 (B0.2 Ia). The peculiarity of this target is the presence of an unusually strong Hα, which was reported by Evans et al. (2004a). Compared with the old archival UVES spectrum, we see that the strong Hα is still there. However, the data for 2001 seems to have an extra increase on the blue part of the line, which is mirrored by the Hβ profile. See Fig. D.1.

In our modeling we were not able to have the intensity of Hα and Hβ compatible with the observed spectrum even by applying very different clumping settings (i.e., very low f_∞ = 0.02), enforcing clumping to start deeper in the wind and using a clumping-profile law that produces a smooth wind at terminal speed (Najarro et al. 2011). The exploration of other clumping settings was also motivated by the presence of an unusual behavior at the IR spectrum. For example, we noticed that the Brackett series is also in strong emission with possible double-peaked structure (see Fig. D.2). This double-peaked feature appears to happen with He I λ10830 as well.

In general, double peak profiles are sign of disk emission, which could also explain the presence of the super strong Hα emission. Moreover, we also notice a slight IR excess in this star (see first panel of Fig. F.1 in the appendix F available at Zenodo), which can also be due simply to higher wind-density via free-free emission. Another possible explanation for the Hα is the presence of a circumstellar nebula. However, there is to date no clear signs of nebula region around this star.

Fig. D.1 Time variation between 2001 (dark blue line) and 2020 (red line) of the spectrum of AzV 235 in the regions of Hα and Hβ.

AzV 104 (B0.5 Ia). AzV 104 is found to be a peculiar target because its UV P Cygni profiles are considerably narrow, indicating a very low wind speed. Additionally, the profiles of Si IV and C IV, which usually are well developed in stars as early as B0.5, do not have emission components. We could only partially reproduce the lack of emission component by using a very high wind microturbulence velocity. Evans et al. (2004b) and Parsons et al. (2024) also had problems in fitting the profiles of Si IV using SEI.

Fig. D.2 Observed (black line and dots) and model (red line) spectrum of AzV 235 around He I λ10830 (upper panel) and Part of the Bracket series (remaining panels).

In Fig D.3 we show the corresponding profiles and N V λ1240 for AzV 104 and AzV 488 as a comparison with a “typical” BSG of the same spectral type. The N V λ1240 which also seems pretty weak in AzV 104 likely indicates a low X-ray luminosity, as we needed a much reduced L_X/L to match it.

Sk 191 (B1.5 Ia). For this star, we did not manage to obtain successful simultaneous fits for Hα and Hβ, which both appear as quite developed P Cygni profiles. BSGs that display P Cygni/emission profiles in their Balmer series (at least Hβ) are considered hypergiants (e.g., van Genderen et al. 1982; Lennon et al. 1992; Clark et al. 2012), and therefore Sk 191 could be classified as a BHG under this criterion.

Typically, when bluer Balmer lines present wind signatures, Hα flux is considerably high, as is the case of AzV 78. However, the Hα profile of Sk 191 seems rather weak as one would expect from its Hβ – the reason behind the difficulty for a simultaneous fit. The best compromise we obtained was by adopting a very high v_cl = 400 and a very low amount of clumping aiming to limit the strength of Hα. Additionally, other lines of the Balmer series, as well as certain He I lines, seem to display some wind signatures (either slightly displaced blueward or showing a P Cygni profile) – see Fig. D.4.

Fig. D.3 Observed (black line) and model (red line) spectrum of AzV 104 in the regions of N V λ1240, Si IV λ1400, and C IV λ1550. The UV spectrum of AzV 488 is shown as a gray dashed line for comparison.

Interestingly, the UV spectrum was much easier to fit, and even normally problematic lines as Si IV λ1400 and Al III λ1855 were very well reproduced. The variability in the UV is not extensive (i.e., without changes to the overall morphology of the line), as it is possible to see in the Si IV λ1400 in Fig. D.5. However, in the same figure, we can also notice that the width of the absorption trough of C IV λ1550 appears to vary about ~100 km s⁻¹, which might also be indicative of variability in the wind terminal speed, or important ionization changes in the outer wind.

Fig. D.4 Optical time variability of Sk 191. The figure shows the Balmer lines (upper and lower rows of plots) as well as He I λ4471 and Si III λ4552,67,78 in the middle row.

Fig. D.5 Time variability in the UV for Sk 191. The upper and lower panels shows, respectively, Si IV λ1400 and C IV λ1550. The blue curve is the ULLYSES spectrum.

AzV362 (B3 Ia). The spectra in the literature of this star show “typical” well-developed P Cygni Hα (Trundle et al. 2004; Evans et al. 2004c). This indicates a quite variable behavior. In the UV, one does not observe strong wind lines, in principle incompatible with the strength of the Hα – which could be due to clumping, though. Such UV-Hα incompatibility happens with AzV 22 as well, albeit there is nothing particularly unusual about its line spectrum. Both AzV 362 and AzV 22 are the only cool BSGs in the sample with a higher degree of clumping.

Fig. D.6 Time variation of Hβ and Hα of AzV 362. Both lines seem to show a prominent blueshifted emission in 2001 and 2020. The Hα in the recent X-Shooter observation even appears as a completely inverted P Cygni profile, which is also shown in Fig 2. Higher-order Balmer lines do not show significant variability. The C II lines are shown for alignment purposes.

Fig. D.7 Observed (black line and dots) and model (red line) spectrum of AzV 362 around He I λ10830.

Looking at the infrared X-Shooter data, we can also notice a (possibly) inverted P Cygni in He I λ10830, as shown in Fig. D.7. Variability (and even momentarily unusual shapes) of Hα among OBSGs is not something uncommon (e.g., Martins et al. 2015; Kraus et al. 2015), however, given the apparent persistence of it and the fact that it might be appearing in helium lines, it is possible that there is something special going on with this object.

The flux of the FUSE spectrum for this star is considerably higher than the HST UV spectra (cf. Fig. F.3 in the appendix F available at Zenodo). This is likely due to an incorrect flux calibration of the FUSE range as other parts of the SED are well captured by the model.

Appendix E Fitting diagnostic lines

In Fig. E.1, we provide an illustrative example for the spectroscopic fitting process to determine the stellar and wind properties. Different stellar atmosphere models with varying parameters are created with Fig. E.1 showing the observed spectrum compared to models with different values for T_eff and $\dot{M}$. Guided by the knowledge on how different lines are affected by variations of the atmosphere properties, we obtain the model which marks the best compromise in fitting the many diagnostic lines for the different stellar and wind properties.

Similarly to Fig. E.1, in Fig. E.2 we illustrate the process of obtaining the surface abundance of CNO for one of our targets.

Fig. E.1 Illustration of the determination of parameters of AzV 215. The solid gray lines are the observed spectra. The black filled line is the best-fit model. The blue dotted(red dashed) line represents a model with higher(lower) temperature. Likewise, the green dotted(purple dashed) line shows models with lower(higher) mass-loss rates.

Fig. E.2 Illustration of the determination of the CNO abundance determination for AzV 18. The observed spectrum is the gray line. The thick black line is the model we found to have the best compromise between the diagnostic lines for each element. Each thin line represents the best-fit abundances for the respective element multiplied or divided by factors of 2, 3, and 5. This also guided us to adopt an uncertainty of 0.3 dex.

References

All Tables

All Figures

	Fig. 1 Positions of BSGs in the SMC from ULLYSES. The filled diamond symbols indicate the BSGs analyzed in this study, whereas the empty diamonds indicate the remaining BSGs of the total ULLYSES sample.
In the text

	Fig. 2 UV and optical spectra used to derive stellar parameters of the sample stars. This figure illustrates the variation in some spectral features with the spectral type. From left to right, the columns show spectral regions around Si IV λ1400, C IV λ1550, O+N+Hδ+He I+Si lines, and Hα+C II.
In the text

Fig. 3 Comparison between observed spectra of AzV 210 (B1.5 Ia) vs. AzV 472 (B1.5 Ia), and AzV 22 (B3 Ia) vs. Sk 179 (B3 II). For the comparison between the B1.5 stars, their optical and photospheric spectra are very similar to each other, but their wind lines at the UV and Hα are considerably different, despite the similar luminosity and temperature. In the juxtaposition between B3 and B8, the profiles in their UV are similar, while their optical spectra differ, especially due to broadening and different surface gravities. In each comparison, we list the old and the updated spectral types (Crowther, priv. comm.).

In the text

	Fig. 4 Spectral fits of exemplary subsets of diagnostic lines for the sample BSGs. The black curves are the X-Shooter data for each star, whereas the colored curves are their respective CMFGEN models.
In the text

Fig. 5 Comparison between effective temperatures and spectral types. The red diamond symbols represent our spectroscopic determinations, whereas the hexagons are literature values from Evans et al. (2004a), Trundle et al. (2004), and Trundle & Lennon (2005). The plotted curves are the calibrations from Schootemeijer et al. (2021) for the SMC stars (solid blue curve), Searle et al. (2008, dotted green curve), and de Burgos et al. (2023, dashed dark green curve) for Galactic BSGs.

In the text

	Fig. 6 Comparison between the derived effective temperatures and luminosities of our sample SMC BSGs with results from the literature, namely: Evans et al. (2004a, fuchsia hexagons), Trundle et al. (2004, pink hexagons), Trundle & Lennon (2005, orange hexagons), and Schootemeijer et al. (2021, black squares).
In the text

Fig. 7 HR diagram of the sample stars (HRD). Tracks of Georgy et al. (2013) (G13, thick orange lines) and Brott et al. (2011) (B11, thin gray lines) with initial rotation of veq/vcrit ~ 0.4 (i.e., ~300 km s−1). The small triangular and circular points on the tracks mark intervals of about 50 kyr for G 13 and B 11 respectively, and the red diagonal lines represent the iso-radii lines (1, 10, and 100, R⊙ from left to right).

In the text

	Fig. 8 HR diagram showing helium enrichment (Y/Yinit, where Yinit is the baseline value of 0.27) in the B11 and G13 SMC tracks compared to our sample BSGs (diamonds). The tracks are the same as in Fig. 7.
In the text

Fig. 9 Comparison of effective temperature, rotation, and CNO enrichment–depletion between G13 (orange tracks with triangle symbols, Georgy et al. 2013), B11 (gray, blue, and pink tracks, Brott et al. 2011), and our spectroscopic determinations (white diamonds). In this figure, we include tracks with initial masses from 12 to 50 M⊙ (to 40 M⊙ for the G13 set). The tracks’ equatorial rotation velocities are scaled by π/4 to account for the stars’ inclinations (e.g., Hunter et al. 2009). As in Fig. 7, the points along the tracks mark intervals of about 50 kyr. Baseline CNO values of B11 models are respectively 7.37, 6.50, and 7.98, while for G13 those are 7.56, 6.95, and 7.83. The dashed horizontal lines mark the CNO average baselines listed in Vink et al. (2023, Table 2) determined as the average values from previous studies. Panel A and B show the evolution of the rotation against Teff while panels C, D, and E show respectively the evolution of the surface C, N, and O vs. the rotation.

In the text

Fig. 10 N/C vs. N/O diagram. The dashed lines show the CN and CNO equilibrium curves for the SMC (in fuchsia) and LMC (in black) using the respective baseline abundances from XShootU Paper I (Vink et al. 2023) and from Hunter et al. (2009), which was the same used by the LMC merger models of Menon et al. (2024, circles). The star symbols with green contours represent the ratios determined for SMC stripped stars from Pauli et al. (2022, six-pointed star), Rickard & Pauli (2023, seven-pointed star), and XShootU Paper VIII (Ramachandran et al. 2024, five-pointed star). The diamonds show the ratios of our sample BSGs and the fuchsia hexagons represent ratios determined by Evans et al. (2004a); Trundle et al. (2004); Trundle & Lennon (2005). Those do not provide He abundance determination, whose enrichment relative to the baseline is represented by the color scale applied to the symbols.

In the text

Fig. 11 Comparison between evolutionary and spectroscopic masses. The color bar represents the effective temperature of the targets with lighter colors reflecting lower Teff. The fuchsia hexagons correspond to the data of Trundle et al. (2004) and Trundle & Lennon (2005) and the dotted blue lines connect the same stars. The darker horizontal orange and gray bands represent, respectively, the intervals of 1σ obtained for G13 and B11. Likewise, the lighter colored bands are for 2σ.

In the text

	Fig. 12 Comparison between models of AzV 488 (B0.5 Ia) with and without X-rays. The green dashed curves show the model without X-rays and the green full curves are the model with X-rays included. For this model, the log(LX/L) −7.4 with a TX parameter of 1.0 MK. All the other properties are identical between each model.
In the text

	Fig. 13 Comparison between models of AzV 187 (B3 Ia) and AzV 343 (B8 Iab) with and without X-rays. The dashed curves show the model without X-rays and the full curves are the model with X-rays included.
In the text

Fig. 14 X-ray luminosity compared to effective temperature. The thick diamonds with numbers represent the SMC targets AzV (or Sk for 179 and 191), while the circles with letters, in alphabetical order, represent the Galactic targets ε Ori (B0 Ia, Puebla et al. 2016), κ Ori (B0.5 Ia, Huenemoerder et al. 2011; Haucke et al. 2018), 9 Cep (B2 Ib), 55 Cyg (B2.5 Ia), o2 CMa (B3 Ia), and 64 Oph (B5 Ib/II, Bernini-Peron et al. 2023, as for the latter three). The colors of the circle and broad diamond symbols represent the wind terminal velocity. The purple thin diamonds represent the alternative values considering a TX = 0.5 MK for the corresponding early target connected by gray dashed lines. Likewise, the fuchsia crosses represent the alternative value for AzV 104 if considering a TX = 1.0 MK.

In the text

Fig. 15 Comparison between models with different relative X-ray luminosity and shock temperature parameters. The thin gray line is the normalized UV spectrum and the colored lines are the CMFGEN models. The blue filled line is the model with TX = 0.5 MK and log(LX/L) = −7.4. The purple dashed line has the same temperature, but with a lower log(LX/L) = −8.0. The red dot-dashed line has the same log(LX/L) = −7.4, but a higher TX = 1.0 MK. The different LX are obtained by varying the X-ray emission filling factor. All the additional X-ray parameters are kept the same for all the models.

In the text

Fig. 16 Panel A: Comparison between derived terminal velocities (v∞) by this study and those from Hawcroft et al. (2024, small blue diamonds) and Parsons et al. (2024, green squares). The fuchsia hexagons are data from Evans et al. (2004a), which also obtained wind properties via optical and UV analysis with CMFGEN. Panel B: Comparison of the wind microturbulent velocities (ξmax) between this work and Hawcroft et al. (2024). Panel C: Comparison between the sum of v∞ and ξmax, which would roughly translate the width of the absorption component of the UV P Cygni profiles. Panel D: Comparison of the Γ-corrected escape velocities (vesc,Γ) between this study and Hawcroft et al. (2024).

In the text

Fig. 17 Terminal velocities as a function of effective temperature: The results from this work (red diamonds) are compared to the values from Hawcroft et al. (2024, blue circles). The thick solid red and solid blue lines show the resulting relation when considering only supergiants (BSGs and OBSGs). The thick red dotted line shows the relation for only BSGs, but excludes the peculiar AzV 104. The blue dot-dashed line shows the linear relation derived for all SMC OB stars (of all types and classes) in the Hawcroft et al. (2024) paper also considering the literature results. Conversely, the dark blue dashed line refers to the LMC.

In the text

Fig. 18 Escape-normalized terminal velocity vs. effective temperature. The dashed black lines are the relation presented by Kudritzki & Puls (2000) compiling previous data. The dot-dashed line is the relation presented by Krtička et al. (2021). The red and blue thick solid lines are the average value of our determinations and their respective shaded regions are the standard deviation. The thin dashed fuchsia lines are the Kudritzki & Puls (2000) relation scaled by (Z/Z⊙)0.1 metallicity dependence from Leitherer et al. (1992) and the thin dashed green line below it is the metallicity scaling found by Hawcroft et al. (2024). The brown horizontal lines indicate the velocity scaling proposed by Vink & Sander (2021).

In the text

Fig. 19 Mass-loss vs. luminosity (upper panel) and modified wind-momentum vs. luminosity (lower panel) diagrams of SMC OB stars of different sources and luminosity classes. The different symbols mark different studies and different colors/sizes mark different luminosity classes. Error bars are not shown for visual clarity. The literature values are slightly translucent whereas the diamonds represent our sample BSGs. The error bars in the lower right corner indicate the typical associated errors for our sample. The scatter points with downward black triangles denote upper limits.

In the text

Fig. 20 Upper panel: Mass-loss rates vs. effective temperature. The scatter points are color-coded by their Γe values. The three points outlined in fuchsia in each plot are the BSGs whose Balmer line profiles were peculiar, and thus have less reliable empirical values for the massloss rates. Lower panel: Transformed mass-loss rates (Gräfener & Vink 2013) vs. effective temperature. The color-coding is the same as in the upper panel.

In the text

Fig. 21 Upper panel: different mass-loss prescriptions applied to a generic star with the same properties but varying Teff for illustrative purposes, to facilitate the visualization of the predicted $\dot{M}$(Teff) behaviors. Vink et al. (2001) rates are represented by white circles, Krtička et al. (2024) by white squares, and Björklund et al. (2023) by white triangles. Middle panel: comparison between different mass-loss rate prescriptions and inferred values for the sample stars (black diamonds) and OSGs from Bouret et al. (2021, gray pentagons). The fuchsia-outlined diamonds are the targets where we did not find a satisfactory fitting of Balmer lines. The applied theoretical values have the same symbols as the upper panel. The stars are ordered according to their Teff in the upper axis. The lower axes show the spectral types of each star. Lower panel: similar to the previous panel, but showing the ratio of the predicted according to each recipe ($\dot{M}$th) to our empirical determinations. The symbols follow the same encoding from the previous panels.

In the text

	Fig. C.1 Output of IACOB Broad tool for AzV 18 applied to Si III λ4552.
In the text

	Fig. C.2 Output of IACOB Broad tool for AzV 18 applied to He I λ4713.
In the text

	Fig. C.3 Output of IACOB Broad tool for AzV 187 applied to Mg II λ4481.
In the text

	Fig. C.4 Model spectra of AzV 187 with different rotational broadening illustrating the fitting process after the unsuccessful application of IACOB Broad tool for this target in particular. The macroturbulence is kept fixed (vmac = 48 km s−1) in all the spectra.
In the text

	Fig. D.1 Time variation between 2001 (dark blue line) and 2020 (red line) of the spectrum of AzV 235 in the regions of Hα and Hβ.
In the text

	Fig. D.2 Observed (black line and dots) and model (red line) spectrum of AzV 235 around He I λ10830 (upper panel) and Part of the Bracket series (remaining panels).
In the text

	Fig. D.3 Observed (black line) and model (red line) spectrum of AzV 104 in the regions of N V λ1240, Si IV λ1400, and C IV λ1550. The UV spectrum of AzV 488 is shown as a gray dashed line for comparison.
In the text

	Fig. D.4 Optical time variability of Sk 191. The figure shows the Balmer lines (upper and lower rows of plots) as well as He I λ4471 and Si III λ4552,67,78 in the middle row.
In the text

	Fig. D.5 Time variability in the UV for Sk 191. The upper and lower panels shows, respectively, Si IV λ1400 and C IV λ1550. The blue curve is the ULLYSES spectrum.
In the text

	Fig. D.6 Time variation of Hβ and Hα of AzV 362. Both lines seem to show a prominent blueshifted emission in 2001 and 2020. The Hα in the recent X-Shooter observation even appears as a completely inverted P Cygni profile, which is also shown in Fig 2. Higher-order Balmer lines do not show significant variability. The C II lines are shown for alignment purposes.
In the text

	Fig. D.7 Observed (black line and dots) and model (red line) spectrum of AzV 362 around He I λ10830.
In the text

	Fig. E.1 Illustration of the determination of parameters of AzV 215. The solid gray lines are the observed spectra. The black filled line is the best-fit model. The blue dotted(red dashed) line represents a model with higher(lower) temperature. Likewise, the green dotted(purple dashed) line shows models with lower(higher) mass-loss rates.
In the text

Fig. E.2 Illustration of the determination of the CNO abundance determination for AzV 18. The observed spectrum is the gray line. The thick black line is the model we found to have the best compromise between the diagnostic lines for each element. Each thin line represents the best-fit abundances for the respective element multiplied or divided by factors of 2, 3, and 5. This also guided us to adopt an uncertainty of 0.3 dex.

In the text