© The Authors 2025

1 Introduction

B hypergiants (BHGs) are a rare type of evolved massive stars. To date, slightly fewer than 40 objects are known in the whole Local Group (see Table A.1). Farther away, B stars, with luminosities comparable to BHGs, are detected at high redshifts (z ≳ 6.2) via gravitational lensing (Welch et al. 2022). However, only five BHGs have so far been properly analyzed with quantitative spectroscopy, and the physical mechanisms responsible for their stellar wind driving are not yet understood.

In the Hertzsprung-Russell diagram, BHGs are located near B supergiants (BSGs). However, their usually high luminosities place them closer to luminous blue variables (LBVs). Spectroscopically, they differ from the “typical” BSGs by the presence of at least Hβ in emission (pure emission or a P-Cygni profile; see, e.g., Lennon et al. 1992). However, it is also common to observe many more H and He lines appearing as P-Cygni profiles (e.g., Clark et al. 2012; Walborn et al. 2015; Mahy et al. 2016).

In recent years, several pioneering efforts have been made to model radiative-driven winds with hydrodynamically consistent atmosphere models (i.e., no longer relying on the use of ad hoc velocity fields like the β law). These studies significantly advanced our understanding of Wolf-Rayet (WR) stars (Sander & Vink 2020; Poniatowski et al. 2021; Sander et al. 2023), O stars (Krtička & Kubát 2017; Sander et al. 2017; Sundqvist et al. 2019; Björklund et al. 2021), and early-B stars (Sander et al. 2018; Krtička et al. 2021; Björklund et al. 2023). A few wind investigations were also performed for the later B-star regime (e.g., Petrov et al. 2016; Vink 2018; Venero et al. 2024), albeit with more approximate wind dynamics. Yet, for BHGs and LBVs, which are closer to the Eddington limit and may play crucial roles in high-mass stellar evolution, there is a severe lack of dedicated investigations. In this study we present the first hydrodynamically consistent atmosphere modeling of an early-type BHG¹.

Using PoWR^HD (see Sander et al. 2017), we can predict the mass-loss rate and velocity field directly from a given set of stellar parameters. To validate and benchmark our dynamically consistent model, we compared the obtained synthetic spectra with observed spectra of ζ¹ Sco (B1.5Ia+) from the UV to the IR. ζ¹ Sco is currently the brightest BHG in the sky (V = 4.78, Zacharias et al. 2004) and was chosen as a prototype as it is the Galactic BHG with the broadest spectral coverage. We further compared our results with previous work using “traditional” (i.e., non-hydrodynamically consistent) atmosphere models (e.g., Crowther et al. 2006; Clark et al. 2012).

In Sect. 2, we briefly describe PoWR and the assumptions adopted in our modeling. In Sect. 3 we present and discuss the observational data and compare them to our model spectra in different wavelength regions. In Sect. 5, we discuss the obtained stellar properties, and in Sect. 6 we discuss the wind properties, focusing on the dynamics, mass-loss rates, and wind clumping. In Sect. 7 we discuss how BHGs compare to Of/WNh stars and how their winds are presumably transitioning to the “optically thick wind regime”. We present our conclusions in Sect. 8.

2 Hydrodynamically consistent atmosphere model

For our investigation, we made use of the state-of-the-art PoWR stellar atmosphere code PoWR (e.g., Gräfener et al. 2002; Sander et al. 2015). More specifically, we used PoWR^HD branch (Sander et al. 2017, 2018, 2023), which solves the momentum equation consistently with the radiation field, temperature stratification, and the ionization/excitation rate equations at each layer of the atmosphere/wind. During the model iteration the velocity/density field is updated so that the acceleration terms (from radiative and pressure forces) are equal to the deceleration terms (due to gravity and inertia) in each atmosphere layer (see Sander et al. 2017, 2018, 2023 for a detailed explanation of the implementation).

One key aspect of consistent stellar atmosphere codes is that the wind structure and terminal velocity v_æ is neither an ad hoc parametrization (as the traditionally used β law) nor an input property. Rather, it is calculated directly from the fundamental properties and inner-boundary conditions of the model (e.g., radius, luminosity, and mass). Likewise, the mass-loss rates, M (or alternatively the mass, M) is also computed consistently. However, given the 1D and stationary nature of stellar atmosphere codes, physical effects such as clumping and X-ray emission still need to be described parametrically.

Clumping description. Both theory and observation reveal that the winds of hot stars are inhomogeneous. The stellar atmospheres are expected to have regions of over-densities (clumping) and rarefaction in their outflows. Clumping affects the formation of lines in various direct and indirect ways. For instance, (i) the strength of recombination lines, such as Balmer lines, dependent on ρ², and (ii) scattering resonance lines, such as UV P-Cygni profiles, are affected by changes in the ionization structure caused by the alterations in the density structure. We applied the exponential parametrization of the clumping structure introduced by Hillier et al. (2003),

$$f\left(r\right)=f_{\infty }+(1-f_{\infty })\exp (-\frac{v\left(r\right)}{v_{\mathrm{cl}}}),$$(1)

for the volume filling factor f, which is the inverse of the clumping factor (f_cl = 1/f) if clumps are assumed to be optically thin and the inter-clump medium is void. The same has previously be used for BHGs in Clark et al. (2012). In this description, the clumping profile depends only on two parameters f „ and -_cl, where the former controls the degree of inhomogeneity at an “infinite” distance, and the latter is the onset characteristic velocity, which controls how abruptly clumping will start growing in the wind, essentially working as a scale height.

X-rays. As discussed in Bernini-Peron et al. (2023, 2024), the inclusion of X-rays in cool BSG models (T_eff ~ 21 kK) is necessary to produce the so-called superionization (e.g., Cassinelli et al. 1981; Odegard & Cassinelli 1982; Krtička & Kubát 2016), characterized by the presence of N V, C IV, and Si IV lines in the UV with well-developed P-Cygni profiles, which are not expected to appear at such low effective temperatures. Conversely, the extra source of ionization is also necessary to quench the profiles of C II and Al III, which would be saturated otherwise. The neglect of X-rays would likely yield to an underestimation of the mass-loss rate, as the saturation of C II and Al III would then have to be done by assuming thinner winds. For BHGs, this is not an option, as they also show well-developed NV, CIV, and Si IV profiles, indicating that superionization needs to happen in their atmospheres despite their comparably slow winds (v_æ < 500 km s^–1). Hence, we included X-rays in our models.

PoWR implements X-rays as free-free emission (bremsstrahlung), following the formalism of Baum et al. (1992). While not including specific line contributions, this formalism is sufficient to provide the necessary radiation field at high frequencies to produce superionization, via either the Auger process or direct ionization due to a hot thermal component. Isolated BHGs have not yet been detected in X-rays. Thus, superionization provides the only means to constrain the X-rays in the winds of these stars and investigate their effect on ionization structure and wind driving.

Atomic data. For models without hydrodynamic consistency, usually, only the ions that produce measurable signatures and/or sufficient blanketing are included in atmosphere calculations. In the context of B stars, these are H, He, C, N, O, Mg, Al, Si, P, S, and Fe. Following Sander et al. (2018), we also explored elements that could possibly contribute to the wind driving. Besides the previously mentioned elements, we additionally included Na, Ne, Cl, K, Ca, and Ar. However, our test calculations revealed that including these elements did not significantly alter the derived wind properties nor the output spectrum.

In PoWR, Fe is treated conflated with Sc, Ti, V, Cr, Mn, Co, and Ni (see Gräfener et al. 2002, for details). Thus, the abundance of the “generalized iron” is given by the sum of the listed iron-group elements. As Fe is by far the most abundant among them (see, e.g., Asplund et al. 2009), and we do not aim for precise Fe abundance determination in this study, we refer to the “generalized iron” simply as Fe in the rest of this work. The list of elements, ionization levels, and the respective number of transitions are shown in Table B.1.

3 Observational data

