© The Authors 2023

1 Introduction

Understanding massive stars (M_ini ≳ 8 M_⊙) is paramount to comprehend the visible Universe. Throughout the cosmic history, these objects inject energy, momentum, and processed material into the interstellar medium (ISM) of their host galaxies, both through their powerful stellar winds and their violent deaths, for example, as core-collapse supernovae (e.g., Smith 2014). The evolution of massive stars not only determines their final fate, but also the impact they have on their surrounding environment. Therefore, the quantitative study of the properties and behaviors of massive stars in different stages gives us important constraints and insights on their still quite uncertain evolution -especially regarding post-main-sequence evolution (e.g., Martins & Palacios 2013).

B-type supergiants stars (BSGs) are particularly interesting as they mark an evolved stage of massive stars beyond the zero-age main sequence (ZAMS) where both core-H and core-He burning objects can be found (e.g., Kraus et al. 2015; Martin et al. 2018). Moreover, stars in the BSG regime may evolve blue-wards or redwards (Georgy et al. 2014; Haucke et al. 2018). Determining the parameters of BSGs could be useful in distinguishing between these different stages and investigating the connection between OB main-sequence stars, BSGs, and other evolved objects such as yellow and red supergiants and hyper-giants, Wolf-Rayet stars, and luminous blue variables (Crowther et al. 2006; Clark et al. 2012).

The winds of BSGs are thought to be radiatively driven, similar to their hotter O-star counterparts. However, particular phenomena inherent to their parameter space, such as an observed sharp decrease in the terminal velocities (Lamers et al. 1995) and superionization (e.g., Walborn & Nichols-Bohlin 1987), turned out to be challenges for their theoretical understanding. Following early simulations by Pauldrach & Puls (1990) that revealed a bi-stability of solutions for the B1 hyper-giant P Cyg, the term "bi-stability jump" has been introduced to refer to the sudden decrease in υ_∞. Subsequently, Monte Carlo simulations by Vink et al. (1999) found a corresponding increase in the mass-loss rate¹, giving this jump an important evolutionary relevance. However, newer modeling efforts yield conflicting results (cf. Krtička et al. 2021; Björklund et al. 2022), challenging the paradigm of a mass-loss increase for cooler B stars.

Using the stellar atmosphere code CMFGEN (Hillier & Miller 1998), Crowther et al. (2006) and Searle et al. (2008) modeled in detail the UV and optical spectra of 31 galactic supergiants, covering the spectral types B0 to B5. Although they managed to obtain physical properties for their sample, they reported problems in reproducing important UV profiles (e.g., N V, C IV, and Al III). Both studies did not include wind clumping and X-ray emission in their modeling, thereby leaving behind open questions about the actual wind structure and some derived properties.

More recent papers, such as those of Martins et al. (2015) and Puebla et al. (2016), analyzed BSGs including the aforementioned effects and were able to achieve better spectral fitting in the UV. Nevertheless, those papers studied only B0- and B0.5-type supergiants, which seem to adjust to the expected behavior of late-O supergiants (OSGs) with regard to their clumping and X-ray properties, namely: highly inhomogeneous winds and high relative X-ray luminosity (L_x/L), (Driessen et al. 2019; Berghoefer et al. 1997).

For these stars, the value of log(L_x/L) ~ −7 is typically employed in their modeling (Bouret et al. 2012; Puebla et al. 2016, e.g.). These values are within the observational constraints of Berghoefer et al. (1997). However, for BSGs later than B1, Berghoefer et al. reports no detection of emitted X-rays, having only upper limits available with significant scattering. For O supergiants, Nebot Gómez-Morán & Oskinova (2018) found a large spread of log(L_x/L), from −5 to −7.5, raising doubt whether the value of −7 is universally applicable to supergiants. Moreover, Nazé (2009) reported X-ray detections for some cool BSGs with log(L_x/L) > −7, albeit using only bolo-metric corrections based on spectral type. Nonetheless, these different findings reveal that hot supergiants in general still lack constraints on their X-ray emission properties.

In this work, we perform detailed atmosphere modeling of cool BSGs (later than B1) including both clumping and X-rays. Using both high-resolution UV and optical spectra, we analyze a small sample of BSGs in detail, namely: HD206165 (B2Ib), HD198478 (B2.5Ia), HD53138 (B3Ia), and HD164353 (B5Ib/II). In our study, we employed the CMFGEN atmosphere code to obtain the stellar and wind parameters. Additionally, motivated by Prinja & Massa (2010) and Petrov et al. (2014), we also model HD53138 using PoWR (Gräfener et al. 2002; Hamann & Gräfener 2003; Oskinova et al. 2007) to include an approximate treatment of optically thick clumps, and test whether such "macroclumping" can improve the agreement of lines such as Si IV λ1394–1403 and Hα.

The rest of the paper is structured as follows: Sects. 2 and 3 present our methodology, first explaining how we obtained and treated our observational data and then detailing how we modeled the atmosphere and the output spectra in order to obtain the relevant physical parameters. In Sect. 4 and 5 we present our results, the former focusing on the impact of clumping and X-rays in our sample stars, and the latter as an additional discussion of their evolutionary status. As closure, we draw our conclusions in Sect. 6.

2 Observations and data reduction

For our study, we used publicly available data of a sample of cool BSGs, presented in Table 1. We used (i) high-resolution UV and optical spectra to obtain the wind and photospheric physical properties; (ii) magnitudes from UV to IR; and (iii) parallax distances (more details below) to estimate the stellar luminosities by reproducing the spectral energy distribution (SED).

The UV spectra were taken with IUE, using the long-wavelength (LWP) and short-wavelength (SWP) spectrographs, with the R ~ 0.2 Å resolution (echelle grating and cross-disperser). This is the region where most of the information on the wind is encoded for OB stars. We also used FUSE spectra available from the MAST database to better constrain the SED in the UV.

In order to have a better signal-to-noise ratio and clearer line profiles in the UV we computed an average spectrum built from all available IUE spectra for our sample stars with high resolution and large aperture (i.e., spectra with HIGH and LARGE flags at the MAST archive). When not available, we used the spectra with the small aperture setting (spectra flagged with SMALL), as the absolute flux calibration is not important for the line profiles evaluation. Albeit this blurs any imprints of variability, the broad structures of the overall wind profiles (which are the focus of our study) do not change drastically over time (e.g. saturating and de-saturating or changing width noticeably), as we can see in Fig. 1. The star with the highest observed variability is HD164353, which may suggest that the BSG was subjected to some variation on the wind properties². Nevertheless, such behavior does not negatively influence the focus of this study, namely, the wind analysis and the obtained general stellar properties. The spectra were manually normalized by fitting cubic splines in order to remove any spurious oscillations. For computing the SED, we used the LWP spectra as well as large-aperture spectra, regardless of the resolution. Table 2 lists the file number for each of the spectra obtained from the MAST archive.

The optical spectra were acquired from the POLARBASE database³ (Donati et al. 1997; Petit et al. 2014) - sourcing from ESPaDOnS (R = λ/Δλ ~ 70000) and NARVAL (R ~ 65000) – and the ESO Archive – using HARPS (R ~ 120000) and FEROS (R ~ 48000) spectra. For each star and instrument, we selected the spectra with the highest signal-to-noise ratio. These spectra, which practically cover the whole optical region (3900 to 7000 Å), were normalized with the same procedure as the UV spectra.

The BSGs in our sample show consistent spectra over time and a "well behaved" SED⁴. This means that the photospheres of these stars are not subject to significant large-scale variability that would spoil determinations of the main properties by modeling their spectra. The only major variability in the spectra of these stars are seen in the Ha profiles, which significantly vary over the years as illustrated in Fig. 2. This variability is more pronounced in HD198478 (B2.5Ia) and HD53138 (B3Ia).

For HD198478, Kraus et al. (2015) also tracked the variations of the Hα profile for about one month and showed that its variability happens quite quickly. They concluded that this is unlikely to reflect sharp mass-loss rate variations. Instead, Kraus et al. (2015) suggested large-scale inhomogeneities propagating through the wind, reinforcing the importance of including clumping in the analysis of these objects. For HD53138, the Hα variability was studied by Morel et al. (2004), who reported no periodicity in the variability behavior. Martin et al. (2018) pointed out this variability could be due to the presence of magnetic fields, but to date, the literature offers no evidence for their presence in the star's atmosphere.

The photometric magnitudes were obtained from multiple sources: FUV and NUV are taken from GALEX (Beitia-Antero & Gomez de Castro 2016), UBVI from the XHIP catalog (Anderson & Francis 2012) and R is obtained from Morel & Magnenat (1978); Gaia's Gbp, G, and Grp were taken from Gaia eDR3 (Fabricius et al. 2021); JHK were obtained from 2MASS (Cutri et al. 2003); and WISE s W1 and W2 were obtained from the AllWISE catalog (Cutri et al. 2013).

For HD198478 and HD206165, we use the distances from Bailer-Jones et al. (2021) employing Gaia eDR3 parallaxes. The distances for HD53138 and HD164353 were derived from the revised HIPPARCOS parallax (van Leeuwen 2007). We opted not to use Gaia s parallax for these stars because they are considerably bright (V = 2.94 and 3.96 respectively), which is close to the allowed limit of Gaia (Lindegren et al. 2018). As Golovin et al. (2023, Sect. 3) also discusses, for such bright stars the HIPPARCOS parallaxes tend to be more reliable. Moreover, the "renormalized unit weighted error" (ruwe) for these stars amounted to 3.62 and 2.51, values that are much greater than the cut-off of 1.4 suggested by Fabricius et al. (2021).

Fig. 1 UV spectra of our sample stars – CIV 1550 and Si IV 1400 region. Note: variability is absent or mild in the available IUE data as the profiles do not present drastic changes (e.g. appearance/destruction of P-Cygni profile or variations in the width). The black profiles are the computed average IUE spectra which were used for the model comparison. For HD164353, the green lines represent small aperture IUE data (re-scaled to match the flux of the large aperture acquisition). The black thick lines are the averaged spectra which were used for model comparison (see Sect. 3).

Fig. 2 Time sequences of Ha profiles for our sample BSGs. The high variability in HD53138 and HD198478 makes the profile not very suitable for obtaining reliable wind properties without auxiliary diagnostics.