We compared our synthetic spectra with archival observed data of the prototypical star ζ¹ Sco in the UV, optical, and IR. The UV observations were obtained using archival data from the Far Ultraviolet Spectroscopic Explorer (FUSE), whose observation PI is T.P. Snow (program ID: P116), and from the International Ultraviolet Explorer (IUE). In the optical, we chose to compare with the spectrum obtained with the Echelle SPectrograph for Rocky Exoplanets and Stable Spectroscopic Observations (ESPRESSO) by the PI A.deCia (program ID: 0102.C-0699) due to its higher resolution among the available spectra. The near-IR spectrum was obtained with X-Shooter at the Very Large Telescope (VLT). We also used spectrum from the Infrared Space Observatory (ISO) obtained via the ISO SWS Atlas (Sloan et al. 2003) in the mid- and far-IR. ζ¹ Sco presents some line variability across the whole spectrum, which we show in Appendix D. As the available observed spectra have an absolute flux incompatible with the retrieved magnitudes, we normalized the optical and near-IR spectra by dividing the observed spectra by a curve obtained from linearly interpolating individual points in the continuum. The UV range was normalized using our elected best model’s synthetic continuum.

The flux-calibrated FUSE, IUE, and ISO spectra were used to build the spectral energy distribution (SED) for constraining the luminosity. Given the available multi-epoch IUE data, we used an averaged spectrum for the SED. Additionally, we also used the magnitudes from the optical to the far-IR: The Johnson-Cousins magnitudes are retrieved from Paunzen (2022) and Zacharias et al. (2004). The Gaia magnitudes used are those listed in Gaia Data Release 3 (Gaia Collaboration 2016, 2023). We adopted the distance from Bailer-Jones et al. (2021). The employed IR magnitudes are the JHK bands from the Two Micron All-Sky Survey (2MASS) and bands from the Wide-field Infrared Explorer (WISE) (both sets retrieved from Cutri et al. 2013). In TableE.1 we list the nominal values of all employed magnitudes.

4 Spectroscopic validation of the model

To create a dynamically consistent atmosphere model that reasonably reproduces the spectrum of ζ¹ Sco, we started with the derived parameters from Clark et al. (2012). We then varied the different stellar properties, initially the luminosity (L), the inner-boundary temperature² (T,), and the stellar mass (M). After finding the best suitable parameter set, we then adjusted the helium mass fraction (X_He), microturbulent velocity (v_urb), clumping, and X-rays. When necessary, the initial parameter set was also readjusted. As we do not aim at a precise abundance determination, we initially adopted the values of Clark et al. (2012) for the He and CNO content. We tested adjustments within their respective error bars, but this did not lead to significant improvements in the spectral reproduction. The elemental abundances of all other metals are fixed to the values of Asplund et al. (2009).

By design, hydrodynamically consistent models have a lower number of free parameters and it is not possible to tweak the stellar and wind parameters separately to improve the spectral fit. Instead, as we discuss in more detail in a related work on very massive stars (Sabhahit et al. 2025), changes in one input parameter can have multiple affects on the spectrum. The main purpose of the current work is therefore not to reproduce every diagnostic in the observed spectrum and fine-tune the stellar parameters, but capture the overall spectral appearance of a BHG like ζ¹ Sco and test if our physical assumptions in this regime are sufficient to describe the wind of an early-type BHG. A more detailed description of our spectral analysis procedure and the employed diagnostics is given in Appendix C.

Fig. 1 Comparison of our final hydrodynamically consistent model (blue line) with the far-UV spectra (FUSE and IUE; gray lines) of ζ1 Sco.

4.1 Ultraviolet spectrum

The overall UV spectrum of ζ¹ Sco remains unchanged over time. We inspected IUE data from 1980 to 1995 and no significant changes in most of P Cygni profiles were observed. Small changes could be seen in the saturation level of Si IV λ1400 and Al III λ1670. This points to a quite stationary global outflow, at least closer to the terminal velocity regime. However, we find that N V λ1240, which in B star models requires X-rays, changes its shape and strength meaningfully. In Fig. 1 we show the comparison between our final PoWR^HD model and observations from FUSE and IUE (average spectrum).

To obtain a spectrum that matches the absorption troughs of the UV P Cygni lines, we applied a maximum wind microturbulence of ξ_max = 180 km s^-1 in the computation of the observer’s frame spectrum. Otherwise, the obtained terminal velocity r_æ from the PoWR^HD solution would not produce wide enough P-Cygni profiles. We discuss this aspect further in Sect. 6.

Despite the general similarities between the UV spectra of BSGs and BHGs (e.g., with respect to the presence of superionization), there are also noticeable differences: A prominent example is the prominent presence of Al II λ1670 and Mg II λ2800 in the few available UV spectra of B1Ia+ and B1.5Ia+ hypergiants³, which in BSGs appears only for B8 and later spectral types (e.g., Walborn & Nichols-Bohlin 1987). This is likely a consequence of the higher wind densities in BHGs. However, we were not able to reproduce these low ionization features.

To reproduce the observed superionization P-Cygni profiles in the UV, we assumed a “hot plasma component”, that is, X-rays with T_X = 0.35 MK and an onset radius of R_X = 500 R.. This means that X-rays only affect the outermost wind layers, where the velocity is similar to r_æ. Our parameters mark the best compromise between reproducing CIV λ1550 and N V λ1240. Other lines were not meaningfully affected. Our final model has L_X/L = −8.5, in line with the upper limit from Berghoefer et al. (1997) of L_X/L < −7.4 (0.2 to 2.4keV). In Fig.2 we show the impact of the X-ray emission on the model spectrum. As evident, the inclusion of X-rays made N V and C IV appear as P-Cygni profiles, a direct consequence of the change in the population numbers. In contrast to BSGs (see Bernini-Peron et al. 2023, 2024), the profiles of CII, Si IV, and Al III in BHGs seem barely affected by X-rays. In other wavelength regimes (optical and IR), the inclusion of X-rays has no meaningful impact.

The X-ray parametrization of our model is very hard to explain under the assumption of shock-induced X-rays. ζ¹ Sco has a slow wind, which would not have enough kinetic energy to produce X-rays via high-velocity collisions in the plasma (Cohen et al. 2014). Moreover, the variability observed in CIV λ1550 and N V λ1240, more intense than in the other lines, would imply a very drastic variation in the X-ray emission over a very extended wind region. Hence, it is possible that the additional ionization source parametrized as X-ray emission simply mimics other physical conditions.

A plausible alternative to explain the superionization is the presence of a “hot component” in the wind, which would not necessarily have X-ray temperatures. Recent multidimensional simulations of O and WR stars further predict a distribution of velocities and temperatures for the wind material (see Moens et al. 2022; Debnath et al. 2024). Although not tested yet, one could speculate that the observed superionization is the manifestation of the hotter part of such a distribution. This might then also explain the different behavior and morphology of CIV 41550 and NV 41240.

Fig. 2 Upper panels: comparison of models’ spectra with and without X-rays (dashed red and solid blue curves, respectively). Middle panel: contribution of Fe III, C IV, and N V for the wind acceleration (black, green, and blue lines, respectively). The lines with thick (thin) crosses indicate the model with (without) X-rays. Lower panel: population number of N V and CIV following the same color and symbol code as the middle panel. in the middle and lower panels, the dashed blue line indicates the photosphere (R2/3) and the solid gray line indicates the critical radius (Rclit). The onset of X-rays (RX) is represented by the dashed brown line.

4.2 Optical

The optical region of BHG spectra is characterized by a rich set of absorption lines, P-Cygni lines, and emission lines. Thereby, it encodes important information about the photosphere as well as the inner layers of the stellar wind (Clark et al. 2012). However, as the Balmer absorption lines are blended with wind P-Cygni features, the determination of surface gravities (log g) via the evaluation of pressure-broadened wings is challenging or in some cases simply impossible. Therefore, one has to rely on higher-order series members, such as He and Hζ – whose widths are slightly underestimated by our models. Additionally, the combined imprint of photospheric and wind features marks a particular challenge in the context of dynamically consistent modeling. The spectral appearance is very sensitive to changes in the stellar properties, which here do not only affect the photospheric features but also the wind parameters and resulting spectral imprints (see also Appendix C).

In Fig. 3, we show a comparison between the observed and synthetic spectrum, focusing on specific diagnostic lines in the optical regime. Given the abovementioned sensitivity, the calculated model atmosphere does an excellent job in reproducing the overall spectral features that are characteristic of ζ¹ Sco and other early BHGs, namely the Balmer lines and some He I lines in P Cygni. The model further reproduces also the bulk of the metal lines across the optical range well. In Fig. D.2, we compare our model with all the available optical spectra from the European Southern Observatory (ESO) archive. It is noticeable that ζ¹ Sco has a lot of line variability, associated with changes in the innermost layers of the wind. While we cannot account for variability explicitly in a stationary approach, our model captures the overall line behavior within the observed range.