3 Stellar atmosphere models

In order to model the BSG atmospheres and analyze their spectra, we used the stellar atmosphere codes CMFGEN (Hillier & Miller 1998; Hillier et al. 2003) and PoWR (Gräfener et al. 2002; Hamann & Gräfener 2003; Sander et al. 2015). Both codes compute the radiative transfer in a co-moving frame formalism and solve the statistical (non-LTE) and radiative equilibrium equations in an expanding atmosphere assuming spherical symmetry and stationarity. They also account for the line-blanketing effect arising from the millions of line transitions in heavier elements such as iron. Once the atmospheric stratification is converged, the synthetic spectrum is computed via a formal solution of the radiative transfer in the observer's frame (see Mihalas 1978, for details). While the codes have very different architectures, which might produce small variations of inferred parameters, the underlying physics within each are similar and throughout the past years were used to successfully model spectra of hot stars (cf. Oskinova et al. 2007; Bouret et al. 2012, for a comparison in modeling ζ Pup as an example).

In this work, we did not solve the equation of motion consistently as we aim to obtain empirically derived parameters for the wind properties. Instead, as typically done in the literature, a fixed wind velocity structure was adopted in the form of a (modified) β-law (see, e.g., Bouret et al. 2012), which smoothly connects the wind regime with the quasi-hydrostatic inner layers. The connection velocity in our models is set to 10 km s⁻¹, slightly lower than the sound speed.

The winds of hot massive stars are line-driven and inherently unstable (Lucy & Solomon 1970). Therefore, time-dependent models predict contrasts in terms of density, velocity, and temperature that could give rise to X-ray emission in the outflow (see, e.g., Owocki et al. 1988). These effects can have an impact on the emerging spectra and have been extensively studied over the last decades (see, e.g., Bouret et al. 2005; Sundqvist et al. 2014, and references therein). In the following sections, we present how the effects of clumping and X-rays are taken into account in the CMFGEN and PoWR models.

3.1 Clumping and X-Rays in CMFGEN

CMFGEN is able to account for wind inhomogeneities via the "microclumping" formalism, where the wind is considered to consist only of optically thin clumps with a void interclump medium. In addition, the clump sizes are assumed to be smaller than the photons' mean free path. The volume filling factor is defined as

$$f_{\text{V}}\left(r\right)=f_{\text{V,}\infty }+\left(1-f_{\text{V,}\infty }\right)e^{-\upsilon \left(r\right)/{\upsilon }_{\text{cl}}},$$(1)

where f_v,∞ is a free parameter denoting the value at r → ∞. υ(r) is the velocity at radius r and υ_cl, the second free parameter, is a characteristic velocity describing the onset of clumping (i.e., where it starts to become relevant). A newer version of CMFGEN presented by Flores & Hillier (2021) allows the inclusion of clumping with arbitrary optical depths, which provides a description with fewer assumptions on the inhomogeneities. However, only a few objects had been analyzed in such a way and this version is not currently public.

X-ray emission can be included in CMFGEN as a product of thermal Bremsstrahlung following a similar approach as Baum et al. (1992) and Pauldrach et al. (1994). This formalism follows from the paradigm that X-rays are produced by shocks within the wind, aiming to parameterize this type of emission. As Hillier & Miller (1998) explains, the X-ray is added during the calculation of the opacities as an additional emissivity term. In the version used in this paper (2017 release) the X-ray emissivities are further making use of the APEC table (Smith et al. 2001) for solar metallicity to account for X-ray line opacities⁵. In the current description, three parameters describe the wind-intrinsic X-rays, namely: (i) the X-ray filling factor f_x, which can be interpreted as the amount of shocks producing X-rays in the wind, (ii) T_x, which represent the temperature of the shocks, defining the X-ray distribution in frequency space, and (iii) the X-ray characteristic velocity υ_x, which determines the spatial distribution of the X-ray emission. From this parametrization, the code computes the wind emitted and observable X-ray luminosity L_x (with energies >0.1 keV) as an output.

3.2 Clumping and X-Rays in PoWR

Similarly to CMFGEN, PoWR allows the inclusion of clumping with the same description. In addition, PoWR has the option to include an approximated formalism for optically thick clumping in the formal integral (Oskinova et al. 2007). In this so-called "macroclumping" formalism, the effective opacity, K_eff, is written as:

$${\kappa }_{\text{eff}}\left(r\right)=\kappa \left(r\right)\left(\frac{1-e^{-{\tau }_C\left(r\right)}}{{\tau }_C\left(r\right)}\right),$$(2)

where κ is the opacity considering microclumping only. τ_C (r) describes the optical depth of the clumps and is defined as

$${\tau }_C\left(r\right)=\kappa \left(r\right)L_0{\left(\frac{r^2}{R_{*}^2{f}_{\text{V}}^2\left(r\right)}\frac{\upsilon \left(r\right)}{{\upsilon }_{\infty }}\right)}^{1/3},$$(3)

where L₀ denotes the typical separation between the clumps.

Concerning shock-generated X-rays parametrization, PoWR includes these as emissions from thermal Bremsstrahlung as described in Baum et al. (1992). While this is similar to CMFGEN, the main difference is that PoWR so far does not include any emissivity tables - see Appendix D for more details on the difference in implementation between codes. Consequently, the effect of X-rays is treated as additional continuum emissivities only. However, this difference is not important for the present work, as we employ the PoWR models only to study the impact of macroclumping on specific line profiles (see Sect. 4.2).

3.3 Atomic Physics

For the model atoms in our CMFGEN models, we included 16 different elements with a total of 46 different ions. This marks a considerable extension compared with Crowther et al. (2006) and Searle et al. (2008), who used 11 elements and 27 ions. Moreover, we extended the number of so-called "superlevels", a concept introduced by Anderson (1989, 1991) to combine several similar-energy levels into a single "superlevel" and reduce the number of equations that need to be solved for the statistical equilibrium. We summarize our underlying atomic data in Table 3. A more detailed description of the atomic data and their implementation in CMFGEN can be obtained from the CMFGEN manual⁶ as well as Hillier & Miller (1998) and Dessart & Hillier (2010).

3.4 Stellar parameter determination

In order to obtain the physical properties, we aim to reproduce all the diagnostic lines in the UV and optical regimes simultaneously. This procedure consists of the following major steps:

Effective temperature and surface gravity. We determine the effective temperature T_eff using the ionization balance of Si II λ4128-30, Si III λλ4552-68-75, as well as He I λ4471 and Mg II λ4481. Given that Si IV λλ4089,4116 appears weakly for most of the targets we also use them as auxiliary T_eff diagnostics. For the surface gravity log g, we use the wings of Hβ, Hγ, and Hδ. The typical errors are 1 kK and 0.1 dex following a conservative estimate – similarly to Crowther et al. (2006) and Searle et al. (2008) – and to account for uncertainties due to microturbulence, which affects most of the optical lines, especially Si III λλ4552-68-75.

CNO abundances. Once these main properties are established, we derive the C, N, and O abundances using mainly C II λ4267, N II λ4447, N II λ/4601–4630 and O II, λ4069, O II λ4590–96,4661. The typical error is 0.3 dex, as we did not perform any fine-tuning – in contrast, for example, to what Puebla et al. (2016) did for HD37128 (B0Ia). The other elemental abundances are assumed to be solar (Asplund et al. 2009) and the He/H number fraction is assumed to be 0.1. This value represents a slightly enriched abundance (corresponding to a mass fraction of X_He = 0.28) since it is expected that BSGs are evolved massive stars with some surface He enrichment (Maeder 2009). Determining the abundance of He is not trivial and, according to Crowther et al. (2006), would not produce significant differences in the lines.

Line broadening. The rotational and macroturbulent broadening for our sample, υ sin i and υ_mac respectively, are taken from the IACOB survey (Simón-Díaz et al. 2017)⁷. The values are shown in Table 4.

Luminosity and extinction. The luminosity, L, is then obtained by reproducing the SED from the UV to the infrared, directly using the flux-calibrated IUE spectra and the flux derived from the available photometry and applying the corresponding distances and the reddening, which is obtained by fitting the 2175 Å absorption bump. The extinction was determined following the reddening law from Cardelli et al. (1989) with the values obtained for each star shown in Fig. E.1. For HD53138 and HD164353, whose distances came from HIPPARCOS (van Leeuwen 2007), we consider an error in L of 35%, which is the largest error in the distance for this part of the sample. Since these two objects have low reddening, the uncertainties associated with E(B – V) are negligible. Conversely, for HD198478 and HD206165, whose distances are sourced in Gaia eDR3, we have an error of ≲ 12%. However, as these objects are subject to higher extinction (E(B – V) ~ 0.5), we consider a L uncertainty of 20% in order to account for that.

Masses and radii. Using L, log g, and T_eff we can obtain the masses (M) and stellar radii (R). We considered that the uncertainties in the three aforementioned parameters follow a Gaussian distribution centered on the determined values, with their error estimates as the standard deviations (see Appendix A for details). From that, we computed the mass probability distribution for each star. To estimate the uncertainty in mass, we take the standard deviations of these resulting distributions. We obtained an error of 40% in the masses for the stars whose distances are based on Gaia parallaxes and 50% for those using HIPPARCOS data. Concerning the radii, we obtained an error of 15% and 20% respectively.

3.5 Wind parameter determination

After determining the photospheric parameters, we then obtained the wind properties, using mainly the UV diagnostics:

Mass-loss rates. We infer the clumping-corrected mass-loss rate ($\dot{M}/\sqrt{f_{\text{V},\infty }}$) from the Ha strength and the "main" UV P-Cygni profiles, namely: N V λ1238-42, C II λ1335-36, Si IV λ1394-1403, C IV λ1548-50 and Al III λ1855-63. Due to the variability in Hα, we did not attempt to fit this line precisely. Instead, we focus on obtaining models with compatible Hα strengths, namely, allowing flux discrepancies not greater than 10% relative to the observations (see Fig. 3 for an illustrative example).

Thus, an error of 30% in Ṁ was adopted based on how the synthetic spectra start to clearly deviate from the observations.