4.3 Infrared

The IR and lower energy spectral regions contain information about the wind density ρ(r). Namely, the continuum is dominated by free-free processes, which scale with ρ² (e.g., Barlow & Cohen 1977; Najarro et al. 2009; Rubio-Díez et al. 2022). In the IR, one also finds lines such as He I λ10830 that can be used as an indicator of the terminal velocity in hot supergiants and hypergiants (see Fig. 4). Our model closely reproduces the width of this specific line, indicating that our derived v_æ is realistic.

In contrast to the UV, but similar to the optical, we observe more pronounced line variability in the IR regime, especially in wind-affected H and He lines. This again suggests changes in the photospheric properties as well as in the base of the wind.

In the mid-IR regime, covered by the ISO observations shown in Fig. 5, we find mostly low-energy transitions of H and He. In this range, the vast majority of the predicted emission lines are stronger than what is observed, similar to what we obtain in the near-IR. Nonetheless, the higher order of the Paschen series in the near-IR and the Pfund series in the mid-IR are in good agreement between the model and the observation.

However, we obtain a clear trend of overestimated emission of the Brackett, Paschen, and Pfund lines, indicating too much recombination in the outer wind. This may indicate that the mass-loss rates are slightly too high or that the clumping structure in the outer wind is not well reproduced. As discussed by Najarro et al. (2011) in the context of LBVs, a better fit in the IR can be obtained considering a volume-filling factor that returns to unity in the outer wind. In our specific case, applying this clumping law does not produce a sufficient reduction of the emission to account for the discrepancy. Still, when considering models with lower clumping, we are able to solve the discrepancy, at the cost of spoiling the fits in other regions as another hydrodynamical solution is obtained (see the discussion in Sect. 6.3). This implies that the findings by Najarro et al. (2011) of a decrease in clumping are qualitatively in line with our modeling results, but the specific description of the clumping decrease in the outer wind is insufficient for our particular case.

In contrast to the spectrum of P Cyg (Najarro et al. 1996), we do not observe forbidden lines in ζ¹ Sco, which would have been evidence of circumstellar material.

Fig. 3 Comparison of our final PoWRHD model (blue line) with the observed optical spectrum (ESPRESSO) of ζ1 Sco (gray line).

Fig. 4 Comparison of our best-fit model (blue line) with the near-IR spectral features (X-Shooter; gray lines).

Fig. 5 Comparison of our best-fit model (blue lines) with the mid-IR spectral features (ISO; gray lines).

5 Comparison to previous works

Overall, the properties of our final PoWRHD match relatively well the properties deduced in quantitative spectroscopic analyses by Crowther et al. (2006) and Clark et al. (2012), performed with CMFGEN (Hillier & Miller 1998) and using atmosphere models without hydrodynamical consistency. In Table 1 we list the basic stellar properties we obtained for our final model and compare them to the aforementioned studies. As discussed in Sect. 4 and Appendix C, the stellar mass obtained with PoWR^HD may be underestimated and affected by a certain degeneracy with the assumed clumping, as we discuss further in Sect. 6.3. Nonetheless, our obtained L/M ratio is essentially the same as obtained by Clark et al. (2012) – differing by less than 2%.

Despite some mismatches between the observed and synthetic spectra, the obtained stellar properties are in good agreement with previous determinations in the literature. Namely, the differences from our obtained T_eff and log g compared to literature values are well within the typical errors associated with different methods of quantitative spectroscopy (see Sander et al. 2024).

Despite the high luminosity and dense wind, the Eddington parameter of our model is “only” Γ_e = 0.51. Our finding aligns with what Clark et al. reported for the late BHG Cyg OB2 #12 (B3Ia+), which, despite possibly surpassing the Humphreys- Davidson limit (Humphreys & Davidson 1979), also has Γ_e = 0.40. Their value for ζ¹ Sco is Γ_e = 0.47, very close to our findings. When considering the total Γ_rad in the layers below the critical point $\left({\mathrm{\Gamma }}_{\mathrm{rad}}^{\mathrm{in}}\right),$, which takes into account the total line opacity, we find ${\mathrm{\Gamma }}_{\mathrm{rad}}^{\mathrm{in}}~0.9$. This illustrates that due to the line and bound-free opacities, ζ¹ Sco is actually quite close to the Eddington limit. Our model also shows a small super-Eddington regime in the subsonic region, which might indicate the presence of a turbulent layer (see Debnath et al. 2024) and could possibly explain some of the observed line profile variations. Currently, there is no self-consistent treatment of the resulting turbulent pressure in 1D atmosphere, but it can be included in a parameterized way (e.g., González-Torà et al. 2025). In this work, we adopted a value of v_turb = 14 km s^–1, which is similar to Clark et al. (2012) but actually might be an underestimation. In Appendix C we briefly discuss the impact of varying v_turb in the model calculations.

6 Wind dynamics

6.1 Mass-loss rate and wind acceleration

Our obtained PoWRHD solution has a mass-loss rate of $\stackrel{·}{M}$ = 10^–5.27M_⊙ yr^–1 = 5.3 M_⊙ Myr^–1. This result lies within the range found in previous studies by Crowther et al. (2006), Clark et al. (2012), and Rubio-Díez et al. (2022), albeit with some variations by factors of more than 3 (see Table 2). However, as also evident from Table 2, when we consider $\stackrel{·}{M}\sqrt{f_{\infty }^{-1}},$, we see an agreement within 10% between the different studies for this quantity.

On the other hand, when we consider the transformed mass-loss rates ${\stackrel{·}{M}}_t=\stackrel{·}{M}\cdot f_{\infty }^{-1/2}\cdot \left({10}^3\mathrm{km}{s}^{-1}/v_{\infty }\right)\cdot {\left({10}^6{L}_{\odot }/L\right)}^{3/4},$, introduced by Gräfener & Vink (2013), we obtain bigger disagreements. Compared to Crowther et al. (2006) a difference of a factor of 2 is obtained, whereas for Clark et al. (2012) the difference is about 20%. The reason for our obtained high Mt is primarily linked to our obtained lower terminal velocity.

In Fig. 6, we show the wind acceleration normalized by gravity for different ions in our dynamically consistent wind model for ζ¹ Sco. The acceleration for each Fe ionization stage is shown in the upper panel. In the subsonic regime, Fe IV and higher ionization stages dominate the acceleration. However, toward the sonic point, Fe III takes over and the contribution of Fe IV drops to ~1%. Notably, at the critical point R_crit, the contribution of Fe III to the total radiative acceleration, being comparable with the electron scattering acceleration, is the most important. Hence, this is the ion that effectively launches the stellar wind. In the outer wind, the Fe III also maintains dominance.

This behavior is similar to the hydrodynamically consistent model of the early-BSG donor of VelaX-1 (B0.5Ia) computed by Sander et al. (2018) using the same code. In their models, Fe III is the ion that effectively launches the wind whereas Fe IV is secondary. However, in their case, the contribution of Fe IV is about 10% near R_crit due to the higher temperature in the B0.5Ia star. In the case of ζ Pup (O4Iaf⁺), analyzed and discussed in Sander et al. (2017) following the same methodology, Fe is also the leading element but the Fe IV ionization stage is the dominant source of acceleration at the critical radius.