Wind velocity. The velocity gradient parameter β is fixed at 2.0 based on the value obtained by Crowther et al. (2006) and on the fact that according to Crowther et al. and Searle et al. (2008) a higher β is preferred in the case of cooler BSGs⁸. Moreover, a more recent study on BSGs conducted by Haucke et al. (2018), who analyzed the optical spectra of 19 stars using FASTWIND (Puls et al. 2005) and found β = 2.0 for HD53138 and a systematic trend of β > 2.0 for other cool BSGs. The terminal velocities were initially taken from Howarth et al. (1997). However, in order to properly fit the main UV lines' absorption widths, we adjusted them for each star if necessary. The adopted error associated with υ_∞ is 100 km s⁻¹.

Microturbulence scaling. The microturbulence in the photosphere ξ_min was constrained using the depth of Si III λλ4552-68-75 and photospheric UV absorption lines. In CMFGEN, the microturbulence ξ_turb in the formal integral is assumed to grow monotonically outward until it reaches its maximum value ξ_max, which is assumed to be 0.1υ_∞ as commonly done in the literature, following

$${\xi }_{\text{turb}}={\xi }_{\text{min}}+\frac{{\xi }_{\text{max}}-{\xi }_{\min }}{{\upsilon }_{\infty }}.$$(4)

Wind inhomogenieties. Regarding clumping, the characteristic velocity υ_cl is fixed to 30 km s⁻¹ while assuming that clumping starts already at the innermost layers of the wind – for all models, as it is commonly done in OB star analyses (e.g., Marcolino et al. 2009; Puebla et al. 2016). Moreover, Torrejón et al. (2015) analyzed X-ray emission in QV Nor, a B0I + neutron star (NS) binary, and provided empirical evidence that clumping starts at r < 1.25R», namely, very close to the photosphere. This would correspond to low wind velocities, in line with our assumption. For estimating the volume filing factor (f_v,∞), we tested models with different values of f_v,∞ = 0.1, 0.2, 0.5, and 1.0, the lattermost corresponding to a smooth wind, always keeping $\dot{M}/\sqrt{f_{\text{V},\infty }}$ constant.

X-rays. For the X-ray parameters we keep the characteristic velocity fixed at around 0.7 υ_∞. This choice is motivated by the fact that shock-induced X-rays would require a high dispersion of velocities in the wind and thus could only emerge further away from the star where the wind achieved higher speeds (~400 km s⁻¹). In O- and early BSGs, such values can be reached closer to the photosphere, ~ 0.2υ_∞ (Lagae et al. 2021). In Appendix D we discuss how this manifests in the spatial distribution of the X-ray emission of our models. This is also in line with expectations from hydrodynamical simulations by Driessen et al. (2019).

Once υ_x is fixed, we test how different combinations of T_x and L_x lead to a good enough fitting of the UV wind lines. Using HD53138 (B3 Ia) as the primary target, we performed a detailed investigation of the possible X-ray parameters before testing the validity of the results on the other stars of the sample.

Fig. 3 Comparison between synthetic and observational spectrum of HD53138 in the Ha region. The color scheme is the same as in Fig. 2 while the black filled line represents the preferred model while the and dotted lines represent models with an increase/decrease in M by a factor of 30%.

Fig. 4 Optical spectrum of HD53138 (black solid line) compared to our best-fit model (red-solid line). Here we depict the main diagnostic lines used to obtain the stellar photospheric properties.

4 Wind clumping and X-ray results

By fitting the observed UV and optical spectra of the sample stars, we obtained the main stellar properties. Using HD53138 as a showcase, we depict the best fit of the main (photospheric) diagnostic lines in Fig. 4. From the figure, we can notice that a satisfactory fitting is achieved for most of the lines, especially the diagnostic features.

The main stellar properties of our sample BSGs are presented in Table 4 and the results of the spectral fitting can be found in Appendix E. Additional information about the individual BSGs can be found in Appendix C. In Table 5, we compare our results for HD53138 with various studies from the literature, generally finding a good agreement despite the partially different underlying methodologies. A more in-depth discussion about the evolutionary status and consequences of the stellar properties obtained is undertaken in Sect. 5.

After the determination of the photospheric properties, we analyze the winds via the UV spectra of our sample stars. A particular emphasis here is to reproduce the so-called "superionization", a phenomenon denoting the presence of highly ionized elements which are not expected from the temperature regime reflected by T_eff (Cassinelli & Olson 1979).

In previous works by Crowther et al. (2006) and Searle et al. (2008), the biggest problem in reproducing the UV spectrum occurred due to the presence of both C II 21335-36 and C IV λ21548-50 with the former typically showing an overly strong Ρ Cygni profile while the latter did hardly show any wind signature.

We illustrate this in Fig 5, where we depict a model for HD53138 without clumping and X-rays, reproducing the findings of Searle et al. and Crowther et al.. In the same figure, we also show that clumping alone is not enough to address the observed superionization.

In Fig. B.1, we show our new results for HD53138, where we tested different combinations of clumping and X-ray parameters. As illustrated, we achieve a good agreement for the two carbon as well as most of the other UV lines, even when varying some of the X-ray and clumping parameters. The overall good agreement in Fig. 6 has already indicated that, despite the general lack of clear X-ray detections, the inclusion of X-rays in modeling cool BSG atmospheres appears to be vital. A more detailed view is shown in Appendix B.

However, with the present X-ray implementation, there is no clear preference for a particular X-ray parametrization, as the fit quality for all lines is somewhat similar. We therefore present and take into account different solutions obtained in our analysis instead of presenting only one, singular solution. Among those, we assign formally one of them as preferred to guide some of our later discussions. Concerning Si IV λ21394-1403, the model profiles are always stronger than observed and even saturated in the models. Given that Crowther et al. (2006) and Searle et al. (2008) obtained weaker lines than in the observations, we can conclude that this line is highly affected already by a moderate amount of X-rays: For every T_x and f_v,∞ tested, Si IV 21394-1403 becomes a strong P-Cygni profile. To avoid saturation, the resulting L_x would need to be reduced so much, that it would spoil the fitting of C IV λ21548-50.

For Al III λ21855-63, we obtain a similar overprediction, but with an opposite trend. The profile only de-saturates at higher L_x. This is in line with the overprediction of Al III λ21855-63 reported by Crowther et al. and Searle et al. when using models without X-rays.

Due to the conflicting behavior of Si IV and Al III, we put more emphasis on the reproduction of the C profiles. A more in-depth discussion of the Si lines in the context of optically thick clumping is given in Sect. 4.2. Exhibiting the broadest disagreement in previous studies, the carbon profiles provide a better idea of the overall ionization structure of the cool BSG winds as both profiles belong to the same element. Therefore, we consider a model as "good" if it recovers both lines simultaneously as unsaturated (albeit not too weak) P Cygni profiles.

Another line to consider is N V λ1238-42: Because of its very weak profile, we also considered models which do not predict N V at all, but still reproduce both C lines adequately, as potential solutions. Whether or not N V λ1238-42 is taken into account has direct consequences on the predicted X-ray properties for these stars.

4.1 Limits of the X-ray luminosity

For HD53138, we are able to find a simultaneous fit of C II λ1335-36 and CIV λ1548-50 for models with log(L_x/L) ~ −13 and T_x = 0.1 MK (largely independent⁹ of the adopted clumping). Such a value of T_x, representing a shock temperature (T_shock), is roughly aligned with the relation

$$T_{\text{shock}}\sim {\left(\frac{\text{Δ}\upsilon }{300\text{\hspace{0.17em}km\hspace{0.17em}}{\text{s}}^{-1}}\right)}^2\cdot {10}^6\text{K,}$$(5)

between T_shock and the wind velocity dispersion (Δυ) obtained by analysing O and early-B stars in Cohen et al. (2014) if we considers a wind velocity dispersion of ~70 km s⁻¹ (cf. the BSG model results by Driessen et al. 2019). However, none of the models with T_x = 0.1 MK (and the resulting very low L_x) were able to yield any trace of N V λ1238-42, which seems to appear to be a very weak profile in the observed spectrum of HD53138 and in the other sample BSGs as well. In short, this means that if we cannot detect N V λ1238-42 in the spectrum of cool BSGs, log(L_x/L) as low as ~ −13 cannot be discarded.

For T_x = 0.5 MK, both C lines as well as N V λ1238-42 can be reproduced simultaneously. In the solution with no clumping (f_V,∞ = 1.0), we require a lower log(L_x/L) ≈ −8.3, whilst for very clumped winds we obtain values up to log(L_x/L) ≈ −7.5 for f_v,∞ = 0.1.

Reasonable spectral fits are also obtained when adopting T_x = 1.0 MK. However, this temperature is hard to motivate as a shock temperature (e.g., Cohen et al. 2014; Driessen et al. 2019). For HD53138 with υ_∞ = 680 km s⁻¹, a velocity dispersion of ~ 300 km s⁻¹ would be required, which seems implausible. Later-type BSGs in general show even lower wind velocities. Thus, taking into account the wind-shock paradigm, our model set with T_x = 0.5 MK would be considered as the preferred description. The good fits with T_x = 1 MK would instead require a different theoretical justification. Among the T_x = 0.5 MK models, the solutions with f_V,∞ = 0.2 and 0.5 are scarcely different. For our later comparison discussions, we pick the 0.5 model as our preferred solution. The model yields a slightly stronger C II profile and thus is marginally better in comparison to the observations.

Combining all model efforts we can conclude that a minimum log(L_x/L) ≈ −8.3 is required to reproduce the UV spectrum of HD53138. This is in line with the current observational upper limit of −7.4 (Berghoefer et al. 1997).

HD198478. Similarly to HD206165, this star was best fitted by models with lower clumping, as more inhomogeneous wind models predict overly weak C lines. Furthermore, both stars have very strong Si IV, also largely recovered by our models. Still, one can notice a slight oversaturation, similarly to HD53138. The profile starts to de-saturate when f_V,∞ is decreased, but for the models with T_x = 0.5 MK, they immediately fail to produce a P-Cygni for C II and predict an overly strong Ν V. For the case with T_x = 1.0 MK, CIV is not fitted at all with lower f_v,∞. Thus, among the solutions, the model with f_v,∞ = 1.0 and T_x = 0.5 MK is the preferred option.

HD 164353. For our target with the latest spectral type, all models reproduce C II quite well, but only the models with lower clumping yield acceptable C IV profiles. For T_x = 0.5 MK, the model with no clumping (and, thus, the highest Ṁ) produces overly strong Ν V lines. For an increasing amount of clumping, Ν V decreases. Thus, our model with f_V,∞ = 0.5 yields the best fit for this line, making it our preferred solution. Models with other values of T_x are not sufficient as they fail to produce noticeable Ν V profiles.

In HD164353, the observed Si IV line is weaker than in our other targets, showing almost a "photospheric" appearance. In this case, the models with lower clumping yield an overpredic-tion of the Si IV lines. The observed unsaturated Si IV could only be reproduced with highly clumped models, implying an extremely low mass-loss rate (~10^-9), which would unfortunately spoil the reproduction of the C IV profile. This pattern of not being able to obtain an unsaturated Si IV together with a reasonable profile for C IV is similar to HD198478.

Overall X-ray limits. From the spectral fitting of our sample stars, we can obtain estimates of the X-ray luminosities intrinsically emitted by the winds of these stars. We present them in Table 6, where we also show the upper limits given by Berghoefer et al. (1997). If we allow for a very low shock temperature, which is not sufficient to produce N V, we derive values as low as log(L_x/L) ~ −13. Acknowledging that the N V line is barely noticeable, we consider this value to be a "hard" lower limit for the X-ray luminosities of these stars. However, if we take N V into account, we need to consider higher T_x and then we find values for log(L_x/L) between −7.0 and −8.3 for reasonable spectral fits. We thus see the latter value of ~ −8.3 as the more likely actual lower limit for the intrinsic L_x/L ratio of cooler BSGs.

Some of our models also yield values of ~ −7.0, thus matching the "canonical" value known for hotter massive stars (e.g., Sana et al. 2006). However, the results of our study do not favor these models as they require a high degree of clumping, which in general tends to weaken (and spoil) the carbon lines. Additionally, considering the whole sample, the best fits were obtained by models with f_V,∞ = 0.5 and 1.0 – namely, with mild and no clumping in the wind. Therefore, we consider lower values of clumping to be preferred for describing the winds of these cooler BSGs, which is in line with recent simulations and observational studies (Driessen et al. 2019; Rubio-Díez et al. 2022). Thus we provide additional empirical-based evidence that the winds of cool BSGs are smoother than those of O- and early BSGs.

In summary, we suggest that cooler BSGs must have X-ray luminosities that are below the "canonical" value of −7.0, namely, a value in the region between −7.3 and -8.3. These findings are in line with the observational constraints of Berghoefer et al. (1997) for our sample stars.

In appearent contrast to our results, Nazé (2009) measured X-ray (0.5 keV < hv < 10 keV) fluxes of OB stars with XMM-Newton and detected X-rays in seven cool BSGs. For some of their targets, they obtained log(L_x/L) = log F_x/F > −7, where F_x and F are their labels for the X-ray and bolometric unreddened flux respectively. However, none of their targets overlap with our sample. Moreover, their sample contains colliding wind binary candidates, as well as stars that have not been thoroughly investigated yet, and their bolometric fluxes were computed solely with a standard recipe based on bolometric corrections and intrinsic colors derived from the spectral type. Thus, further studies are required to check whether all cooler BSGs really have an intrinsic log(L_x/L) < −7.0. In any case, both our and their results underline that X-rays seem to be intrinsic to at least a significant fraction of the cool BSGs and need to be considered when modeling their atmospheres.

Fig. 5 Comparison between models (red dash-dotted and blue dashed lines) and observed (black full line) spectra of HD53138. The red model includes neither clumping nor X-rays in the wind and illustrates the current status in the literature regarding the UV spectral fitting of cool BSGs. The blue model includes clumping (fv,∞= 0.1) but no X-rays.

Fig. 6 Best fitting models for HD53138 in the UV region, focusing on the main P-Cygni profiles, namely N v λ21238-42, C II λ21335-36, S iiv λ1394-1403, C iv 21548-50 and Al III λ21855-63, depicted in each row. Each column of panels shows models with the same shock temperature. For each panel the red-dashed stands for f∞ = 1.0, orange/full for f∞, = 0.5, green-dash-dotted for f∞ = 0.2 and blue-dotted for f∞ = 0.1.

4.2 Effects of optically thick clumping

Previously, we demonstrated that the inclusion of (micro-)clumping, and especially X-rays substantially improved the fitting of the UV carbon profiles for our sample stars. However, all of our solutions for the different BSGs predict too strong Si IV λ1394-1403 and Al III λ1855-63. Moreover, the absorption part of the derived Hα profiles is in every case stronger than the observed spectrum, which may indicate a systematic problem in the ionization structure.

Analyzing BSG winds with the SEI method (Lamers et al. 1987), Prinja & Massa (2010) found that the lines in the Si IV λ1394-1403 doublet have optical depth ratios smaller than 2. They concluded that the profile must be affected by optically thick clumps in the wind. In BSGs, also Ha could behave like a resonance line, as discussed by Petrov et al. (2014). Hence, this profile could be affected by optically thick clumps as well.

Given this motivation, we tested the inclusion of macro-clumping as described in Oskinova et al. (2007). In this simplified initial approach, the macroclumping plays no role in impacting the density or ionization structures. Instead, as explained in Sect. 3.2, it changes the opacity in the computing of the observer's frame spectrum.

Keeping the properties we obtained with our CMFGEN model for HD53138, we calculated PoWR models with and without macroclumping for this star. For this task, we selected input options for PoWR that were as close as possible to the CMFGEN description, such as the radial clumping stratification or the considered X-ray luminosity. Despite the technical differences in the codes (cf. Appendix D), the PoWR models produce a spectrum similar to the CMFGEN models, reflecting the similar physical framework of the two codes. However, we did not aim for a precise benchmark between the two codes. Thus, we did not perform further iterations for a fine-tuned reproduction of the observation with the PoWR models employed for this task; rather, we focused on the effect on the main wind lines instead.

The synthetic spectra of our PoWR models with and without macroclumping are shown in Fig. 7. In order to obtain this outcome, we employed a typical separation between clumps L0 = 0.5 R_★. Notably, the most relevant impact of the macroclumping treatment is the decrease of the emission part of the P-Cygni profiles, except for Ha and Al III λ1855 to a small extent, since they are more affected in their absorption component.

The weakening of the absorption component in Hα and Al III λ1855 can be attributed to the reduction of the effective opacity. Due to the optically thick clumps, the number of available absorbing atoms¹⁰ is lower as a fraction of the material is shielded, and thus not contributing to the opacity of resonance lines (Oskinova et al. 2007). Likewise, the reduction in the effective opacity tends to reduce the emission components.

For Si IV λ1394-1403, despite improving the fit in the emission part, the major discrepancy (i.e., the non-observed saturation of the absorption) was not substantially affected and thus it remains unresolved by this approximate macroclumping treatment. The carbon lines, especially C II λ1335-36, only show minor effects in their emission parts. For these lines, as they were already closer to full saturation in the models without macroclumping, the impact of the reduction of the absorbers is limited.

In short, we obtain improved profiles for Ha and, to a smaller degree, in Al III and Si IV. However, in order to properly quantify the impact of macroclumping, we would need to integrate the effects of optically clumps into the density and opacity structure – for instance, an alteration in the opacity of strong profiles will affect the radiative transport within the wind, which could change the ionization structure. However, despite the approximate treatment of macroclumping, we provide evidence that optically thick clumping likely plays a role in shaping some of the wind-affected line profiles in later-type BSGs, which is in line with the findings of Prinja & Massa (2010) and Petrov et al. (2014).

Fig. 7 Impact of the inclusion of optically-thick clumps on the resolution of the formal integral on the stronger wind profiles. Both models have the same parameters, except for the presence of macroclumping in the formal integral. For the model that includes macroclumping (purple-solid line), we use L0 = 0.5 R★.

4.3 Mass-loss rates of cool BSGs

After deriving and discussing the clumping and X-ray properties of our sample, we now have sufficient information to compare our results with theoretical predictions and discuss potential consequences. Our wind properties are presented in Table 7. After testing different values for f_V,∞, we were able to confirm that the product $\dot{M}/\sqrt{f_{\text{V},\infty }}$ has to stay constant in order to yield good spectral fits. Thus, solutions with different f_V,∞ also have different M. In Table 7, we list only the preferred solutions according to the tailored analysis described on Sect. 4.1 with respect to the choice of volume filling factor. The f_V,∞ is indicated as a subscript to the value of Ṁ.

When considering the mass-loss rates re-scaled by the respective obtained radii, our analysis in general points towards lower mass-loss rates compared to previous work. For HD206165 and HD164353, the differences are most prominent, while the results for the two other targets are within the derived uncertainties of M (about 0.2 dex). Compared with previous works, the most significant distinction in our models are the presence of X-rays, which allowed us to successfully obtain the observed super-ionization and reproduce multiple UV wind lines – which were used as Ṁ diagnostics. This contrasts with the previous studies in the literature which needed to rely basically only Ha as a spectroscopic diagnostic. Additionally, we find that clumping was not a decisive factor for influencing the derived Ṁ because we found values that would imply (almost) smooth winds, close to the description employed in previous studies (e.g., Crowther et al. 2006; Searle et al. 2008). Moreover, Krtička et al. (2021), who produced cool BSGs models that include clumping but not X-rays, did not reproduce higher ions such as C IV (see their Fig. 9). In that line, we also noticed that in our set of models, those without X-rays, regardless of the clumping values were not able to produce superionization.