The leading acceleration ion is important in the context of the so-called bi-stability jump. The idea of a bi-stable solution for the mass-loss rates stemmed from Pauldrach & Puls (1990), who computed the wind solution for P Cyg (B1Ia+/LBV), whose temperature was estimated ~21 kK in their study. Later on, studies such as Lamers et al. (1995) and Markova & Puls (2008) found drops in the terminal velocity roughly at that temperature, reinforcing the scenario of a sudden change in the wind properties due to the recombination of Fe IV to Fe III within a narrow range of stellar temperatures around ~22kK. Utilizing the observed terminal velocity change, this led to a significant jump in the mass-loss rates in the Monte Carlo calculations by Vink et al. (1999), which also entered the widely used Vink et al. (2001) mass-loss recipe including its recent update (Vink & Sander 2021). For stars with winds driven (mainly) by Fe III, the massloss rates obtained with the scheme from Vink et al. (1999) are about an order of magnitude larger than those of winds driven by Fe IV. A similar conclusion was reached in the analysis by Petrov et al. (2016). However, in the past two decades, empirical studies on BSGs across that region (e.g., Crowther et al. 2006; Bernini-Peron et al. 2024; Verhamme et al. 2024) and new theoretical calculations (cf. Krtička et al. 2021; Björklund et al. 2023) challenge the scenario of a sharp M increase associated with the dominance of Fe III as the main wind driver. Our modeling results place ζ¹ Sco on the cool side of the (hot) bi-stability jump. Ignoring the second bi-stability jump, the recipe from Vink et al. (2001) yields log $\left(\stackrel{·}{M}\left(M_{\odot }{\mathrm{yr}}^{-1}\right)\right)=-4.42,$, while our obtained value of −5.27 is almost an order of magnitude lower. Interestingly, both the descriptions from Björklund et al. (2023) and Krtička et al. (2024) yield results that are significantly too low, namely −5.93 and −5.94, respectively. Assuming our mass is actually an underestimation, the result from Björklund et al. (2023) would even be lower. However, neither the studies by Björklund et al. (2023) nor Krtička et al. (2024) were designed for BHGs. Nonetheless, we can conclude that despite all similarities to BSGs, the mass-loss rates of BHGs are higher than recently predicted, though still considerably lower than estimated by Vink et al. (2001).

Fig. 6 Wind-driving contribution of different ions to the total wind acceleration, ${\mathrm{\Gamma }}_{\mathrm{rad}}^{\mathrm{total}}$. The upper panel focuses on the contribution of different Fe ions. The acceleration by electron scattering (Γes) and by H I are shown. The lower panel shows the breakdown of the contributions of CNO and Si ions. Fe III and the electron scattering acceleration, the most important contributors to wind launching, are also shown for comparison. The vertical lines represent the photospheric, sonic, and critical radii (dotted gray, dashed gray, and thin black lines, respectively). The dotted brown line marks the point where X-ray emission is switched on.

6.2 Wind velocity field

BHGs and LBVs are characterized by their slow winds, which can often be slower than 300 km s⁻¹ and, from most of the β-values required in earlier spectral analysis, are expected to have a quite shallow velocity gradient. These characteristics make it hard to explain the superionization by shock-induced X-rays and defy the typical m-CAK picture of a velocity field with β values around and below unity (e.g., Curé & Araya 2023).

Wind terminal and turbulent velocities. The stellar properties used as input for the hydrodynamically consistent atmosphere model yield a terminal velocity of v_æ = 300 km s⁻¹. This value is about 90 km s⁻¹ lower than what previous non- hydrodynamically consistent spectroscopic analyses find (e.g., Crowther et al. 2006; Clark et al. 2012). As mentioned in the discussion of the UV results (cf. Sect. 4.1), require a microtur- bulent velocity of ξ_max = 180 km s⁻¹ in the formal integral to reconcile the model with the observation. Crowther et al. (2006) adopted ξ_max = 50 km s⁻¹, whereas Clark et al. (2012) considered a value of ξ_max = 75 km s⁻¹. While their values are lower than ours, they are higher than what is typically used in quantitative spectroscopy of OB stars in general (ξ_max ~ 0.1v_æ) and higher than the velocity dispersion obtained in line-deshadowing instability (LDI) hydrodynamic simulations of a BSG wind by Driessen et al. (2019). While our derived terminal velocity might be slightly underestimated, this generally higher microturbulence could point to a more fundamental difference between the winds of BSGs and BHGs, potentially related to the previously discussed turbulent pressure arising from the proximity of the hot iron bump to the Eddington limit.

Validity of the β law approximation. To first order, the resulting velocity field in the supersonic regime is qualitatively similar to a β law with a value of β ≳ 3 as this would lead to a comparably gradual acceleration of the wind. However, important differences appear when examining the field in detail, especially the layers further in the wind. Deviations from the β law are also found in previous studies with hydrodynamically consistent models for O-type stars and BSGs (e.g., Sander et al. 2017, 2018; Björklund et al. 2021).

In Fig. 7, we compare the velocity fields of a hydrodynamically consistent model with a set of models that use a hydrostatic stratification connected to a β law, hereafter referred as “classical PoWR models”. All the models have the same stellar properties. Their work ratio⁴ is also about unity, implying a global, but not necessarily local consistency with respect to the force balance. The classical models differ among themselves in terms of how the connection between the (quasi- )hydrostatic and the wind profile is described. We calculated one model with a smooth subsonic-to-supersonic gradient connection ((dv/dr)_hys = (dv/dr)β) and models where the velocity of the connection point v_con is fixed to different fractions of the sound speed v_sonic (0.9, 0.7, and 0.5). The choice of β = 3.3 was obtained by employing nonlinear least squares to fit a β law to the obtained hydrodynamically consistent solution.

The comparison reveals that the subsonic and transition regions agree remarkably well between the hydrodynamically consistent PoWRHD and the classical PoWR models with v_con = 0.9 v_sonic. Their velocity gradients also match considerably well. Such an agreement is to be expected for the inner subsonic layers as long as the atmosphere is nearing (quasi-)hydrostatic equilibrium, which is ensured even for PoWR models without full dynamical consistency (Sander et al. 2015). However, when forcing a smoother connection by imposing equal gradients at the transition point (cf. Fig. 7), the differences increase. The smoothconnection model also produces spectra that differ more from the observations (see, e.g., the last panel of Fig. 7, where the Hα of the smooth model is the weakest). When inspecting this model, we find that the resulting connection point is significantly more inward than in the “abrupt” transition models and connects to the hydrostatic gradient during its large increase. Considering both the hydrodynamic and the 0.9 v_sonic models, we can conclude that the additional subsonic peak in the gradient seems to mark an important feature for correctly predicting the spectrum.

In the supersonic regime, the velocity fields start to differ more distinctively, which is also noticeable in the velocity gradient. Up to ~10 R_*, the wind accelerates faster than the β law profile but afterward the trend reverses. This result suggests that winds of BHGs have a shallower acceleration than the typical β laws applied to OB stars (i.e., 0.5 > β > 1.0). This is in alignment with a series of investigations in the past decades that systematically infer β ≳ 2.0 for late BSGs and BHGs. Many studies modeling LBVs, which have denser winds, also find β ~ 3 (e.g., Najarro et al. 2009; Kostenkov et al. 2020) – though values of β~1 have also been found (e.g., Najarro et al. 2009; Maryeva et al. 2022).

We also tested models with different β parameters. We do not find meaningful differences in the UV spectral appearance as most of the lines are already saturated. However, the lines in the optical and near-IR are significantly affected, as the velocity field around the connection between the quasi-hydrostatic and wind regime is considerably altered. High-velocity gradients (β < 2.3) in the inner wind layers produce atmospheres with weaker wind features. This underlines the need for β > 2.0 for such stars. Additionally, even if adopting extremely high values of β (e.g., β > 10), there is no good match with the wind velocity profile of the hydrodynamic model at r ≳ 10R_*. One possibility to circumvent this problem is to adopt a double-β law velocity profile, which can account for an inner acceleration different from the acceleration in the outer parts of the wind. This has for example been done for the hydrodynamical solutions in Björklund et al. (2021). Using the same approach, we obtain a parametric profile overall more similar to our hydrodynamical model up to ~4 R_* (see Fig. 8). Still, we obtain rather high values of the two parameters, namely β₁ = 2.7 and β₂ = 2.0. In contrast, Björklund et al. (2021) obtained a value of only 0.8 for their inner β₁ -parameter in their O star models describing a thinner wind regime.

As discussed for WR stars in Lefever et al. (2023), there may be degeneracies between the β and effective temperature for optically thick winds. This is especially an issue if the whole spectrum is formed in the (outer) wind, which is the case for a fraction of the Galactic WR population (e.g., Hamann et al. 2019; Sander et al. 2019) or objects like η Car (e.g., Hillier et al. 2001). While also some BHGs, such as the more extreme P Cyg, might be affected by this degeneracy, ζ¹ Sco seems far from that regime. Therefore, any choice of the β value should not significantly influence its temperature. Still, different β values alter the connection regime, and thus the density structure, which would impact the many diagnostic lines formed at the base of the wind in typical BHGs as discussed above.