In Fig. 8, we compare the derived mass-loss rates values for each BSG with the mass-loss recipes of Vink et al. (2000), Krtička et al. (2021), and Björklund et al. (2022). Despite the preference for less clumped winds, we also include values of M for the other f_V,∞ tested as solutions with less clumped winds cannot be discarded (cf. Sect. 4.1). For comparison, we also included stars on the hot side of the bi-stability jump that were analyzed by different authors considering clumping and/or X-rays. To account for the different radii arising from different luminosities, we further compare the mass-loss rates scaled by their respective stellar radii following Puls et al. (2006). For the cool stars, all our obtained values are located between the predictions of Vink et al. (2000) and Björklund et al. (2022), independent of the clumping correction. In general, the Vink et al. (2000) predictions are about a factor of 2 or 3 too high, while the values from Björklund et al. (2022) are systematically lower than our derived results, usually by a factor of ~ 6. Instead, the predictions from Krtička et al. (2021) fit remarkably well with our results, except for the B2Ib star HD206165. For the O- and early B- supergiants from the literature, the Vink et al. (2000) discrepancy remains while the predictions of Björklund et al. and Krtička et al. start to coincide and reasonably reproduce the observations, except for the B0.5 supergiant HD167264. It is important to mention though that this star is part of a binary system (Sana et al. 2014; Mahy et al. 2022), and the parameters derived by Martins et al. (2015) were obtained without taking that into account.

The discrepancy between Krtička et al. (2021) and Björklund et al. (2022) for cool BSGs is striking. In contrast to the results from Krtička et al. (2021) and as well as the predictions by Vink et al. (1999), the models from Björklund et al. (2022) do not yield any increase of Ṁ when transitioning to lower T_eff values at a fixed luminosity. Empirical determinations of Ṁ (e.g., Crowther et al. 2006; Benaglia et al. 2007; Markova & Puls 2008) across this bi-stability jump region (BSJR) so far could not provide clear constraints for or against a jump in Ṁ, given also the spread in luminosities that complicate direct comparisons.

Our new results now prefer the Krtička et al. (2021) prescription, which includes an increase in Ṁ albeit in a quantitatively different way than the Vink et al. (1999) prescription commonly applied in stellar evolution models. However, the Krtička et al. (2021) recipe predicts the increase to start at lower T_eff than Vink et al. (1999, 2000, 2001), which seems to be contradicted by the B2 supergiant in our sample, which show a higher mass-loss than predicted by Krtička et al. (2021). It is evident that the biggest discrepancies with Krtička et al. (2021) are approaching around the “dividing” line of the bi-stability jump – T_eff ~ 22 kK – where supergiants of the spectral type B0.5 to B2 are located. We thus conclude that (at least in the case of supergiants) some form of mass-loss increase seems to happen when stars evolve towards cooler T_eff. The precise T_eff regime of this increase is hard to constrain with our limited sample but seems to be closer to the BSJR from Vink et al. (1999) as well as the more recent calculations by Vink (2018), albeit with an increase in Ṁ on the quantitative level found in the models from Krtička et al. (2021). A broader and more coherent sample of both hot and cool BSGs ought to be analyzed in order to draw a more conclusive picture.

Our insights on Ṁ also have potential consequences for the evolution of massive stars due to (i) the direct change of the resulting stellar mass and (ii) its impact on the loss of angular momentum (see, e.g., Keszthelyi et al. 2017). A significant increase of Ṁ in the BSJR, as discussed by Vink et al. (2010), could explain the steep drop trend in the observed rotation of BSGs below ~ 22 kK – the so-called “bi-stability breaking”. However, (Keszthelyi et al. 2017) argued that young O stars could alternatively lose rotational speed already before reaching the BSJR due to internal structural changes, thus questioning the need for a mass-loss induced breaking mechanism. Given that our results point in the direction of the Krtička et al. (2021) recipe, we would expect an increase in the loss of angular momentum, which is qualitatively also in line with the obtained low rotational velocities for our targets (~ 30 kms⁻¹). However, whether such an increase in Ṁ provides a quantitatively sufficient rotational breaking mechanism is yet to be tested by evolutionary models.

Fig. 8 Comparison between our derived mass-loss rates and the predictions according to different recipes, Vink et al. (2000), Krtička et al. (2021), and Björklund et al. (2022), for our sample stars shown in the upper panel. Early BSGs and OSGs, represented with cyan crosses, which were also analyzed using clumping and/or X-rays are included, namely: HD192639 (O7.5Iabf, Bouret et al. 2012), HD188209 (O9.5Iab, Marcolino et al. 2017), HD37128 (ϵ-Ori, B0Ia Puebla et al. 2016) and HD167264 (B0.5Ia/Iab, Martins et al. 2015). Bottom panel shows the same as the upper panel, but with all the stars scaled to the same radii following Puls et al. (2006) scaling relation of OB stars. HD167264 Ṁ values have fainter symbols in order to indicate that the parameters might not be very reliable as the star is a binary (Mahy et al. 2022).

Fig. 9 Position of our sample BSGs on the HR Diagram in comparison with evolutionary tracks (ZZAMS = 0.014, with rotation) from Ekström et al. (2012). The numbers at the beginning of the tracks represent the initial stellar mass and the black thin crosses indicate where the star leaves the main sequence ($X_{\text{H}}^{\text{core}}\to 0$). The cross with the red stroke marks the position of HD53138 (B3Ia).

5 Evolutionary context

5.1 Possible mass-discrepancy in BSGs

From the obtained stellar (photospheric) parameters, we can place our sample BSGs on the Hertzsprung-Russell diagram (HRD), presented in Fig. 9 and compared with evolutionary tracks from Ekström et al. (2012) at solar metallicities (assuming an initial rotation of υ/υ_crit = 0.4).

Both the position in the HRD and their surface CNO abundances place our sample BSGs as evolved stars, displaced from the main sequence. For such objects, we would expect that their surfaces display (significant) chemical alteration due to internal mixing processes (e.g. Ekström et al. 2012). For BSGs that evolve directly off the main sequence towards the red supergiant (RSG) stage¹¹, it is expected that the mixing occurs due to overshooting and rotation-driven large-scale circulation (e.g., meridional circulation Sweet 1950; Maeder 2009). Alternatively, stars can also become BSGs after having experienced a RSG stage. For these objects, one would expect higher surface chemical alteration due to internal convection at the RSG phase, and a reduced envelope due to the mass loss experienced in that extreme stage. This implies that post-RSG BSGs have a higher L/M ratio, and, consequently, different pulsational properties compared to pre-RSG BSGs as well (cf. Saio et al. 2013; Georgy et al. 2014, 2021).

In this paper, we perform an exploratory analysis of their status by comparing them with recent evolutionary tracks (Ekström et al. 2012). We determine the best-matching evolutionary models for our BSGs by performing a χ² fitting procedure taking into account different stellar parameters (see, e.g., Sect. 2 of Schneider et al. 2014, for an application using the Bonn tracks). The absoluteχ² values as such do not yield any judgment about our derived physical quantities, but provide us with a systematic approach to select the best-describing structure model from a given set of evolutionary tracks.

Initially, we considered L, T_eff, and log g for finding the best evolutionary track, as these are the main stellar properties and can directly be obtained via quantitative spectroscopy, alongside knowledge about the distance. The results are shown in Fig. 10. From the χ²-fitting, the estimated evolutionary masses are higher than the spectroscopic masses for HD 164353 and HD206165, albeit still within the mass uncertainty range for the latter. HD 198478 (B2.5Ia) reasonably agrees with the closest available track (M_ZAMS = 32 M_⊙) and HD53138 is best fitted by a 8 M_⊙ post-RSG model (M_ZAMS = 20 M_⊙). The mass discrepancy is most pronounced in the case of HD 164353, in which the evolutionary mass is about three times the spectroscopic determination. However, our derived χ² value for this object is also the largest one among our sample, thus indicating that HD 164353 might not be represented by any of the evolutionary tracks considered so far.

Given the properties of BSGs and their position in the HRD, it is not obvious whether individual objects are pre- or post-RSG stars. From the ;χ² minimization procedure, considering L, T_eff and log g, we find all but one object (HD53138) to be pre-RSGs.

One way proposed to discriminate their evolutionary status is via analyzing the surface CNO abundances (Hunter et al. 2009; Saio et al. 2013)¹². If we also take our derived surface CNO abundances at face value and account for them in the ;χ² fit (cf. the squares in Fig. 10), even HD53138 is considered to be a pre-RSG object. In that case, the spectroscopic mass (13 M_⊙) would be lower than the evolutionary model’s value (19 M_⊙).

Concerning the evolutionary status, when we compare our results with previous investigations, we find that for HD53138, Martin et al. (2018, citing Zorec et al. 2009) argued that there is no evidence that this BSG is a post-RSG object. For HD198478, Kraus et al. (2015) showed evidence that this object would be categorized as a post-RSG star due to its pulsational behavior (a Cyg variable) and the presence of a bow-shock feature in WISE images (4 μπι < λ < 24 μm), expected for stars with a previous phase of intense winds (e.g., Gordon & Humphreys 2019).

If we assume that our targets are pre-RSG objects (as indicated by the ;χ² fitting procedure accounting for the CNO abundance) our findings would suggest that three of our four cool BSGs might be affected by the so-called “mass discrepancy problem” (e.g., Herrero et al. 1992; Mokiem et al. 2007; Searle et al. 2008; Cantiello et al. 2009), which is commonly discussed for OB stars, where systematic differences are obtained between the masses inferred from spectroscopy and those derived from their luminosity via evolutionary models.

In Table 8, we quantify the possibility of a mass discrepancy for our sample by computing the probability of deriving lower spectroscopic masses with respect to the masses obtained from the ;χ² results. Thereby, the uncertainties in log L, T_eff and log g are taken into account (see Appendix A). We integrated the spectroscopic mass distribution for each star up to the value derived from the evolutionary tracks. For two cases, HD198478 and the putative post-RSG HD53138, the obtained evolutionary mass is close to the peak of the distribution. For the other two objects, however, the bulk of the distribution is lower than the mass from the evolutionary tracks (see the P_lower value in Table 8). Among these, HD206165 has a reliable Gaia DR3 distance and thus could be considered as a more robust example of a mass-discrepant BSG.

Among our sample, HD164353 exhibits the most extreme discrepancy with a spectroscopic mass of less than 3 M_⊙, while the corresponding evolutionary track describes a star with M_ZAMS = 9 M_⊙. From Table 8, the obtained value of P_lower indicates that it is extremely unlikely that we obtain this value by chance. From a single-star evolution perspective, a star born with ~ 9 M_⊙ would not be able to lose about 2/3 of its mass in an instrinsic way.

Errors in log g and L could be possible explanations. With higher log g, we could obtain higher masses. For example, with log g = 2.75, as found by Searle et al. (2008), we would obtain 7.4 M_⊙, which is much closer to the evolutionary track.

However, in our study, such a high value would yield significantly too wide Balmer line wings. This is illustrated in Fig. 11, where we compare different models with log g of 2.5 and 2.7, respectively. Alternatively, one could see this as a possible issue with the distance and thus the derived luminosity. If we use the Gaia eDR3 distance, as recently done by Weßmayer et al. (2022), we would obtain a higher luminosity and, thus, a higher mass (M = 13 M_⊙), which would yield properties more similar to other sample BSGs. However, as the Gaia parallax for this object is indicated as unreliable, we chose the HIPPARCOS parallax as discussed in Sect. 2.

Besides all the obtained statistics, there are important caveats to consider when trying to infer the evolutionary status of an object from comparing spectroscopic values to evolutionary tracks. First, we are limited to a discrete set of tracks with predetermined initial masses and initial rotation. Moreover, the masses and surface abundances of post-RSG stars are particularly uncertain as their properties depend on all the physics and assumptions included in the former evolutionary stages (e.g., Martins & Palacios 2013). For example, the lifetimes of massive stars in the BSG regime can depend significantly on the assumptions related to the internal mixing, such as the overshooting parameter (e.g., Schootemeijer et al. 2019; Higgins & Vink 2023). Moreover, the mass-loss rates in the RSG stage are still quite uncertain (see, e.g, the recent review by Decin 2021). Some recent empirical studies yield lower values for M (e.g., Beasor et al. 2021). However, if RSG mass-loss rates would be lower during the whole RSG lifetime, the resulting synthetic populations might be at odds with observed populations (e.g., Massey et al. 2023). Consequently, our CNO abundances are of limited diagnostic value as they could potentially also be reached in a post-RSG stage. Investigation into whether the “mass-discrepancy problem” is manifested in cool BSGs was previously undertaken (e.g., Searle et al. 2008; Trundle & Lennon 2005). Nonetheless, the comparisons between the spectroscopic results and evolutionary models were done “by eye.” From our analysis, given the small size of the sample, it is not possible to affirm as it does for other categories of OB stars. However, we present an approach that allow us to deliberate the existence of the problem in a mode systematic way. Thus, performing such an analysis on a larger sample would yield more concrete answers to this question.

Fig. 10 Individual comparison of each BSG in the HRD with each of the best-fitted points (through an χ2 minimization) at the evolutionary tracks of Ekström et al. (2012). The pink diamonds indicate the best model considering Teff, L, log g, while the white squares indicate the same but with the CNO abundances taken into account. For most stars, except HD53138, the result is exactly the same.

Fig. 11 Comparison between models with different log g for HD164353. The left panels show the Balmer lines for the model with log g = 2.5 (brown lines), while the right panels depict the model with log g = 2.7 (green line).

5.2 Investigating multiplicity

In general, BSGs are difficult to understand in the standard picture of stellar evolution because they are observed in a region of the HRD where massive stars spend very little time (thermal-timescale evolution through the Hertzsprung gap) assuming “typical” mixing settings. For example, the best-fitting singlestar models found here spend ≲ 10⁴ years in that region. Stellar mergers in binary stars can lead to the formation of stars that spend their entire core-helium burning lifetimes as BSGs (e.g., Hellings 1983; Braun & Langer 1995; Podsiadlowski & Joss 1989; Podsiadlowski et al. 1990; Claeys et al. 2011; Vanbeveren et al. 2013; Justham et al. 2014).

To understand whether merging can explain the observed stellar properties of the BSGs studied here, we mimic the merger process by rapid accretion of mass onto stars that just finished core hydrogen burning. In that procedure, we do not account for alteration in the surface chemical composition in the merging process. We employed MESA revision 10398 (Paxton et al. 2011, 2013, 2015, 2018, 2019) and use the same model setup as in Schneider et al. (2021). The merged stars can simultaneously reproduce the observed luminosities, effective temperatures and surface gravities while burning helium in their cores (see e.g., Fig. 12).

As with the single-star models, we tend to find masses higher than the spectroscopic determinations, especially for HD164353, namely, the L/M-ratio of the relaxed merger model is not much different from the L/M ratio of a single-star evolutionary model with a higher initial mass. Therefore, it does not resolve the mass discrepancy problem. However, the merger models spend ~ 10⁵ yr in the location of the HRD where the BSGs are observed, making a merger solution for them much more likely compared to a single-star scenario. This idea can be further investigated by, for example, searching for possible, close binary companions, and having more detailed merger simulations that trace the alterations in surface chemical abundances for comparison with the measured composition. This idea can be further investigated by, for example, searching for possible, close binary companions, and having more detailed merger simulations that trace the alterations in surface chemical abundances, for comparison with the measured composition. Although it does not provide a clear solution to the mass discrepancies, this result provides evidence that many BSGs could be products of stellar mergers.

Alternatively, the discrepant masses could be the result of stellar stripping. There are indications that some partially stripped stars can appear as blue supergiants. As discussed for example in Klencki et al. (2022), such objects would display a lower spectroscopic mass (relative to the “evolutionary mass”), CO depletion plus N enrichment, and slow rotation. All of these characteristics are present in most of our sample stars, especially in HD164353. However, their results were obtained at LMC metallicity and the scenario did not occur at solar metallicity. Moreover, much lower log g values than our determinations are expected. On the other hand, Irrgang et al. (2022) analyzed the B2.5 subgiant γ Col and identified it as a stripped, pulsating core of a previously more massive star. They found a strong nitrogen enrichment of ~ 10 times the solar value. While our sample stars do show N enrichment, our factors of up to 4 are significantly below this value.

For our targets, there is no real evidence for existing companions so far. Of course, the large luminosities of BSG-like stars can easily outshine darker companions, but for systems that interacted in the past one would expect to see radial velocity shifts when comparing spectra from different epochs, as for example in the high-mass X-ray binary (HMXB) QV Nor (B0I+NS), which has a mass ratio of ~15 and a clearly detected orbital motion of the donor star (Reynolds et al. 1992). In the case of such compact companions, the detection of a high X-ray luminosity, caused by the accretion of (wind) material on the compact object, would be a clear indication. However, as we discuss in Sect. 4, our sample does not contain objects that are overluminous in the X-ray regime. Moreover, wind-fed HMXBs typically have B0-B1 type donors (Martínez-Núñez et al. 2017). From the current binary evolutionary scenarios and their predictions, it thus seems unlikely that our sample stars are part of a binary system with a “recent” interaction.

In summary, the derived properties of our BSGs, except for the luminosity and mass of HD164353, are in the range of properties found in other BSG studies (cf. Searle et al. 2008; Fraser et al. 2010; Haucke et al. 2018). Our obtained mass discrepancy is relatively common, which leaves the puzzling question of whether the observed population of BSGs is dominated by immediate post-MS stars or consists of a significant number of more evolved objects such as post-RSG stars, stellar mergers, or partially stripped post-interaction binary products. A more extensive and coherent sample study, such as provided by ULLYSES (Roman-Duval et al. 2020) and its accompanying X-shooter program (XShootU, Vink et al. 2023), as well as (if possible) more precise mass measurements will be key to answering this question.

Fig. 12 Comparison between single evolutionary tracks Ekström et al. (2012) and a MESA model with a merger (purple broader line) for HD206165. The squares with numbers are the ages in units of 105 yr for the Ekström et al. (2012) tracks while the white squares show the same time difference for the merged model. The shaded region in pink shows the core-He burning phase. The arrow indicates the point of the 15 M⊙ evolutionary track where the merging (mass addition of 7 M⊙) happened.

6 Conclusions

In this study, we performed an in-depth spectral analysis of four Galactic cool BSGs with CMFGEN and PoWR, covering the UV and optical regime. For the first time in this regime, we simultaneously reproduced the observed optical spectrum together with major wind-affected UV profiles. We obtained the photospheric and wind properties of our sample stars as well as constraints on their clumping and X-ray properties. From the immediate spectral analysis, we drew the following conclusions:

• By including clumping, and especially X-rays, in the atmosphere models of cool BSGs, we can fit key UV wind lines and generally reproduce the observed superionization. This represents a major improvement in spectral modeling of later type BSGs.

• We find lower relative X-ray luminosities than the “canonical” relation of log(L_x/L) = −7.0. Assuming a shock temperature compatible with the typical wind velocities of our sample stars (T_x = 0.5), we obtain log(L_x/L) in the range of −8.3 to −7.3.

  In some spectra, it is unclear whether N V is actually present. Disregarding N V, we find even lower limits of log(L_x/L) ~ 13.3 and shock temperatures of only 0.1 MK. This happens independently of the adopted clumping values. Revisiting cooler BSGs with new X-ray observations will be necessary to put more robust constraints on the wind-intrinsic X-ray emission of these stars.

• Our results produce further evidence that the winds of BSGs are smoother than the winds of their hot counterparts as recent hydrodynamical simulations predict (Driessen et al. 2019). In the microclumping formalism, our targets are described by volume filling factors of 0.5 to 1.0. Additionally, we found evidence for the presence of optically thick clumping, which improves the Hα fit when accounted for in the spectral modeling.

  From our best-fit models, we investigated the consequences of the derived stellar and wind parameters on their evolutionary status:

• We find stellar properties compatible with evolved stars, namely altered CNO abundances and a clear displacement from the main sequence in the HRD. Our results are in line with previous literature providing properties of cool BSGs.

• For two of our four targets, our derived spectroscopic masses are lower than those expected from immediate post-MS stars assuming single-star evolution (including rotation). Our systematic analysis indicates tentative evidence that the so-called “mass discrepancy problem” that is known for many OB stars might also apply to at least some BSGs. A more in-depth look into this problem with a larger sample is necessary to get a more general overview of the situation in the BSG regime.

• Our results indicate that at least a sub-sample of cool BSGs might not stem from a single-star evolution. Given the intrinsically high binary fraction among massive stars, it is likely that some BSGs are products of binary interaction (e.g., stripping or merger). We find possible merger scenarios to explain some of our objects, but the corresponding evolution models do not resolve the mass discrepancy problem.