Venero et al. (2024) modeled the velocity field and Hα of late BSGs (and the BHG, HD 80077, B2Ia+) using approximative hydrodynamical solutions within the “m-CAK formalism”⁵ (see also Curé et al. 2011; Venero et al. 2016; Vink 2018). They found a steeper velocity gradient in the inner wind compared to a β law profile and a flatter gradient outward. Qualitatively, we find a similar behavior of the velocity field of our early BHG. Their obtained velocity field for HD 80077 belongs to the δ-slow solution regime. In fact, they found that the BHG wind is only viable within this topological domain.

Another possible solution of the m-CAK hydrodynamics in the B star domain is the so-called Ω-slow solutions (see Curé & Rial 2004; Curé & Araya 2023). This type of solution emerges in fast-rotating objects and is also associated with low v_æ, as observed in BHGs. However, this type of solution would only occur for stars rotating ≳75% of the critical velocity⁶. ζ¹ Sco (and BHGs in general) are not fast rotators and far from that threshold - namely, ζ¹ Sco rotates only about 15% critical. Moreover, the obtained PoWR^HD velocity profile does not have a gradient in the inner wind as steep as predicted by the Ω-slow solution.

Fig. 7 Comparison of hydrodynamically consistent, pure β law profiles with classical PoWR models with hydrostatic profiles. Top panel: velocity profile as a function of the stellar radius. The gray lines indicate the pure theoretical β law profiles with β = 2.9, 3.3, and 3.9 (dot- dashed, dashed, and solid lines, respectively). The bold blue line is the hydrodynamically consistent PoWRHD model, whereas the green (continuous velocity gradient connection), purple (forced connection at v = 0.5vsound), and red (forced connection at v = 0.9vsound) profiles indicate the classical (β law + hydrostatic integration) PoWR models with β = 3.3. The light blue vertical line indicates the critical point of the PoWRHD model. Second panel: same as the top panel but dividing all the velocity profiles by the hydrodynamical solution, and hence highlighting where velocity profiles differ more. Third panel: comparison of the velocity gradient (dv/dr) profiles. Bottom panel: Spectroscopic comparison between the PoWR models (following the same color code). The observed optical spectrum (ESPRESSO; light gray) is also shown for reference.

Fig. 8 Comparison of the hydrodynamic consistent solution (blue thick line) with fits of a double β profile (black thick line) as presented in Björklund et al. (2021). We also show for reference the best β model from Fig. 7 as the red line.

6.3 Clumping

In this work, wind inhomogeneities are treated in the “microclumping” approximation, assuming that clumps are optically thin. Given that the winds of BHGs can be optically thick in certain spectral lines, this approximation neglects the effects of so-called macroclumping (Oskinova et al. 2007; Sundqvist & Puls 2018), that is, the effect of optically thick clumps on the spectral diagnostics. So far, PoWR only treats macroclumping in the formal integral, which implies that it does not affect our derived hydrodynamic solutions. While calculations for a BSG in Bernini-Peron et al. (2023) showed that macroclumping has an effect on the spectrum, we refrain from introducing it in this study, as we do not aim for precise reproduction of all spectral lines.

In atmosphere models without dynamical consistency, the volume filling factor f (r) (cf. Eq. (1)) is adjusted to reconcile mass-loss rate diagnostics obtained from resonance and recombination lines. In hydrodynamically consistent models, however, changes in the clumping factor can also directly impact the derived mass-loss rate. In particular changes of f (R_crit) affect the resulting opacity (via a change in the wind density; see, e.g., Bouret et al. 2012) and thus the inferred wind dynamics. Therefore, we investigated the impact of different clumping parameters both for wind driving and the spectral appearance of our early-BHG model.

Amount and onset of clumping. To obtain spectra that resemble the observations of ζ¹ Sco, we required a low clumping value (f_cl = 1.5 or f_∞ = 0.66) with v_cl = 5. This f_∞ value is similar to those inferred by Rubio-Díez et al. (2022) via fitting the IR and radio spectral regions, providing additional evidence that the winds of BHGs are not very clumpy. Our findings further align well with the results from time-dependent hydrodynamical simulations by Driessen et al. (2019). They found that BSGs on the “cool side” of the bi-stability jump should show way less clumping than their O and early-B counterparts. Notably, this result differs from the spectroscopic analysis by Clark et al. (2012), who obtained f_∞ = 0.06.

We explored different characteristic velocities, v_cl, for clumping, including subsonic values. Changing f_∞ while keeping v_cl mainly affects the amount of clumping in the outer atmosphere and thus also changes mostly the terminal velocity. However, even for the same value of f_∞ different values of v_cl affect the amount of clumping in the inner layers, especially at the critical point. In dynamically consistent models, this then impacts the wind driving as well as the spectral formation at the photosphere (see also Sabhahit et al. 2025).

In Fig. 9, we illustrate the change in the position of the critical point different v_cl. Models with v_cl ~ 5 km s⁻¹ or lower produce a clumping stratification where the filling factor at the critical point f (R_crit) is essentially identical to f_∞ On the other hand, models with higher 3cl have reduced or hardly any clumping at the critical point. Spectroscopically, these models behave more like the smooth wind model and do not predict H and He lines as P Cygni profiles, hence indicating a thinner wind.

In Clark et al. (2012), their best-fit model also has a clumping onset of about ~60% of the sound speed. This aligns well with the spectral analysis done by Najarro et al. (2009) on LBVs in the Quintuplet cluster (Pistol Star and FMM362) using CMFGEN. They also obtained a low onset of clumping (v_cl ~ 2 km s⁻¹). However, both Clark et al. (2012, for ζ¹ Sco) and Najarro et al. (2009, for the LBVs) inferred a much more clumped wind with maximum volume filling factors of 0.06 and 0.08, respectively.

Interestingly, when we consider a model whose clumping starts from the inner boundary we find that the wind becomes again slightly optically thinner, associated with a lower $\stackrel{•}{{\rm M}}$. This may suggest that a variation of clumping within the subsonic regime is important to form the characteristic spectrum of ζ¹ Sco, and perhaps other early BHGs.

Spectral impact of clumping. In Fig. 10, we compare models with different (maximum) clumping factors while leaving all other parameters untouched. More clumped winds produce spectra with more (pronounced) P-Cygni features, a consequence of the resulting higher mass-loss rates. On the other hand, if clumping is switched off, most of the optical lines remain completely in absorption, yielding a BSG-like spectrum.

The differences in wind density arising from different clumping are also reflected in a change of the critical point. As evident from Fig. 11, the models with smoother winds (f_∞ = 0.73–1.0) have R_crit > R_2/3 ≡ R(τ = 2/3) whereas models with higher clumping (f_∞ = 0.35–0.2) have R_crit < R_2/3, meaning that the wind is launched below the formation of the observed spectral continuum. Interestingly, when considering progressively more clumped winds (f_∞ → 0), η reaches a peak close to 0.5 for f_∞ ~ 0.5, but then decreases again. The spectra of these models with a high clumping, however, do not look like a typical BHG anymore but instead resemble the spectral appearance of known LBVs in the quiescent phase – also referred to as “P-Cygni supergiants” from a spectroscopic point of view, such as PCyg, R81 and AG Car (Clark et al. 2012). The eventual decrease in η is the result of a rapid decrease in the terminal velocity v_∞ while $\stackrel{·}{M}$ still increases (see the first panel of Fig. 10). Assuming the same energy budget, a denser wind would in principle move more slowly, but this reciprocal behavior between v_∞ and M is notably only seen for very high or very low clumping in our model sequence, while for intermediate clumping values (0.3 ≲ f_∞ ≲ 0.8), v_∞ increases with $\stackrel{·}{M}$. This behavior switch should not be expected generally but depends significantly on the specific parameter regime and the choice of the clumping stratification and onset.

Changes in the velocity gradient. When inspecting the velocity gradient structure for the models with different clumping, it is noticeable that the models with the highest and the lowest clumping values have a single peak on their respective gradients, which correlates well with the position of their critical point - see Fig. 12. Curiously, our best model of ζ¹ Sco has a double-peak structure, which as discussed in Sect. 6.2 is necessary to get the proper spectral appearance of BHG winds as a transition between BSGs and quiescent LBVs.

Clumping-mass degeneracy. Sabhahit et al. (2025) found a certain degeneracy between the clumping at the critical point and the inferred stellar masses when analyzing WNh stars of ~50kK with hydrodynamically consistent models. For a given reference model, lowering (increasing) clumping at the critical point and lowering (increasing) the mass will cause the model to keep a similar spectral appearance and also yield similar wind properties.