Besides the individual status of each object, we can also draw more general conclusions about massive stars and their evolution from the derived wind parameters:

• Our derived mass-loss rates are systematically lower than values from previous studies, albeit usually within the uncertainty range. We attribute this to the effect of the included X-rays on the UV P Cygni profiles used for mass-loss diagnostics.

• When comparing our findings to different mass-loss recipes, we find that our cool BSGs are significantly lower than predicted by Vink et al. (2001), but generaly align remarkably well with the predictions by Krtička et al. (2021). Including early-B and O supergiants for comparison, we observe the same agreement. For these earlier-type stars, there is a good agreement between the predictions by Krtička et al. (2021) and Björklund et al. (2022). This, however, does not happen for the late BSGs, where the Björklund et al. rates seem to be too low.

• In contrast to the Björklund et al. (2022) predictions, our findings provide strong evidence that an increase in the mass-loss rates does occur when a supergiant crosses the bi-stability jump region during its evolution. The increase, however, seems to be less pronounced than predicted by Vink et al. (1999, 2000, 2001), but rather on the level predicted by Krtička et al. (2021). Interestingly, the characteristic temperature of the switch seems to align more with the predictions of Vink et al. (1999) than the Krtička et al. (2021) recipe. A detailed wind-driving investigation on the bi-stability region will be necessary to draw more robust conclusions about where and under which circumstances the mass-loss rate increases in this regime.

Overall, our study reveals the importance of accounting for X-rays in the study of cooler B supergiants, despite a so far quite limited observational constraints. Moreover, our findings underline that the puzzling status of B supergiants is far from being solved and a larger sample study with multi-wavelength observations is required. Additional information from variability studies and asteroseismology as well as theoretical and numerical studies – such as dynamically-consistent atmospheres, 3D structure modelling, or interacting binary evolution – will be required to fully determine the nature of these crucial objects at the crossroads of massive star evolution.

Acknowledgements

The author would like to thank the referee who provided insightful and constructive comments that helped to improve this manuscript. M.B.P., A.A.C.S., and V.R. 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). M.P.B., A.A.C.S., and V.R. further acknowledge support from the Federal Ministry of Education and Research (BMBF) and the Baden-Württemberg Ministry of Science as part of the Excellence Strategy of the German Federal and State Governments. FRNS received funding from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme (Grant agreement No. 945806) and is supported by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) under Germany’s Excellence Strategy EXC 2181/1-390900948 (the Heidelberg STRUCTURES Excellence Cluster). F.N. acknowledges grant PID2019-105552RB-C41 funded by the Spanish MCIN/AEI/ 10.13039/501100011033.

Appendix A Quantifying uncertainties in mass and radius

In this study, we quantify the errors in M and R by considering that the errors in luminosity (L), effective temperature (T_eff) and surface gravity (log g) follow a Gaussian distribution around their nominal values for each star. Based on this, we then compute the mass distribution for each BSG. The mean value of each distribution corresponds to the nominal value of M (or R) and the standard deviation is adopted as the corresponding uncertainty. The results of this procedure are illustrated for M in Fig. A.1, where we show the mass distribution for each star together with the predicted evolutionary mass from the χ²-flt (see Sect. 5).

Fig. A.1 Histograms depicting the probability distribution of the masses of each BSG in the sample. The blue-filled and red-dashed lines indicate, respectively, the spectroscopic and evolutionary masses obtained. The grey dotted line for HD53138 indicates the evolutionary mass obtained accounting for the CNO abundance.

Appendix B UV wind lines of the models with different clumping and X-rays combinations.

In this appendix, we display the models for each clumping and X-rays settings in more detail for HD53138 and for the other sample stars (see Figs B.1 and B.2 respectively).

Fig. B.1 Best-fitting models for HD53138 in the UV region, focusing on the main P-Cygni profiles, namely N V λ1238-42, C II λ1335-36, Si IV λ1394-1403, C IV λ1548-50 and AI III λ1855-63, depicted in each column. Each group of four rows has the same shock temperature, in which each row has a model with an assigned clumping and X-ray luminosity.

Fig. B.2 Best-fitting models for the main UV lines of HD206165, HD198478 and HD164353. Each row of plots shows the lines of the same model with a given shock temperature Tx. The different colors/line styles represent the different adopted values of clumping, namely: red-dashed stands for f∞ = 1.0, orange/full for f∞ = 0.5, green-dash-dotted for f∞ = 0.2 and blue-dotted for f∞ = 0.1. The associated log(Lx/L) and Ṁ for each clumping and shock temperature value are displayed in Tables 6 and 7.

Appendix C Additional information and comments about the sample BSGs

Appendix C.1 HD206165 (B2Ib)

We note that 9 Cep is a star that has not been studied in great detail previously. Given its spectral type, it was assigned a temperature of ~ 20 kK, very close to the bi-stability jump, by McErlean et al. (1999). Its spectrum was modeled by Searle et al. (2008) and Markova & Puls (2008), using CMFGEN and FASTWIND respectively. Searle et al. found a temperature and surface gravity similar to our values while Markova & Puls found a similar log g but higher temperature (19.6kK).

The Combined General Catalogue of Variable Stars (CGVS Samus et al. 2004) shows that the star has a magnitude variation of ~ 0.1 and classified it as an α-Cygni variable (but with no evidence of periodic pulsation). However, we did not notice any meaningful alteration in the spectral lines between different available data (2007 and 2011).

For the wind terminal velocities, both authors following Howarth et al. (1997) used υ_∞ = 640 km s⁻¹. However, we could only obtain a good fit for the UV P-Cygni lines – except for C II, which is even broader – when applying υ_∞ = 900 km s⁻¹.

Appendix C.2 HD198478 (B2.5Ia)

HD198478 is a BSG that has been analyzed in greater detail in the past decades by different authors using multiple techniques, namely: Crowther et al. (2006) and Searle et al. (2008) who used CMFGEN to model its optical+UV spectrum, Kraus et al. (2015) who employed FASTWIND to model tens of optical spectra collected in a five-year campaign, and Gordon et al. (2019) who measured some of its properties using interferometry from the CHARA array. This makes HD198478 one of the most well-studied cooler BSGs so far.

Despite the multiple analyzes performed for this, there is quite a broad discrepancy between these different works regarding the obtained properties for HD198478. For instance, Crowther et al. (2006, similarly to our work) found T_eff = 16.5 kK and log g = 2.15 using CMFGEN, while Kraus et al. (2015) used FASTWIND to find T_eff ≈ 19.0 kK and log g ≈ 2.4. Markova & Puls (2008), also using FASTWIND, found an intermediate value of 17.0kK.

Regarding the wind velocity, the situation is equally discrepant with terminal velocities found in the range between 200 and 250 km s⁻¹ (Markova & Puls 2008; Kraus et al. 2015) and 550km s⁻¹ (Searle et al. 2008). In this work, we derived υ_∞ = 570km s⁻¹.

The main physical parameters obtained for HD198478 are notably close to those of HD53138 (B3Ia) – which is not totally surprising given that they have a very similar spectral type and their optical spectra are almost identical.

Appendix C.3 HD53138 (B3Ia)

24-CMa or o²-CMa is a relatively well-studied BSG, hence, it is our choice for a more in-depth analysis. Its atmosphere was modeled by several authors using CMFGEN and FASTWIND (e.g., Crowther et al. 2006; Lefever et al. 2007; Searle et al. 2008; Haucke et al. 2018).

As with many other BSGs, this star displays spectral variability in Hα (as discussed in Sect. 2). Additionally, it also possesses weak photometric variability detected by TESS and HIPPARCOS (Bowman et al. 2022; Lefèvre et al. 2009). From that variability and its relatively elevated macroturbulent velocity, Haucke et al. (2018) also classified this star as an α-Cygni variable.

HD53138 was further scrutinized by the LIFE project (Martin et al. 2018), which seeks for magnetic fields in OB stars. However, no evidence of magnetic activity was detected.

Appendix C.4 HD164353 (B5lb/ll)

67 Oph is a BSG on the edge of the bright giant classification – which is in line with its larger Balmer lines and higher inferred surface gravity. However, HD164353 has not been extensively studied in the literature. Zorec et al. (2009) analyzed this star looking at broader features in the spectrum and reported a T_eff of 15.4 kK.

Currently, the only works which analyzed this star using a modern atmosphere code are Searle et al. (2008) and Weßmayer et al. (2022). Searle et al. derived a temperature (T_eff = 15.5 kK), luminosity (log(L/L_⊙) = 4.3) and radius (16 R_⊙) that are similar to the values we obtained (see Table 4). However, their log g is 0.3 dex higher than ours, which results in a stellar mass twice as large as ours. Weßmayer et al. found nearly identical parameters to our study. However, they considered the Gaia EDR3 distance, which is twice the value from HIPPARCOS and we disregarded it due to its high RUWE value. Consequently, Weß-mayer et al. obtained a higher value for the luminosity and mass of HD164353.

Vaiana et al. (1981) reported the detection of X-rays (>0.2 keV, using the Einstein X-ray telescope) with a log(L_X/L) = −6.9. Employing ROSAT, Berghoefer et al. (1997) only obtained the same ratio, but only as an upper limit for the X-ray flux.

Appendix D Further details on the X-ray parameters

Appendix D.1 Choice of X-ray onset velocity

In CMFGEN, wind-intrinsic X-rays are implemented via an additional X-ray emissivity η_x calculated via

$${\eta }_{\text{x}}\left(r,v\right)\propto f_{\text{x}}^2{e}^{-{\upsilon }_{\text{x}}/\upsilon \left(r\right)}\frac{e^{-hv/k_{\text{B}}{T}_x}}{\sqrt{T_x}}{\displaystyle \sum _j^{\text{all ions}}{z}_j^2{g}_{\text{ff}}\left(v,r\right)n_{\text{e}}\left(r\right)n_j\left(r\right)},$$(D.1)

where υ_x is the characteristic velocity and f_x is a filling factor, which translates the amount of shocks producing X-ray in the wind. T_x is a characteristic temperature reached in the shocks, which controls the X-ray distribution in the frequency space. The quantities n_e and n_j are the density of electrons and ions respectively, g_ff is the free-free Gaunt factor and k_B and z are, respectively, the Boltzmann constant and the charge of the ion j.