For our studied BHG regime, we find a similar effect when comparing a clumped (f_∞ = 0.67, our final model) with a smooth-wind model. This comparison is illustrated in Fig. 13, where we note a different behavior between the smooth-wind and clumped model regarding the driving and spectral appearance. Namely, the smooth-wind model has lower mass-loss rates, higher terminal velocities, and different wind-driving morphology – all of which translates to a spectrum with fewer emission and P-Cygni lines. However, by decreasing the mass (only by 5%) we obtain a model that is very similar to the clumped model in terms of structure, wind properties, and spectral appearance. The smooth-wind model thereby also sets a lower limit of 36 M_⊙ for the derived mass of ζ¹ Sco, which coincides with the value from Clark et al. (2012), albeit for a different luminosity.

Fig. 9 Variation in the clumping structure and output wind properties according to different characteristic onset velocities, vcl. Left panels: photosphere, i.e., where the Rosseland-mean optical depth, τ, is 2/3, as marked by the diamond symbol with its vertical line. Right panels: corresponding output spectra compared to the observed spectrum of ζ1 Sco.

Fig. 10 Upper panel: Change in the mass-loss rates, $\stackrel{·}{M}$, and terminal velocities, v∞, as a function of different values of the clumping volume filling factor, f∞. The colors trace f∞ purple indicates the homogeneous-wind model (f∞=1.0) and red is most clumped wind tested (f∞ = 0.2). The final model is represented by the black diamonds. Middle panels: changes in the wind efficiency parameter (η) and Teff as of function of f∞. Lower panels: different lines’ spectral morphology changes as a function of different clumping values. The color scheme is the same as the upper panels, with the thick black line indicating the best-fitting model spectrum.

7 Are BHG winds in the transition regime?

O stars with a spectrum dominated by P Cygni-type profiles in the optical regime are close to the onset of the regime of optically thick winds, where the wind onset point has an optical depth larger than unity. Vink & Gräfener (2012) demonstrated that in the case of $\stackrel{·}{M}$, the mass-loss rate can be directly inferred via $$\stackrel{·}{M}=\frac{L}{v_{\infty }c}\equiv {\stackrel{·}{M}}_{\text{trans}}.$$(2)

where the τ_F is the flux-weighted optical depth, including line opacities, at the critical point. This is called the “transition massloss rate”, ${\stackrel{•}{M}}_{\text{trans}}$⁷. Detailed modeling with both Monte Carlo (Vink & Gräfener 2012) and PoWRHD (Sabhahit et al. 2023) has shown that often η = 1 does not exactly correspond to τ_F = 1 and a correction factor f has to be applied, typically with f ≈ 0.6 (Vink & Gräfener 2012; Sabhahit et al. 2023).

The concept of this “transition mass-loss rate” was successfully used to predict the M of very massive Of/WN-type stars (“slash stars”; see Sabhahit et al. 2022, 2023). In the optical, these stars show photospheric profiles mixed with clear wind features (see, e.g., Crowther & Walborn 2011). Thus, the question that may arise is whether BHGs are cooler counterparts of the slash stars.

While BHGs are morphologically similar to the slash stars in terms of the presence of wind features contaminating the photospheric lines, our models suggest that early BHGs have a very low η, < 0.1, which is rather similar to that of BSGs. For τ_F, we derive a value of 0.4 for our preferred model, meaning that ζ¹ Sco is not that close to the transition regime. To quantify this, we calculated ${\stackrel{·}{M}}_{\text{trans}}$ according to Eq. (2) using the luminosity and terminal velocity of ζ¹ Sco. The computed value of ${\stackrel{·}{M}}_{\text{trans}}$ = 10^−4.37 is almost an order of magnitude above our determined value of 10⁻⁵-²⁷ (and actually very similar to the derived value when using the recipe from Vink et al. 2001).

When considering models with higher clumping (cf. Sect. 6.3), we observe an increase in η, which approaches ~0.5. However, when the wind density reaches a certain point, the terminal velocity starts to decrease, reducing η as a consequence. While further studies with dynamically consistent atmospheres for individual BHGs and P Cygni-like supergiants are necessary, it seems currently unlikely that typical BHGs reach the transition regime and thus their mass-loss rates cannot be determined via the ${\stackrel{•}{M}}_{\text{trans}}$-formalism.

Fig. 11 Upper panel: changes in the position of the critical (Rcrit; dark gray diamonds), sonic (Rsonίc; light gray squares), and photospheric (R2/3; colored stars connected by the blue line) radii as a function of different values of the clumping volume filling factor, f⊙ The bold black symbols indicate the best-fit model. Lower panel: same as the upper panel but showing the sonic and critical points in terms of the Rosseland-mean optical depth, τ. The blue horizontal line indicates the photosphere, where τ = 2/3.

Fig. 12 Velocity gradient (dv/dr) vs. Rosseland-mean optical depth (τ). Each model is shifted vertically by an arbitrary amount for visualization purposes. The diamond symbols indicate the critical point, whereas the blue vertical line indicates the photosphere (τ = 2/3). The thick black curve and diamonds refer to the final model (f∞ = 0.66).

Fig. 13 Degeneracy between clumping and stellar mass. In all panels, the black lines indicate the final model, the light blue lines the model without clumping, and the dark blue lines the model without clumping with a reduced mass. Upper panel: radiative acceleration vs. Rosseland optical depth. The total radiative acceleration, ${\mathrm{\Gamma }}_{rad}^{total}$, is indicated by the thick line with dots, and the thinner lines below indicate the contribution of Fe III and Fe IV. The thin vertical lines indicate the position of the critical points. Middle panel: same as the upper panes but showing the contribution of the main ions of C, N, S, and Si. Lower panels: comparison between the models’ output spectra. Mass-loss rate units are M⊙ yr−-1.

8 Conclusions

We used PoWR^HD to produce the first hydrodynamically consistent atmosphere model of ζ¹ Sco - the BHG with currently the most extensive spectroscopic coverage, considered a prototype of its kind. We successfully reproduce most of the characteristic spectral features of BHGs from the UV to the IR, which indicates that the atmosphere and wind structure are well described by our model. From our modeling efforts and observational comparisons, we derive the following conclusions:

• We derive a mass-loss rate for ζ¹ Sco of 10^−5.27 M_⊙ yr⁻¹, which is about one order of magnitude below the prediction from Vink et al. (2001) but almost one order of magnitude higher than predicted by Björklund et al. (2023) and Krtička et al. (2024). The radiative acceleration in the wind is dominated by Fe III. This is also true at the critical point, which is located close to, but still above, the photosphere. The derived ionization and acceleration is coherent with objects located on the cool side of the so-called bi-stability jump;

• Despite BHGs being spectroscopically similar to Of/WNh in terms of the presence of wind profiles in Hβ (and higher- energy Balmer lines), their wind efficiency is too low to place them in the “thick wind” regime. This shows that for such low-temperature stars, the optical spectrum is capable of producing wind features while remaining relatively optically thin. This prevents the application of the transition mass-loss rate concept (Vink & Gräfener 2012) in this regime;

• The radiative acceleration produces a velocity profile that differs from the widely used β law and has a subsonic as well as a supersonic peak in the velocity gradient dv/dr. When restricted to β laws, our derived structure can be partially reproduced by using β = 3.3 connected to a hydrostatic regime at 90% of the sound speed;

• When deriving the wind dynamics, the amount of clumping and its onset play an important role in the wind launching and in shaping the optical spectrum. Lower clumping values and a subsonic onset are necessary to reproduce the main spectral features of ζ¹ Sco. On the one hand, in the absence of clumping, the BHG features are lost and the spectrum instead resembles that of a BSG. On the other hand, when clumping is increased further, the spectral appearance of quiescent LBVs is eventually reached. In addition, our BHG model shows the presence of a subsonic super-Eddington layer, which will likely cause radiatively driven turbulent pressure and could also help trigger the necessary clumping;

• The terminal velocity obtained from the self-consistent wind solution is lower than obtained in previous quantitative spectroscopy studies with prescribed velocity fields. To reproduce the observed UV P-Cygni absorption troughs,

we needed to employ a high microturbulent velocity in the formal integral. Future investigations of more targets will be required to test whether such turbulent velocities are reasonable in this regime.

Acknowledgements