Focusing on the space dependency (i.e., the exp[−υ_x/υ(r)] term), υ_x determines the spatial distribution of the X-ray emissivity. The characteristic velocity describes the point where X-rays become meaningful in the wind, where the profile function would reach ~0.37, namely, ~37% of the maximum reachable value.

In Fig. D.1, we show an example of the velocity-dependence of η_x. A low value of υ_x describes a wind where the X-ray emissivity is essentially ubiquitous in the wind. In contrast, a high υ_x severely reduces the emissivity in the innermost wind regions (i.e., at low velocities) and describes a wind where X-rays are progressively getting more relevant towards the outermost regions.

Fig. D.1 Velocity dependence of the X-ray emissivity in CMFGEN for different onset velocities: υx/υ∞ = 5 % in red and υx/υ∞ = 70 % in blue. The black dotted horizontal line marks the point where −υ∞/υ(r) = 1. which can be physically interpreted as where the shocks in the wind would roughly become relevant to produce X-rays.

Fig. D.2 Comparison between cool BSG models with low and high υx/υ∞ respectively (left and right columns of panels respectively).

Appendix D.2 CMFGEN vs. PoWR descriptions of X-rays

As mentioned in Sect. 3, PoWR includes X-rays as thermal Bremsstrahlung emission similarly to CMFGEN. The resulting expression for η_x is expressed as:

$${\eta }_{\text{x}}\left(r,v\right)\propto f_{\text{x}}^{2}H\left(r;r_0\right)\frac{e^{-hv/k_{\text{B}}{T}_x}}{\sqrt{T_x}}{\displaystyle \sum _j^{\text{all ions}}{z}_j^2{g}_{\text{ff}}\left(v,r\right)n_{\text{e}}\left(r\right)n_j\left(r\right)},$$(D.2)

where the repeated symbols from Eq. (D.1) have the same meaning. A major difference is the spatial distribution of the emissivity, which is modeled by a Heaviside step function¹³ centered in r₀ H(r; r₀) in PoWR. This describes a situation where the X-rays start sharply in the wind beyond a certain radius (or velocity), which differs from the gradual build-up in CMFGEN.

Given the different implementations, transitions affected by X-rays could produce slightly different UV line profile shapes despite having similar L_x, T_x and “onset” descriptions (i.e., υ_x or corresponding radius). In this study, our reported parameters are all taken from CMFGEN and solely employ PoWR to study potential impacts of optically thick clumping, so the difference in the descriptions is not relevant. Nonetheless, it implies that X-ray luminosities derived only from UV line profiles have a non-negligible dependence on the employed atmosphere code. Detailed observations and multi-dimensional hydrodynamic simulations could prove helpful in gaining quantitative insights on the spatial origin of wind-intrinsic X-rays in order to improve the parametrization in atmosphere models.

Appendix E SEDs, Optical and UV spectra

In this section, we present combined figures showing the spectral energy distribution (SED) and the normalized spectra. In Fig. E.1, the different panels show the different targets studied in this work. The model SEDs are compared to observed flux-calibrated UV spectra (IUE and if available also FUSE) as well as the flux values obtained from the available magnitudes spanning from the UV to the near IR (see Table E.1 and Sect. 2 for details). Then, Figs. E.2 to E.5 show the normalized observations in the UV and optical range compared to the best-fitting model spectra.

Fig. E.1 Spectral Energy Distributions of the sample cool BSGs.

Fig. E.2 UV and optical spectra of HD206165.

Fig. E.3 UV and optical spectra of HD198478.

Fig. E.4 UV and optical spectra of HD53138.

Fig. E.5 UV and optical spectra of HD164353.

References

All Tables

All Figures

Fig. 1 UV spectra of our sample stars – CIV 1550 and Si IV 1400 region. Note: variability is absent or mild in the available IUE data as the profiles do not present drastic changes (e.g. appearance/destruction of P-Cygni profile or variations in the width). The black profiles are the computed average IUE spectra which were used for the model comparison. For HD164353, the green lines represent small aperture IUE data (re-scaled to match the flux of the large aperture acquisition). The black thick lines are the averaged spectra which were used for model comparison (see Sect. 3).

In the text

	Fig. 2 Time sequences of Ha profiles for our sample BSGs. The high variability in HD53138 and HD198478 makes the profile not very suitable for obtaining reliable wind properties without auxiliary diagnostics.
In the text

	Fig. 3 Comparison between synthetic and observational spectrum of HD53138 in the Ha region. The color scheme is the same as in Fig. 2 while the black filled line represents the preferred model while the and dotted lines represent models with an increase/decrease in M by a factor of 30%.
In the text

	Fig. 4 Optical spectrum of HD53138 (black solid line) compared to our best-fit model (red-solid line). Here we depict the main diagnostic lines used to obtain the stellar photospheric properties.
In the text

	Fig. 5 Comparison between models (red dash-dotted and blue dashed lines) and observed (black full line) spectra of HD53138. The red model includes neither clumping nor X-rays in the wind and illustrates the current status in the literature regarding the UV spectral fitting of cool BSGs. The blue model includes clumping (fv,∞= 0.1) but no X-rays.
In the text

Fig. 6 Best fitting models for HD53138 in the UV region, focusing on the main P-Cygni profiles, namely N v λ21238-42, C II λ21335-36, S iiv λ1394-1403, C iv 21548-50 and Al III λ21855-63, depicted in each row. Each column of panels shows models with the same shock temperature. For each panel the red-dashed stands for f∞ = 1.0, orange/full for f∞, = 0.5, green-dash-dotted for f∞ = 0.2 and blue-dotted for f∞ = 0.1.

In the text

	Fig. 7 Impact of the inclusion of optically-thick clumps on the resolution of the formal integral on the stronger wind profiles. Both models have the same parameters, except for the presence of macroclumping in the formal integral. For the model that includes macroclumping (purple-solid line), we use L0 = 0.5 R★.
In the text

Fig. 8 Comparison between our derived mass-loss rates and the predictions according to different recipes, Vink et al. (2000), Krtička et al. (2021), and Björklund et al. (2022), for our sample stars shown in the upper panel. Early BSGs and OSGs, represented with cyan crosses, which were also analyzed using clumping and/or X-rays are included, namely: HD192639 (O7.5Iabf, Bouret et al. 2012), HD188209 (O9.5Iab, Marcolino et al. 2017), HD37128 (ϵ-Ori, B0Ia Puebla et al. 2016) and HD167264 (B0.5Ia/Iab, Martins et al. 2015). Bottom panel shows the same as the upper panel, but with all the stars scaled to the same radii following Puls et al. (2006) scaling relation of OB stars. HD167264 Ṁ values have fainter symbols in order to indicate that the parameters might not be very reliable as the star is a binary (Mahy et al. 2022).

In the text

Fig. 9 Position of our sample BSGs on the HR Diagram in comparison with evolutionary tracks (ZZAMS = 0.014, with rotation) from Ekström et al. (2012). The numbers at the beginning of the tracks represent the initial stellar mass and the black thin crosses indicate where the star leaves the main sequence ($X_{\text{H}}^{\text{core}}\to 0$). The cross with the red stroke marks the position of HD53138 (B3Ia).

In the text

	Fig. 10 Individual comparison of each BSG in the HRD with each of the best-fitted points (through an χ2 minimization) at the evolutionary tracks of Ekström et al. (2012). The pink diamonds indicate the best model considering Teff, L, log g, while the white squares indicate the same but with the CNO abundances taken into account. For most stars, except HD53138, the result is exactly the same.
In the text

	Fig. 11 Comparison between models with different log g for HD164353. The left panels show the Balmer lines for the model with log g = 2.5 (brown lines), while the right panels depict the model with log g = 2.7 (green line).
In the text

Fig. 12 Comparison between single evolutionary tracks Ekström et al. (2012) and a MESA model with a merger (purple broader line) for HD206165. The squares with numbers are the ages in units of 105 yr for the Ekström et al. (2012) tracks while the white squares show the same time difference for the merged model. The shaded region in pink shows the core-He burning phase. The arrow indicates the point of the 15 M⊙ evolutionary track where the merging (mass addition of 7 M⊙) happened.

In the text

	Fig. A.1 Histograms depicting the probability distribution of the masses of each BSG in the sample. The blue-filled and red-dashed lines indicate, respectively, the spectroscopic and evolutionary masses obtained. The grey dotted line for HD53138 indicates the evolutionary mass obtained accounting for the CNO abundance.
In the text

	Fig. B.1 Best-fitting models for HD53138 in the UV region, focusing on the main P-Cygni profiles, namely N V λ1238-42, C II λ1335-36, Si IV λ1394-1403, C IV λ1548-50 and AI III λ1855-63, depicted in each column. Each group of four rows has the same shock temperature, in which each row has a model with an assigned clumping and X-ray luminosity.
In the text

Fig. B.2 Best-fitting models for the main UV lines of HD206165, HD198478 and HD164353. Each row of plots shows the lines of the same model with a given shock temperature Tx. The different colors/line styles represent the different adopted values of clumping, namely: red-dashed stands for f∞ = 1.0, orange/full for f∞ = 0.5, green-dash-dotted for f∞ = 0.2 and blue-dotted for f∞ = 0.1. The associated log(Lx/L) and Ṁ for each clumping and shock temperature value are displayed in Tables 6 and 7.

In the text

	Fig. D.1 Velocity dependence of the X-ray emissivity in CMFGEN for different onset velocities: υx/υ∞ = 5 % in red and υx/υ∞ = 70 % in blue. The black dotted horizontal line marks the point where −υ∞/υ(r) = 1. which can be physically interpreted as where the shocks in the wind would roughly become relevant to produce X-rays.
In the text

	Fig. D.2 Comparison between cool BSG models with low and high υx/υ∞ respectively (left and right columns of panels respectively).
In the text

	Fig. E.1 Spectral Energy Distributions of the sample cool BSGs.
In the text

	Fig. E.2 UV and optical spectra of HD206165.
In the text

	Fig. E.3 UV and optical spectra of HD198478.
In the text

	Fig. E.4 UV and optical spectra of HD53138.
In the text

	Fig. E.5 UV and optical spectra of HD164353.
In the text