We thank Joris Josiek, Wolf-Rainer Hamann, Helge Todt, and Maurício Gómez-González for their insightful comments and discussions on the modeling process as well as in observational aspects of this work. MBP, AACS, and RRL are supported by the German Deutsche Forschungsgemein- schaft, DFG in the form of an Emmy Noether Research Group - Project-ID 445674056 (SA4064/1-1, PI Sander). GGT is supported by the German Deutsche Forschungsgemeinschaft (DFG) under Project-ID 496854903 (SA4064/2-1, PI Sander). GGT further acknowledges financial support by the Federal Ministry for Economic Affairs and Climate Action (BMWK) via the German Aerospace Center (Deutsches Zentrum für Luft- und Raumfahrt, DLR) grant 50 OR 2503 (PI: Sander). FN acknowledges support by PID2022-137779θB-C41 funded by MCIN/AEI/10.13039/501100011033 by “ERDF A way of making Europe”. GNS and JSV are supported by STFC (Science and Technology Facilities Council) funding under grant number ST/Y001338/1. DP acknowledges financial support from the Deutsches Zentrum für Luft und Raumfahrt (DLR) grant FKZ 50OR2005 and the FWO junior postdoctoral fellowship No. 1256225N. This work has made use of data from the European Space Agency (ESA) mission Gaia (https://www.cosmos.esa.int/gaia), processed by the Gaia Data Processing and Analysis Consortium (DPAC, https://www.cosmos.esa.int/web/gaia/dpac/consortium). Funding for the DPAC has been provided by national institutions, in particular the institutions participating in the Gaia Multilateral Agreement. Based on data obtained from the ESO Science Archive Facility with DOIs, https://doi.org/10.18727/archive/21, https://doi.org/10.18727/archive/33, https://doi.org/10.18727/archive/50, and https://doi.org/10.18727/archive/71

References

Appendix A OB hypergiants in the Local Group

In Table A.1 we provide a list of stars classified as BHGs (and O hypergiants) in the Local Group. In total, we counted 35 objects. The list also include objects that did not formally receive the “Ia+” but show the spectroscopic characteristics. In this table, we sought to exclude stars confirmed LBVs such as PCyg and HD 168607.

Appendix B Atomic data

In Table B.1 we list the elements and ions with the number of levels and line transitions considered in our atmosphere models. For Fe, the levels and line numbers refer to superlevels and superlines (see Gräfener et al. 2002, for more details).

Appendix C Obtaining the final model

We started our hydrodynamically consistent PoWR^HD modeling with a model based on the properties of the best-fit CMFGEN model from Clark et al. (2012). While their study currently marks the most precise spectral fitting of ζ¹ Sco, employing their basic parameters in a dynamically consistent model does not automatically yield the desired spectroscopic similarity. This was expected because not only are there differences between PoWR and CMFGEN, but there is an inherent coupling of stellar and wind properties in PoWR^HD, which affects the derived spectrum.

Consequently, we varied the stellar properties aiming to alter the wind parameters and improve the similarity with the observed UV, optical, and IR spectra. The primary focus is to reproduce the main BHGs spectral signatures (van Genderen et al. 1982), namely the P-Cygni profiles (full or partial) of Hβ, Hγ, Hδ, and He I lines. We further sought to reproduce metal lines in pure absorption as they inform us about the star’s more fundamental properties and adequate boundary conditions (e.g., temperature, mass, and photospheric turbulence).

Fig. C.1 Comparison of the hydrodynamically consistent model of ζ1 Sco (gray line) with the star’s SED and flux-calibrated UV and mid- IR spectra. The UV spectra are retrieved from FUSE and IUE (fuchsia, purple, and blue lines). The mid-IR spectrum was obtained with ISO (brown line). The cyan triangles are the flux points of the UBVI photometry. The orange diamonds are from the 2MASS JHK photometry, and the red squares are from WISE mid- and far-IR photometry. The inlet plot shows the millimeter/radio regime; the black crosses are the flux points obtained with the Very Large Array (VLA). The truncation in the synthetic spectrum is an artifact of the code output. See Appendix E.

To obtain the luminosity, we aimed to reproduce the SED from the far-UV to the millimeter (or radio) regime. We applied the extinction law by Cardelli et al. (1989) considering initially the reddening and extinction parameters from Clark et al. (2012), namely E(B − V) = 0.66 and R_v = 3.3. To obtain a slightly better fit, these were adjusted to the slightly higher E(B − V) = 0.70 and a “standard” R_v = 3.1. In Fig. C.1 we show the resulting fit.

The final model nicely reproduces the observed SED across the whole considered range, indicating a well-constrained luminosity. However, there are remaining offsets of about 0.1 dex in the far-UV regime (IUE-sws) and ~0.2 dex, in the near-UV regime (IUE-lws). While this could potentially imply an underestimation of the extinction, higher values did not improve the overall SED fit.

To determine the effective temperature T_eff, we sought to find the best compromise between Si II 24128,32, Si III 24550,69,76, Si IV λ4090,4116, and Mg IIλ4481. Other lines such as O I λ7772, and He I lines were used as secondary diagnostics. We notice that while Si III λ4550, 65, and 75 are well fitted, SiIV λ4090,4116 is overestimated while Sill λ4128,32 is underestimated, hinting to an overestimation of T_eff. We also observe a similar result for the ionization balance of O II and O I, where O II lines (e.g., O II 24070 and O II 24590,95) seem slightly overestimated, but overall well reproduced, while O I λ7772 is underestimated. Yet, as illustrated in Fig. C.2, lower choices of T_* would have worsened the spectral reproduction due to the consequences on the hydrodynamic wind parameters coupled to the stellar parameters. As depicted, a slight change in T_* (ΔT_eff ~ 1 kK) causes the M to increase by a factor of 2. As log L and M are kept fixed, decreasing(increasing) T_* is equivalent to increas- ing(decreasing) the stellar radius – via the Stefan-Boltzmann relation –, thereby reducing the surface gravity.

Fig. C.2 Resulting variation in the spectral features and wind properties when varying the T* input for PoWRHD calculations. The observed spectrum is shown in black.

To fix the surface gravity, we relied on He and Hζ as they were the cleanest Balmer lines without meaningful wind features - even though models with higher mass-loss rates eventually show wind features in this line as also illustrated in Fig. C.2. In general, we could not obtain models with the wings as broad as the observed spectra, which indicates that the obtained surface gravity, and likely mass, are still too low for this particular star. Due to the inherent coupling between wind and stellar parameters, increasing the mass would spoil the agreement of the main characteristics of ζ¹ Sco.

In Fig. C.3 we illustrate the change in the model’s properties and spectrum by varying the v_turb. Values lower than v_turb ~ 14 km s^–1 yield similar wind properties as output, but weaken the absorption lines, including temperature and gravity diagnostics (e.g., Si and H lines). Since in general our final model, whose v_turb = 14 km s^–1, underestimates the absorption features of most of the spectral lines, it is possible that v_turb ~ 15 km s^–1 underestimates the microturbulence in the star. Moreover, the presence of super-Eddington layers and results from multidimensional simulations of hot stars (albeit in the context of O stars; Debnath et al. 2024) point in that direction.

Fig. C.3 Resulting variation in the spectral features and wind properties when varying the vturb input for PoWRHD calculations. The observed spectrum is shown in black.

Still, an increase in the v_turb beyond ~25 km s⁻¹ would cause meaningful increase in the $\stackrel{•}{M}$ and v_∞, as well as an increase in the absorption features. This value of v_turb is higher than microturbulence values usually used in classical modeling of BSGs (e.g., Crowther et al. 2006; Haucke et al. 2018). We also notice that models with high turbulence have lower T_eff (linked to a higher $\stackrel{•}{M}$), and regarding many wind lines behave similarly to cooler models, despite having the putative temperature diagnostics SiIII and Si IV lines stronger. This illustrates some degeneracies that may arise between certain parameters, which also highlights how challenging it is to obtain final spectral fitting based on the stellar properties within the hydrodynamically consistent framework. In that context, having a consistent treatment of the turbulence would lift such a degeneracy, as v_turb would also emerge from the input physics.

Appendix D Spectral variability

A thorough examination of the spectral variability of ζ¹ Sco, and their causes is out of the scope of this study. Nonetheless, to give a general impression of the observed variability, we present the most relevant diagnostic lines in the UV (Fig. D.1), optical (Fig. D.2) and IR (Fig. D.3). For reference, we also compare with our final PoWR^HD model. More information on the observational dataset can be found in the supplement material at Zenodo⁸.

In the UV range, we notice only mild variability in the P-Cygni profiles. N V λ1240 is the line that presents the most variable behavior. In the optical, we notice the bulk of the variability happening in the emission component of the P-Cygni lines, especially the Balmer lines. However, we noticed a distinctly strong absorption with a secondary component in the HARPS spectra from the year 2006. The metal lines do not suffer much variation. In the X-Shooter IR spectra, which in general are much noisier than the optical, we do not observe meaningful variation, except for the emission components of the He I λ10830 and Paschen-Hβ. We also see variations on the NII λ16164,256, which switches between emission (reproduced by our PoWR^HD model) and absorption.

Fig. D.1 IUE spectra of ζ1 Sco. The blue line is the observation, and the black line is our final PoWRHD model. The shaded area marks the interval of ±500 km s−1 around the reference wavelength.

Fig. D.2 Variability in the spectral features in the optical of ζ1 Sco. Each column refers to a different instrument.

Fig. D.3 Comparison of multi-epoch X-Shooter IR observations of ζ1 Sco.

Appendix E Magnitudes used in the SED

Table E.1 provides an overview of all the photometry used to derive the observed SED for ζ¹ Sco. For each magnitude, the source is listed as well.

All Tables

All Figures

	Fig. 1 Comparison of our final hydrodynamically consistent model (blue line) with the far-UV spectra (FUSE and IUE; gray lines) of ζ1 Sco.
In the text

Fig. 2 Upper panels: comparison of models’ spectra with and without X-rays (dashed red and solid blue curves, respectively). Middle panel: contribution of Fe III, C IV, and N V for the wind acceleration (black, green, and blue lines, respectively). The lines with thick (thin) crosses indicate the model with (without) X-rays. Lower panel: population number of N V and CIV following the same color and symbol code as the middle panel. in the middle and lower panels, the dashed blue line indicates the photosphere (R2/3) and the solid gray line indicates the critical radius (Rclit). The onset of X-rays (RX) is represented by the dashed brown line.

In the text

	Fig. 3 Comparison of our final PoWRHD model (blue line) with the observed optical spectrum (ESPRESSO) of ζ1 Sco (gray line).
In the text

	Fig. 4 Comparison of our best-fit model (blue line) with the near-IR spectral features (X-Shooter; gray lines).
In the text

	Fig. 5 Comparison of our best-fit model (blue lines) with the mid-IR spectral features (ISO; gray lines).
In the text

Fig. 6 Wind-driving contribution of different ions to the total wind acceleration, ${\mathrm{\Gamma }}_{\mathrm{rad}}^{\mathrm{total}}$. The upper panel focuses on the contribution of different Fe ions. The acceleration by electron scattering (Γes) and by H I are shown. The lower panel shows the breakdown of the contributions of CNO and Si ions. Fe III and the electron scattering acceleration, the most important contributors to wind launching, are also shown for comparison. The vertical lines represent the photospheric, sonic, and critical radii (dotted gray, dashed gray, and thin black lines, respectively). The dotted brown line marks the point where X-ray emission is switched on.

In the text

Fig. 7 Comparison of hydrodynamically consistent, pure β law profiles with classical PoWR models with hydrostatic profiles. Top panel: velocity profile as a function of the stellar radius. The gray lines indicate the pure theoretical β law profiles with β = 2.9, 3.3, and 3.9 (dot- dashed, dashed, and solid lines, respectively). The bold blue line is the hydrodynamically consistent PoWRHD model, whereas the green (continuous velocity gradient connection), purple (forced connection at v = 0.5vsound), and red (forced connection at v = 0.9vsound) profiles indicate the classical (β law + hydrostatic integration) PoWR models with β = 3.3. The light blue vertical line indicates the critical point of the PoWRHD model. Second panel: same as the top panel but dividing all the velocity profiles by the hydrodynamical solution, and hence highlighting where velocity profiles differ more. Third panel: comparison of the velocity gradient (dv/dr) profiles. Bottom panel: Spectroscopic comparison between the PoWR models (following the same color code). The observed optical spectrum (ESPRESSO; light gray) is also shown for reference.

In the text

	Fig. 8 Comparison of the hydrodynamic consistent solution (blue thick line) with fits of a double β profile (black thick line) as presented in Björklund et al. (2021). We also show for reference the best β model from Fig. 7 as the red line.
In the text

	Fig. 9 Variation in the clumping structure and output wind properties according to different characteristic onset velocities, vcl. Left panels: photosphere, i.e., where the Rosseland-mean optical depth, τ, is 2/3, as marked by the diamond symbol with its vertical line. Right panels: corresponding output spectra compared to the observed spectrum of ζ1 Sco.
In the text

Fig. 10 Upper panel: Change in the mass-loss rates, $\stackrel{·}{M}$, and terminal velocities, v∞, as a function of different values of the clumping volume filling factor, f∞. The colors trace f∞ purple indicates the homogeneous-wind model (f∞=1.0) and red is most clumped wind tested (f∞ = 0.2). The final model is represented by the black diamonds. Middle panels: changes in the wind efficiency parameter (η) and Teff as of function of f∞. Lower panels: different lines’ spectral morphology changes as a function of different clumping values. The color scheme is the same as the upper panels, with the thick black line indicating the best-fitting model spectrum.

In the text

Fig. 11 Upper panel: changes in the position of the critical (Rcrit; dark gray diamonds), sonic (Rsonίc; light gray squares), and photospheric (R2/3; colored stars connected by the blue line) radii as a function of different values of the clumping volume filling factor, f⊙ The bold black symbols indicate the best-fit model. Lower panel: same as the upper panel but showing the sonic and critical points in terms of the Rosseland-mean optical depth, τ. The blue horizontal line indicates the photosphere, where τ = 2/3.

In the text

	Fig. 12 Velocity gradient (dv/dr) vs. Rosseland-mean optical depth (τ). Each model is shifted vertically by an arbitrary amount for visualization purposes. The diamond symbols indicate the critical point, whereas the blue vertical line indicates the photosphere (τ = 2/3). The thick black curve and diamonds refer to the final model (f∞ = 0.66).
In the text

Fig. 13 Degeneracy between clumping and stellar mass. In all panels, the black lines indicate the final model, the light blue lines the model without clumping, and the dark blue lines the model without clumping with a reduced mass. Upper panel: radiative acceleration vs. Rosseland optical depth. The total radiative acceleration, ${\mathrm{\Gamma }}_{rad}^{total}$, is indicated by the thick line with dots, and the thinner lines below indicate the contribution of Fe III and Fe IV. The thin vertical lines indicate the position of the critical points. Middle panel: same as the upper panes but showing the contribution of the main ions of C, N, S, and Si. Lower panels: comparison between the models’ output spectra. Mass-loss rate units are M⊙ yr−-1.

In the text

Fig. C.1 Comparison of the hydrodynamically consistent model of ζ1 Sco (gray line) with the star’s SED and flux-calibrated UV and mid- IR spectra. The UV spectra are retrieved from FUSE and IUE (fuchsia, purple, and blue lines). The mid-IR spectrum was obtained with ISO (brown line). The cyan triangles are the flux points of the UBVI photometry. The orange diamonds are from the 2MASS JHK photometry, and the red squares are from WISE mid- and far-IR photometry. The inlet plot shows the millimeter/radio regime; the black crosses are the flux points obtained with the Very Large Array (VLA). The truncation in the synthetic spectrum is an artifact of the code output. See Appendix E.

In the text

	Fig. C.2 Resulting variation in the spectral features and wind properties when varying the T* input for PoWRHD calculations. The observed spectrum is shown in black.
In the text

	Fig. C.3 Resulting variation in the spectral features and wind properties when varying the vturb input for PoWRHD calculations. The observed spectrum is shown in black.
In the text

	Fig. D.1 IUE spectra of ζ1 Sco. The blue line is the observation, and the black line is our final PoWRHD model. The shaded area marks the interval of ±500 km s−1 around the reference wavelength.
In the text

	Fig. D.2 Variability in the spectral features in the optical of ζ1 Sco. Each column refers to a different instrument.
In the text

	Fig. D.3 Comparison of multi-epoch X-Shooter IR observations of ζ1 Sco.
In the text