© ESO 2022

1 Introduction

Massive stars strongly influence the interstellar medium by means of their winds, which carry momentum, energy, and chemical elements (see reviews from Kudritzki & Puls 2000; Puls et al. 2008). The study of their stellar winds is then an important topic in stellar astrophysics because it connects the physics of the mechanism generating the outflows of massive stars with the observed phenomena at a Galactic level, such as stellar evolution -and subsequently, stellar populations-or gas dynamics at the Galactic centre (Cuadra et al. 2008), and even at a extragalactic level, such as the so-called wind-momentum-luminosity relation (Puls et al. 1996; Kudritzki et al. 1999).

Through the pioneering work of Lucy & Solomon (1970), we know that the winds of hot stars are accelerated by the radiation field emanating from the photosphere of the star. Momentum is delivered to the wind primarily by scattering (absorption or re-emission) of photons in bound-bound (i.e. line) transitions. A quantitative description of line-driven winds was first performed by Castor et al. (1975, called the CAK theory) and was later improved by Abbott (1982), Pauldrach et al. (1986), and Friend & Abbott (1986, called modified or m-CAK theory). The approximations and assumptions underlying m-CAK and their implications for the wind dynamics have been extensively discussed (Schaerer & Schmutz 1994; Puls et al. 2008; Krtička & Kubát 2010); the most important approximation is the Sobolev approximation (Sobolev 1960; Castor 1970), which allows a description of the line-acceleration in terms of local wind quantities and the photospheric flux (Castor et al. 1975; Lamers & Cassinelli 1999).

The beauty of m-CAK theory is that it provides key insights into the relation between the line-force parameters (k,α,δ) and fundamental parameters describing the atomic transitions responsible for driving the flow, such as their gf-values and their line-statistics (see Puls et al. 2000). Moreover, these line-force parameters provide a solution for the wind hydrodynamics, which is beyond the so-called β-law. Another approach for studying stellar winds is to use the Sobolev approximation in conjunction with Monte-Carlo simulations. One of the most well-known studies is the Monte-Carlo simulations performed by Vink et al. (2000, 2001), whose theoretical mass-loss rates have been widely used to calculate evolutionary tracks (Georgy et al. 2012; Ekström et al. 2012). Unfortunately, improved observational studies have shown discrepancies between these Vink rates and observations (e.g., Bouret et al. 2005, 2012).

In parallel, a renewed effort has been made in past years to improve our understanding of line-driven winds beyond the m-CAK. Sander et al. (2017) calculated consistent line-acceleration and hydrodynamics for the particular case of ζ-Puppis using the PoWR code (Gräfener et al. 2002; Hamann & Gräfener 2003), whereas Krtička & Kubát (2017, hereafter KK17) provided new confident theoretical values for wind parameters of O dwarfs, giants, and supergiants by solving the radiative transfer equation and hydrodynamics simultaneously using the Sobolev approximation for transition rates and co-moving frame (CMF) for the calculation of the radiative force. We also mention the recent studies of Sundqvist et al. (2019) and Björklund et al. (2021), who used their full non-local (NLTE) CMF radiative transfer solutions from the code FASTWIND (Santolaya-Rey et al. 1997; Puls et al. 2005; Sundqvist & Puls 2018) to calculate the radiative acceleration responsible for wind driving. In the same line, we highlight the study by Gormaz-Matamala et al. (2021), who also performed a self-consistent calculation of hydrodynamics and radiative acceleration under a full NLTE treatment using CMF-GEN (Hillier 1990a,b; Hillier & Miller 1998). However, all these works require a considerable computational effort, which at the moment prevents the creation of a large grid of self-consistent solutions for a wide range of temperatures and masses.

Based on the m-CAK theory, Gormaz-Matamala et al. (2019, hereafter Paper I) developed a prescription to derive self-consistent solutions given different sets of stellar parameters for hot massive stars, following the quasi-NLTE treatment for atomic populations of Mazzali & Lucy (1993), the radiation field calculated from TLUSTY (Hubeny & Lanz 1995), and the hydro-dynamic solutions performed by HYDWIND (Curé 2004; Curé & Rial 2007). Likewise, Araya et al. (2017) self-consistent velocity profiles calculated from this m-CAK prescription can be introduced in FASTWIND to calculate synthetic spectra without any β-law. Paper I demonstrated that this prescription provides useful theoretical wind parameters for different sets of temperature, surface gravities, abundances, and even the rotational velocity, whose derived synthetic spectra quickly approach a fair fit. Moreover, the line-force parameters have recently been used to cover a broader range of temperatures, although an LTE scenario was assumed and an optical depth was considered that are inconsistent with the hydrodynamics (Lattimer & Cranmer 2021). In this work, we proceed in the usage of the m-CAK prescription performed by Paper I and it to the spectral fitting of a sample of late-type O stars observed by the HERMES spectrograph. From this spectral analysis, we obtain new stellar and wind parameters for the stars of the sample (HD 192639, 9 Sge, HD 57682, HD 218915, HD 195592, and HD 210809), which are compared with previous studies. We also introduce small improvements in the calculation of the line-force parameters, such as a stratified temperature profile and a more robust handling of the errors for our solutions.

This paper is organised as follows. Details of these improvements and a brief description of their implications are given in Sect. 2. Details of the HERMES spectrograph and the spectra of the sample are provided in Sect. 3. A description of the steps required to obtain the synthetic spectra from the self-consistent solutions is given in Sect. 4. The results, with the new stellar and wind parameters measured for our sample, are presented in Sect. 5. Finally, Sects. 6 and 7 offer the discussion and conclusions of our work, respectively.

2 Self-consistent solutions under m-CAK prescription

The study performed in Paper I provided an exhaustive analysis of the line-force parameters k, α, and δ from the CAK theory (Castor et al. 1975; Abbott 1982), calculating them self-consistently with the hydrodynamics derived for the wind by the code HYDWIND (Curé 2004). We call the velocity υ(r) and density ρ(r) profiles that describe the wind wind hydrodynamics. These two profiles are closely coupled by the equation of momentum, $$\upsilon \frac{\text{d}\upsilon }{\text{d}r}=-\frac{1}{\rho }\frac{\text{d}p}{\text{d}r}-\frac{G{M}_{\text{eff}}}{r^2}+g_{\text{es}}k{t}^{-\alpha }{\left(\frac{N_{\text{e}}}{W}\right)}^{\delta },$$(1)

where the term ∝dp/dr is the dependence on the pressure gradient, and the term ∝M_eff/r² is the dependence on the gravitational force. The last term corresponds to the line-acceleration, where the line-force parameters (k, α, δ) describe the force multiplier $$\mathcal{M}(t)=k{t}^{-\alpha }{\left(\frac{N_{\text{e}}}{W}\right)}^{\delta }=\frac{g_{\text{line}}}{g_{\text{es}}},$$(2)

which is given this name because it multiplies the acceleration due to electron scattering, g_es. We note that ℳ(t) is presented as a function of t, the CAK optical depth $$t={\sigma }_{\text{es}}{\upsilon }_{\text{th}}\rho (r){\left(\frac{\text{d}\upsilon }{\text{d}r}\right)}^{-1},$$(3)

with σ_es being the electron scattering opacity and υ_th the mean thermal velocity of the protons of the wind. This t term is defined by the CAK theory (Castor et al. 1975) and it only depends on the wind structure, differently from the usual optical depth τ, which depends on the structure and composition of the wind. The force multiplier also depends on the ionisation density N_e/W, which corresponds to the electron density divided by the geometrical dilution factor W.

Hence, k, α, and δ are calculated by fitting ℳ(t) over a 2D plane formed by the CAK optical depth t and the ionisation density Ne/W, where t is constrained in the range in which the Sobolev approximation can be assumed (i.e. the inner limit is the sonic point of the wind, and the outer limit is the infinite). As previously noted in Paper I, the optical depth t decreases at infinitum, but never reaches zero. The reason is that at larger distances, both density profile and velocity gradient in Eq. (3) become proportional to ~r⁻², so that both factors cancel each other out. The meaning and effects of each of the line-force parameters over the line-acceleration and subsequently over the wind parameters are summarised in Gormaz-Matamala et al. (2019, Sect. 2). In addition to Eq. (1), the velocity and density profiles are also connected to each other by means of the mass-conservation equation, where the term for the mass-loss rate, M, is introduced, $$\rho (r)=\frac{\dot{M}}{4\pi r^2\upsilon (r)}.$$(4)

As a consequence, new theoretical values for the wind parameters (mass-loss rate, terminal velocity) can be determined for any specific set of stellar parameters (effective temperature, surface gravity, stellar radius, abundances). In the particular case of mass-loss rates, self-consistent values have proven to agree with those determined by observations when homogeneous wind is assumed (see Fig. 13 in Gormaz-Matamala et al. 2019); whereas for the case of clumped winds, the clumping factor is a free parameter in the spectral fitting. This is an important result because these spectral fittings were achieved without using the power law that is commonly assumed to describe the velocity profile of stellar winds, namely β-law $$\upsilon (r)={\upsilon }_{\infty }{\left(1-\frac{R_0}{r}\right)}^{\beta }.$$(5)

This description, even when adjusted to the observations, is not derived from hydrodynamic calculations and then is justified a posteriori (Kudritzki & Puls 2000). Therefore, self-consistent velocity fields cannot be described by any β exponent (see Fig. 1), and thus they need to be calculated from the m-CAK equation of motion (Eq. (1)) following the recipe introduced by Curé (2004). This departure from β-law agrees with that observed for the self-consistent solutions calculated from the Lambert-procedure in Gormaz-Matamala et al. (2021), where the new velocity profile could not be recovered by any β exponent. However, it is necessary to remark that the calculated velocity field differs from the line-force parameters, compared with hydro solutions in the CMF: if we compare our Fig. 1 with Fig. 9 from Gormaz-Matamala et al. (2021), we observe a steeper profile at the base of the wind for the case of v(r) calculated from m-CAK prescription. This is because the CMF approach accounts for non-Sobolev effects such as source-function dip, back-scattering, or line overlaps, and line-acceleration calculated in CMF additionally does not present an explicit dependence on velocity gradient or velocity (see discussion from Sect. 5.3 of Sundqvist et al. 2019), and therefore it is very different from the CAK framework. Nevertheless, in spite of these discrepancies between the velocity profile for the subsonic region of the wind, they are not translated into the synthetic spectra. We performed a posteriori self-consistent solutions with a subsonic region based on the CMF approach (following the coupling strategy introduced in Gormaz-Matamala et al. 2021, see their Sect. 2.1.1) and eventually found that the resulting differences in the synthetic spectra are just marginal. Hence, we continued without modifying the subsonic region for our self-consistent velocity profiles.

The spectral fitting for the self-consistent solutions from the m-CAK prescription was performed with the radiative transfer code FASTWIND (Santolaya-Rey et al. 1997; Repolust et al. 2004; Puls et al. 2005), which has the option of using a velocity profile generated from a formal solution of the equation of motion as input instead of a classical β-law (Araya et al. 2017). Thus, it is possible to create synthetic spectra in which the wind parameters are not longer free, but depend on the stellar parameters. This reduces the number of free parameters. Based on this statement, we implement a method for fitting spectral observations in which the wind parameters are subordinate to the fit of the stellar parameters that generate them.

The self-consistent m-CAK solutions provided by Paper I were shown to be reliable for a wide range of O stars. Studies aiming for a self-consistent prescription under a full NLTE treatment, such as Gormaz-Matamala et al. (2021), have provided only slightly different results. However, we considered improvements of our m-CAK prescription necessary, especially to rectify the high terminal velocities provided by it. For this reason, we introduce in the next subsection a stratification for the radial temperature of the wind, which replaces the assumption of an isothermal wind for the calculation of the line-force parameters.

Fig. 1 Velocity profile generated from a self-consistent wind solution using the m-CAK prescription, compared with some v(r) following β-law. The bottom panel corresponds to a zoom into the base of the wind.

2.1 Stratification of temperature

A complete inclusion of a temperature structure requires a sophisticated calculation of the radiative equilibrium and the electron balance through the wind, which is commonly made using radiative transfer codes such as CMFGEN, PoWR, or FASTWIND. Because the solution of radiative equilibrium is beyond the scope of our procedure, we introduce a stratification for the temperature field over the wind based on the simplified prescription used previously by Sundqvist et al. (2019), and based on the Lucy (1971) estimation for a temperature structure in a spherically extended envelope, $$T(r)=T_{\text{eff}}{\left[W(r)+\frac{3}{4}{\tau }_{\text{F}}\right]}^{1/4},$$(6)

with W(r) being the dilution factor $$W(r)=\frac{1}{2}\left(1-\sqrt{1-\sqrt{R_{*}^2/r^2}}\right),$$(7)

τ_F being the flux-weighted optical depth $${\tau }_{\text{F}}(r){\displaystyle \int {\kappa }_{\text{F}}(r^{\prime })\rho (r^{\prime }){\left(\frac{R_{*}}{r^{\prime }}\right)}^2\text{d}{r}^{\prime },}$$(8)

and κ_F being the flux-weighted mean opacity $${\kappa }_{\text{F}}(r)F(r)={\displaystyle {\int }_0^{\infty }{\kappa }_v{F}_v\text{d}v}.$$(9)

This means that the mean opacity, even when it typically is calculated by radiative transfer, is directly connected to the radiative acceleration of the stellar wind (see Eq. (3) in Sundqvist et al. 2019), $$g_{\text{red}}(r)=\frac{{\kappa }_{\text{F}}(r)L_{*}}{4\pi c{r}^2}=\frac{{\kappa }_{\text{F}}(r)F(r)}{c},$$(10)

Then, including the relation between radiative and line-acceleration (Eqs. (3) and (6) from Gormaz-Matamala et al. 2021), it is possible to rewrite the flux-weighted optical depth as $${\tau }_{\text{F}}(r)=\frac{4\pi c{R}_{*}^2}{L_{*}}{\displaystyle \int \left[g_{\text{line}}(r^{\prime })+g_{\text{grav}}(r^{\prime }){\Gamma }_{\text{e}}\right]\rho (r^{\prime })\text{d}{r}^{\prime },}$$(11)

with g_grav(r) = GM_*/r² and Γ_e being the Eddington factor (Gormaz-Matamala et al. 2021, Eq. (4)). For details, see the step-by-step calculation in Appendix A.

Because the entities inside the integral are calculated by our m-CAK prescription, it is possible to compute the flux-weighted mean optical depth simultaneously with the iterative loop calculations, and thus the final stratificated temperature is self-consistent with the final line-acceleration. Like Sundqvist et al. (2019) and Puls et al. (2005), we imposed a ground temperature of T(r) = 0.4 T_eff to avoid excessively low temperatures at larger distances.

2.2 Temperature structure and line-force parameters

The direct effects of the change in temperature over the wind are the alteration of the excitation and ionisation stages for the atomic populations, being the lower levels more present in the more external parts of the wind with respect to the case of a fixed temperature for the wind. Since lower stages of ionisation tend to have more lines (Abbott 1982; Puls et al. 2000), a temperature structure should generate a larger contribution to the line force in the outer region of the wind, which can be translated into higher values for the line-force parameter α. However, this trend is observed only for low values of T_eff: variation in α passes from diminishing its value (in comparison with the fixed-temperature case) for hotter effective temperatures to slightly enhance it for cooler T_eff. At the same time, the variation in δ (directly linked with the ionisation of the wind) seems to be just marginal, with the exception of the 32 kK star.

Another consequence of the stratification of the temperature is the restructuring of the optical depth t (Eq. (3)), which is implicitly dependent on T(r), in particular, v_th, which is the mean thermal velocity of the protons of the wind, $${\upsilon }_{\text{th}}(r)=\sqrt{\frac{2k_{\text{B}}T(r)}{m_{\text{H}}}},$$(12)

with k_B being the Boltzmann constant. Because the temperature through the wind is no longer constant, this thermal velocity decreases as radial distance increases, and subsequently, the optical depth. As a result, we observe from Fig. 2 that the range in which the force multiplier is calculated appears to be shifted to lower values for t. Subsequently, this displacement produces a direct decrement of the line-force parameter k, and therefore also of α and δ due to the iterative loop searching for self-consistency.

The wind parameters calculated from our four analysed standard stars are tabulated in Table 1. When the temperature structure is included, the mass-loss rate is clearly most strongly affected: from a great decrease at high temperatures (45 kK) to a marginal gain at cooler temperatures (34 kK), and another high declination for the case T_eff = 32 kK (which can altogether be explained by the reduction in the line-force parameter δ). For modifications of the terminal velocities are directly linked with the correction of the values for α, and they are again due to the reduction of δ for the coolest case.

This procedure of calculating a stratified temperature was made under the assumption that we can recover opacities and emissivities from our m-CAK line-acceleration. In other words, our implemented T(r) is not calculated by radiative transfer and therefore will differ from a fully calculated temperature structure. Because our synthetic spectra are determined by FASTWΓND, we can compare the stratificated temperature that was used during the iterative procedure calculating our self-consistent parameters and the temperature obtained with opacities and emissivities calculated by radiative transfer. These comparisons for two stars of our sample are shown in Fig. 3, where we observe that differences in T(r) are ~3 kK (except for the external part of HD 210809). It might be argued that the similarity in the results is forced because the two prescriptions use the same boundary conditions (where T(r) = T_eff for the photosphere and T(r) = 0.4 T_eff for the imposed external limit), but even in this case, the replacement of our m-CAK temperature structure by the fully calculated by FASTWIND produces only marginal variations in the self-consistent parameters with a higher computational effort. Therefore, the prescription for the stratification of temperature introduced in this work, which is considered as a general approximation, quickly leads to more reliable results than the case T(r) = T_eff for our self-consistent wind parameters.

2.3 Variance of line-force parameters

The range of validity for the self-consistent solutions from Paper I was established for effective temperatures from 32 to 45 kK and for surface gravities with log g ≥ 3.4 because within these limits, the uncertainties over the final wind parameters resulting from the fitting of the line-force parameters were smaller than the error bars derived from the observational uncertainties of the stellar parameters. This was made on the basis of a qualitative analysis of the behaviour of the force multipliers ℳ(t). We aim to complement this outline by including a quantitative measurement of the fit for the line-force parameters by means of calculating their respective standard deviations. Hereafter, error bars associated with the wind parameter mass-loss rate and terminal velocity come from the variation in (k, α, δ) instead of from the variation in stellar parameters.

It is important to recall that these analyses over the error bars come from the assumption that line-force parameters are constant throughout the wind instead of a function of depth, as has been proposed by authors such as Schaerer & Schmutz (1994) and Kudritzki (2002). In other words, we can consider k, α, and δ as constants only if their uncertainties lies below the uncertainties found for the stellar parameters. This condition is satisfied for the range of temperatures we studied, where log ℳ(t) presents an almost linear shape as a function of log t (see Fig. 2).

From Paper I, the uncertainties on the stellar parameters generated errors of ~7−15% for terminal velocities and of ~20% for mass-loss rates (~0.1 in logarithmical scale). Hence, we could consider these values the threshold of the maximum permisible error on our models, with special emphasis on M. For this reason as well, the values for the mass-loss rate and their uncertainties are always given in logarithmical scale here.

The line-force parameters introduced in Table 1 are given with errors, which come from the standard deviation for each parameter when the fit is performed. Likewise, the terminal velocities and the logarithm of the mass-loss rates contain error bars that represent the new associated uncertainties for our self-consistent wind parameters. These uncertainties do not come only from the standard deviation of (k, α, δ), but also from the discrepancies produced by numerical issues such as the selected number of depth points when the equation of momentum (hereafter e.o.m., Eq. (1)) is calculated by HYDWIND (for details of the numerical calculation of the e.o.m., see Curé 2004). This is calculated by summing them quadratically, $${\Delta }_{\text{wind}}=\sqrt{{\Delta }_{(k,\alpha ,\delta )}^2+{\Delta }_{N_{\text{points}}}^2},$$(13)

with ∆(k,α,δ) and ${\text{Δ}}_{N_{\text{points}}}$ being the error intervals either for v_∞ or log $\dot{M}$ generated by the variation in line-force parameters and the number of points, respectively.

For our selected four standard stars, self-consistent wind parameters were calculated for a set of different numbers of depth points from 200 to 900 (in intervals of 100), plus the maximum number of points allowed by the code (999). An increase in the number of depth points results in a fine grid for the wind structure, but implies higher time consumption. By normalising these values by the terminal velocity and mass-loss rate under the maximum N (v_∞,999 and M₉₉₉), the deviation generated for the use of different number of points is illustrated better, as is shown in Fig. 4. It is observed that only marginal discrepancies are created for mass-loss rates of ~0.03 in logarithmical scale (below the threshold of 0.1 introduced above), whereas a remarkable deviation for the terminal velocity appears for the star with T_eff = 32 kK and log g = 3.4. Therefore, uncertainties over v_∞ increase considerably for lower values of the effective temperature and surface gravity, which reach ~12−13%, as is shown in Table 1, near to the ~15% previously measured from the variation in stellar parameters. However, we emphasise that these are partial trends, and more models are needed to generate a definitive conclusion.

Including these error intervals allows us a more quantitative analysis of the validity range of our self-consistent m-CAK prescription. This gives us the chance to explore and evaluate the existence of wind solutions beyond the rigid thresholds derived in Paper I. Throughout this work, we scrutinise the presence or absence of solutions for effective temperatures cooler than 32 kK and surface gravities below 3.4.

Fig. 2 Comparison of self-consistent fitting of the force multiplier M(t) between the implementation of the temperature structure (blue) and the fixed temperature (red). Calculated line-force parameters (k, α, δ) are displayed in the same colour code. The shaded area represents the t optical depth range where the force multiplier is defined within, from the sonic point outwards.

Fig. 3 Temperature stratification for HD 57682 (leftpanel, Teff = 35 kK) and HD 210809 (rightpanel, Teff = 31.5 kK), either obtained from Eq. (6) (red) or calculated by FASTWIND (blue).

Fig. 4 Comparison of the obtained self-consistent wind parameters v∞ and M, employing different numbers of depth points N to numerically solve the equation of motion. Values are expressed as departures from the self-consistent terminal velocity and mass-loss rate when the maximum number of points is allowed for HYDWIND (999).

2.4 Range of validity

Before we implement the synthetic spectra to fit observed spectra, a quick evaluation of our self-consistent wind parameters with respect to previous studies is necessary. Analogously to Paper I, we compared our resulting self-consistent mass-loss rates, M_sc with the theoretical values provided by the Monte Calo procedure of Vink et al. (2000, 2001), and for this case, we specifically include a comparison with the self-consistent wind solutions for O-type stars calculated by KK17. We selected these two studies because they present a surprisingly well-correlated theoretical mass-loss rate, with the values from KK17 being about three times lower (a difference of ~0.5 in logarithmical scale) than those given by the so-called Vink formulae.

For this purpose, we took the stellar parameters from the sample introduced by KK17 and calculated our own self-consistent wind solutions. The results are displayed in Table 2, and we compare the different obtained mass-loss rates in Fig. 5. Two stars with effective temperatures of 30 kK, 300-V and 300-I (see Table 1 from KK17), were not included in this comparison. The main-sequence star (300-V) was excluded because the values of stellar radius (6.6 R_⊙) and mass (12.9 M_⊙) are too to solve the e.o.m. by HYDWIND, and the supergiant (300-I) was excluded because the low temperature and surface gravity led to an unstable model without a unique self-consistent solution (for details, see Sect. 6.1). For the remaining stars, it is interesting to note that high effective temperatures generate values for M_sc that agree more closely with those provided by KK17, whereas (with the remarkable exception of 325-V, the outlier of our sample) mass-loss rates from cooler temperatures show a better correspondence with Vink’s formulae. This unexpected trend is curious because first, our theoretical self-consistent mass-loss rates differ from the good correlation observed between M_Vink and M_KK17, and second, it would be possible to establish a smooth transition for theoretical M at some point around T_eff 32.5 kK and log g 3.5−3.3 for future evolutionary tracks of massive stars if we were to change the mass-loss prescription from self-consistent values to Vink’s formulae (Gormaz-Matamala et al., in prep.).

Fig. 5 Theoretical mass-loss rates tabulated in Table 2 for the values predicted by Vink’s formula (orange symbols), by KK17 (green symbols), and by this study (blue symbols). Main-sequence stars of the sample are represented by starss, giants by triangles, and supergiants by circles.

3 Observational data

The selected sample is composed of a set of O supergiants (HD 192639, 9 Sge, HD 218915, HD 210809 and HD 195592), which are situated at the lower border of validity for effective temperature and surface gravity as introduced in Sect. 2.4. Additionally, we included a star known by its magnetic field (HD 57682) in order to confirm the quality of a spectral fitting for wind structures departing from spherically symmetry. For the spectral data, the sample was chosen based on the high quality of the signal-to-noise ratio (S/N) and for stars with compelling evidence that they are single stars (i.e. excluding any possible binarity).

3.1 HERMES spectra

The spectra were observed with the HERMES echelle spectro-graph mounted on the Mercator 1.2 m telescope at the Roche de los Muchachos Observatory, La Palma, Spain. HERMES (High Efficiency and Resolution Mercator Echelle Spectrograph) is a high-efficiency prism-cross-dispersed fibre-fed bench-mounted spectrograph that observes the complete wavelength range 3800–9000 Å in a single exposure at a spectral resolution of R = 85000 (Raskin et al. 2011).

Each star was observed using the high-resolution fibre mode in one to five exposures, depending on the V magnitude, in succession throughout the same night (see Table 3). Additional flat-field exposures were observed during the night to avoid introducing systematic noise in the flux calibration process. The spectra were calibrated with the latest version of the HERMES pipeline (release v6.0) developed at the Royal Observatory of Belgium and the University of Leuven, in collaboration with the HERMES Consortium. The typical calibration steps were performed, including spectral order tracing and extraction, average flat-fielding, Th-Ar lamp wavelength calibration, and cosmic ray removal using cross-order profiling.

The calibrated exposures were co-added to produce a single stacked 1D spectrum per star with an S/N of ~ 150–250. The S/N was estimated at four effective wavelengths (U, B, V, and R bands). The typical uncertainty of the spectrograph wavelength calibration is below ± 0.02 Å, or the HERMES wavelength scale has a high accuracy.

3.2 Continuum normalisation

The raw exposures were calibrated with the removal of cosmic rays for the final sum spectra. Dedicated software tools were developed for BRASS (Lobel et al. 2019) to remove cosmic ray flux spikes in HERMES spectra. Moreover, the spectral response function of the instrument shows a strong wavelength dependence and some time variability. However, spectra observed with regular and smooth continuum flux shapes are required to perform automated continuum normalisation calculations. Consequently, all signatures of the instrumental response were removed from the observed spectra. Advanced HERMES calibration pipeline algorithms were previously developed in BRASS to determine the response correction curve at every point in time.

The method for modelling the instrumental response curve is important to perform high-quality continuum flux normalisations, including in the small wavelength regions around selected spectral lines. We automatically normalised the observed spectra to the stellar continuum flux level using a special template normalisation procedure. It searches for variable wavelength points over sufficiently continuous flux regions close to the continuum level in theoretical (template) spectra to fold the HERMES spectra to these (continuum anchor) points. This ensures that local flux normalisation effects are minimised, providing reliably normalised spectra around the H and He lines we select for our detailed line profile modelling.

4 Method

4.1 Selection of FASTWIND spectral lines

The group of 15 lines calculated by FASTWIND is presented in Table 4.

As a first stage, we decided to use only H and He lines for the spectral fitting in the optical and infrared wavelength range for practical purposes, whereas our observational data consist of optical spectra alone.

Fig. 6 Scheme of the self-consistent m-CAK procedure combined with FASTWIND, showing the stage at which the different initial parameters (in rectangles) are set. Ellipses represent the processes or codes to be run, and the diamonds represent the outputs.

4.2 Self-consistent modification of parameters

The scheme showing the steps required to achieve a self-consistent synthetic spectrum is shown in Fig. 6. This procedure is the same as we used to obtain the spectra for Paper I, but now we go one step further: the spectral fit is improved by the fine-tuning of the free parameters (represented as rectangles in the figure) at the different stages of the process.

The upper half of Fig. 6 corresponds to the iterative process that searches for the self-consistent solution under our m-CAK prescription, therefore its only input parameters are the stellar parameters (temperature, gravity, and abundances). When this solution is achieved, that is, when line-force parameters converge and theoretical values for the mass-loss rate, terminal velocity, and velocity profile are obtained, a synthetic spectrum is computed by FASTWIND using the new velocity profile v(r) provided by the solution of (k, α, δ) as an input. As we previously stated in Paper I, a clumping description is introduced at this point, together with the microturbulence υturb. By clumping description we mean not only f_cl, but also the starting and ending zones (rin and rout) in which the inhomogeneities are located. This means that we use a clumping law in which f_cl is constant over a region of the wind that is delimited by rin and rout (see Puls et al. 2006, for details). Finally, convolution for the resulting synthetic spectrum requires inserting a value for rotational velocity υrot and for macroturbulence v_macro as input. We therefore observe that input parameters can be classified as first-order (stellar parameters, whose modification implies a full restructuring of the self-consistent wind solution), second-order (clumping and turbulence, to recalculate the FASTWIND atmospheric model but without modifying the self-consistent wind solution), and third-order parameters (rotation and macroturbulence, to reshape line profiles).

Although one of the advantages of developing synthetic spectra derived from self-consistent wind prescription is that the number of free parameters is reduced, the number of remaining parameters that need to be fit freely is still large enough to require a deeper evaluation of their effects upon the spectral lines. After many empirical tests and comparisons in which only one parameter was let free to vary while the others were held fixed, combined with previous analyses found in the literature, we observed some general trends that helped us to determine the best model fit. We describe them below.

• The effective temperature T_eff is fitted by the ratio between the He II and the He I, given the expected change in ionisa-tion. Because it is a first-order parameter, modifications of the effective temperatures imply an alteration of the self-consistent mass-loss rate (which are directly dependent), so that an increase or decrease in T_eff will also affect the intensity of the hydrogen lines, especially Hα. Moreover, the effective temperature also provides information about the stellar radius, which is constrained with the luminosity (i.e. observed magnitude) by means of the Stefan-Boltzmann law. Therefore, modifications in T_eff require modifications in R_*/R_⊙ for the models to maintain coherence with the V magnitudes tabulated in Table 3.

• The surface gravity log g is evaluated by the absorption intensity observed for Hß, Hγ, Hδ, and He, together with the absorption profiles of the He I lines. Because it also is a first-order parameter, gravity presents an inverse proportionality with the self-consistent mass-loss rate. As a consequence, its increase implies a decrease in the potential emission component of Hα.

• The turbulent velocity vturb, also known as microturbu-lence velocity, is fitted to the H I lines. In concordance with Villamariz & Herrero (2000), we observe that the effect of vturb is only perceptible for hotter values of effective temperature (in our case, T_eff ≳ 33 kK), primarily over the lines He I λ4387 and He I λ4922.

• The clumping factor f_cl is fitted to the hydrogen lines, in particular, it fits the wings of Hα. As was pointed out by Repolust et al. (2004), clumping effects over the wings of Hα are more prominent when they are in emission than in absorption. In addition, we observe an effect over the core and the red wing over the line profiles of Hβ, Hγ, Hδ, and Hϵ as a function of the wavelength. Here, we follow the same notation as we used in Paper I that was taken from Puls et al. (2006), where the clumping factor in FASTWIND is denoted as f_cl and takes values ≥ 1 (where f_cl = 1 is the value for the smooth limit). This f_cl is then the inverse of the volume filling factor, denoted in CMFGEN as f_∞ (Bouret et al. 2012; Gormaz-Matamala et al. 2021). This scale f_cl = 1/f_∞ relies on the assumption that the inter-clump medium is empty (Sundqvist & Puls 2018).

• The macroturbulent velocity vmacro, which corresponds to the non-rotational contributions to the line broadening (see Simón-Díaz & Herrero 2014), is fitted to the width of He I lines, which have been empirically shown to be more sensitive to this broadening.

• The rotational velocity v_rot is fitted to the slope in the shape of the wings of He II λ4200 and He II λ4541. For this purpose, the synthetic spectrum needs to be convolved after the end of the FASTWIND process. Moreover, we decided to make this convolution coherent with the effects of the Ω parameter over the self-consistent wind hydrodynamics,

$$\Omega =\frac{{\upsilon }_{\text{rot}}}{{\upsilon }_{\text{crit}}},$$(14)

where v_crit is the critical rotational speed for a spherical star (see e.g. Araya et al. 2017, for details). The rotational velocity is a third-order parameter, but Ω is a first-order parameter, which means that modifications in the convolution may require a complete restructuring of the self-consistent solution. Because we focus on the theoretical performance of synthetic spectra and because the focus is not on rapidly rotating stars, we prefer to denote rotational velocity simply as vrot instead of the classical v sin i. The convolution to broaden the spectral lines by rotational effects and the inclusion of the Ω factor over the hydrodynamic solution are both assumed to happen in the line of sight to the observer.

As a final remark, we mention that the radial velocity, v_rad, is obtained from the SIMBAD database¹.

5 Spectral results

In this section we present the results of the spectral analysis for the six stars of the sample, HD 192639, 9 Sge, HD 57682, HD 218915, HD 195592, and HD 210809. For each of them, we present the best spectral fit (based on the trends introduced in Sect. 4.2) and a tabulation with the final stellar, wind, and line-force parameters found for these stars. An analysis and comparison of some of these new parameters with the literature, especially for the case of the mass-loss rate, are also considered.

The spectra were fit with a by-eye inspection that searched for the model that globally satisfied all the criteria described above better. We selected this option instead of a χ² reduction or an algorithm procedure (e.g. Mokiem et al. 2005) because an inspection by eye allows a more qualitative discussion of the fits. Of the 15 lines introduced in Table 4, the focus is mostly given to the fit of the wings for the set of hydrogen lines and to the emission peak of Hα. For helium, good-quality fits are expected for all the lines, except for He II λ4686 and He I λ6678, whose appropriate synthetic reproduction has been a challenge for the codes FASTWIND and CMFGEN (Puls et al. 2005; Holgado et al. 2018). The short execution time of the combined m-CAK procedure with FASTWIND (Fig. 6), which takes about one hour, allows developing multiple self-consistent wind solutions and then calculating the best model fit is quickly achieved with ~30 model runs.

Same as for Tables 1 and 2, the self-consistent wind parameters introduced in this section are presented with the error bars derived from the uncertainties of the line-force fitting and those derived from differences in the numbers of depth points used for the wind structure, as outlined in Sect. 2.3. Variance in the wind parameters due to uncertainties in the fitted stellar parameters (T_eff, log g, He/H, etc.), which is a priori assumed to be about ±0.125 for log M, is not included in these error bars.

5.1 HD 192639

The first star we analysed was the blue supergiant O7.5 I, HD 192639, observed with HERMES in 2013. This star was previously studied by Bouret et al. (2012), who fitted a CMFGEN synthetic spectra with parameters T_eff = 33.5 kK, log g = 3.42, log M = −5.68, andβ = 1.3. We used these parameter as a starting point to perform our calculations instead of those derived from the calibrations by Martins et al. (2005) for a standard O 7.51 star (T_eff = 32 kK and log g = 3.36), which differ moderately, but still agree closely. The stellar parameters we derived are tabulated in Table 5 and also lie close to theirprevious values. The synthetic spectra are shown in Fig. 7.

The initial rotational velocity for HD 192639 is slightly increased from v_rot = 90 km s⁻¹ given by Bouret et al. (2012) to 100 km s⁻¹. This rotational velocity corresponds to an angular velocity of Ω = 0.26, which has a real impact on the generated velocity profile and on the line profiles, making them broader. This additional rotational broadening, which comes not only from the final convolution over the spectra, but also from the self-consistent solution, can partly replace the macrotur-bulent velocity, which is set to v_macro = 30kms⁻¹ instead of the 43kms⁻¹ given by Bouret et al. (2012). Another important upgrade: the self-consistent solution provides a higher effective temperature than the initial value, whereas the surface gravity is ~10% lower. The final difference with the CMFGEN fit is our reduction of the He to H ratio from 0.15 to 0.1, which is necessary to reduce the intensity of all helium lines.

The spectral fit we achieved shows a good match and is particularly good for the emission component of Hα, and with the remarkable exception of He II λ 4686 and its inherent variability, the match is of the same order of accuracy as Bouret et al. (2012, see their Fig. A.16). Despite the mentioned variability, this fit gives us confidence for the self-consistent value we found for the mass-loss rate, log $\dot{M}=-5.783$. This is just ~38% higher (+0.14 dex in logarithmical scale) than the value found by Bouret et al. (2012), but we used a clumping factor that is 3.2 times lower.

Fig. 7 Best fit for HD 192639 with the self-consistent solution tabulated in Table 5.

5.2 9 Sge

9 Sge (HD 188001) is a runaway star (Underhill & Matthews 1995) with a spectral type of O 7.5 lab (Sota et al. 2011). Do not confuse 9 Sge (from the Sagitta constellation) with 9 Sgr (from the Sagittarius constellation). This latter one is also an O starm but belongs to the main sequence spectral type O 3.5 V (Sota et al. 2014). 9 Sge exhibits a periodic variation in its radial velocity of P = 97.6 days, but it is still debated whether this period is produced by the presence of a companion (Maíz Apellániz et al. 2019). Based on its spectral classification, the initial stellar parameters taken from the catalogue of Martins et al. (2005) are: T_eff = 34, log g = 3.36 and R_*. = 20.8 R_⊙, which are also the parameters provided by the VizieR Online Data Catalog² (Nebot Gómez-Morán & Oskinova 2018). Alternatively, Martins et al. (2015a) found T_eff = 33, log g = 3.35. The fitted parameters and the spectrum are presented in Table 6 and Fig. 8 respectively.

To achieve a closer agreement for the helium lines, it was necessary to increase the He to H ratio from the solar standard from Asplund et al. (2009, He/H = 0.085) up to He/H = 0.125. Moreover, wide shape of the line profiles led us to require an increase in both macro-turbulence (25kms⁻¹) and rotational velocity (90 km s⁻¹), together with perform wind solutions with a lower gravity. Same as for the case of HD 192639, Hα is fitted by its wings and by the emission component, whereas the absorption is only partially reproduced. The self-consistent mass-loss rate obtained is 2.33 × 10⁻⁶ M_⊙ yr⁻¹, roughly ~l.8 times higher than valued provide by Nebot Gómez-Morán & Oskinova (2018).

Fig. 8 Best fit for 9 Sge with the self-consistent solution tabulated in Table 6.

5.3 HD 57682

We modelled the wind of the 09 star HD 57682, which is the only non O-supergiant star of our HERMES sample. It was originally classified as luminosity-class IV by Walborn (1972), but reclassified as main sequence (luminosity-class V) by Grunhut et al. (2012). The study of Grunhut et al. (2009) found line profile variability (LPV) due to a magnetic field B ~ 1600 G. Despite this high value for B, these studies have successfully achieved accurate spectral fittings using the spherically symmetric code CMFGEN. Based on this study, we started with the following stellar parameters: Γ_eff = 35 kK, log g = 4.0, and = 7R_⊙. In this particular case, the self-consistent solution obtained from this set of stellar parameters quickly leads to a good fit after little tuning of the effective temperature and surface gravity. The spectrum is presented in Fig. 9, and the obtained parameters are shown in Table 7.

Grunhut et al. (2012) found that the most prominent LPV for HD 57682 is observed for Hα, where a strong emission peak in the core of the line arises. Moreover, the circumstellar material surrounding the star also contaminates the core of the other hydrogen lines such as Hγ (Grunhut et al. 2009, see their Fig. 2). For the case of helium lines, we derive a slight increase in the effective temperature of (500 K) and an also slight decrease in surface gravity (0.1 dex), which helps to improve the fit for He I lines in particular (He I λ24471, probably because of its variable behaviour).

For the mass-loss rate, Grunhut et al. (2009) determined a value of log $\dot{M}=-8.85$ from the C IV lines λ21548 and 21551, but they stated that the mass-loss rate derived from the emission component of Hα should be ~ 10⁻⁷ solar masses per year. The peculiar emission profiles for Hα were later reproduced in some detail by Grunhut et al. (2012, see their Fig. 15), based on the magnetohydrodynamic (MHD) simulations performed by Sundqvist et al. (2012), where the mass-loss rate was readapted to log $\dot{M}~-7.73$. This means that the emission component of HD 57682 is only reproduced in 2D models. It is important to note that these values from previous authors were all calculated assuming a β-law for the wind velocity profile using the code CMFGEN (Hillier 1990a; Hillier & Miller 1998). In contrast, the theoretical mass-loss rate $\left(\text{log}{\dot{M}}_{\text{sc}}=-6.759\right)$ we derived is able to fit the wings of Hα down to ~ ±100 km s⁻¹ along the line core, even though it is ~ 120 times higher than the initially established value. This shows the difference between velocity profiles that are calculated hydrodynamically with radiative acceleration and jß-law profiles, although it would be necessary to perform self-consistent MHD calculations to study the full shape of Hα in magnetic stars such as HD 57682 in more detail in the future.

Finally, we also mention the issue of the macroturbulence. For their CMFGEN models, Grunhut et al. (2009) adopted a value of v_macro = 40 km s⁻¹, whereas in Grunhut et al. (2017) the authors adopted an even higher value of 65 km s⁻¹. For our spectra, we fit the wings of He I lines using a lower value for the macroturbulence of v_macro = 20 km s⁻¹ ‚ which does not dispose any overestimation of the rotational velocity (we found v_rot = 15kms⁻¹, which is of the same magnitude order as the 13 km s⁻¹ from Grunhut et al. 2012).

Fig. 9 Best fit for HD 57682 with the self-consistent solution tabulated in Table 7.

Fig. 10 Best fit for HD 218915 with the self-consistent solution tabulated in Table 8.

5.4 HD 218915

HD 218915 is a blue supergiant that was classified as O 9.5 lab by Petit et al. (2019), and as O 9.2 lab by Martins (2018). Like 9 Sge, HD 218915 is a runaway star with a periodicity in its radial velocity of P = 0.89 days (Barannikov 2006). From Leitherer (1988), we can obtain the following parameters to set our initial calculations: M_* = 45 M_⊙‚R„ = 26R_⊙ and log L_*/L_⊙ = 5.62, whereas the study performed by Holgado et al. (2018) provides T_eff = 31.1 kK and log g = 3.21. The final parameters we calculated for the self-consistent solution are tabulated in Table 8, and the spectral fit is shown in Fig. 10.

The most prominent aspect to be noted from the new parameters we derived from our self-consistent procedure is the apparent reduction of the He to H ratio to 0.07: this is ~I8% below the new ratio established by solar metallicity by Asplund et al. (2009). This is a remarkable situation because it is commonly assumed that Galactic stars are born with solar metallicity, and the abundance of helium in their surface and wind will only be enhanced through their evolutionary tracks due to nucleosynthesis (Ekström et al. 2012; Gormaz-Matamala et al., in prep.), implying therefore that HD 218915 might have been born with a peculiarly low He/H ratio. Nevertheless, this modification was absolutely necessary to fit the faint absorption profiles of the helium lines, as the decrement on surface was saturating the lines. This was confirmed after we ran several models to scan the entire parameter space (including v_turb). An alternative explanation for the depleted helium of HD 218915 is binarity: the line might become diluted by the light of a potential companion. However, because this is a runaway star, this hypothesis does not seem to apply in this particular case. The implications of this behaviour deserve a particular study, which is beyond the scope of this work. Regarding other stellar parameters, both T_eff and log g lie below the threshold rigidly stated in Paper I (same as the supergiant 325-I from Table 2), but even so, an adequate spectral fit for the star was achieved, therefore we can establish self-consistent wind solutions under m-CAK prescription remain a valid framework for effective temperatures and surface gravities below 32 kK and 3.4.

Similarly as HD 57682, HD 218915 exhibits an emission component in Hα. This might be an indication for the existence of circumstellar material surrounding the star, again due to magnetic activity. However, studies from Grunhut et al. (2017) and Petit et al. (2019) did not find evidence of a strong magnetic field for HD 218915, so that the origin of this additional emission component in Hα remains unclear. This additional component has not been reproduced before, even using a β-law (Holgado et al. 2018). Therefore, its origin can be attributed to the variability of the star. Although it has been possible to accurately reproduce an Hα emission component in another O9.5 I star using synthetic spectra (see e.g. the case of HD 188209 in Martins et al. 2015b), this was achieved using an extremely high value of β = 2.2 for the velocity profile (far from the β = 0.8−1.2 used by Holgado et al. 2018). In the case of HD 188209, the rotational velocity was also lower (v sin i = 45 km s⁻¹). Because Hß and Hγ also show a weak additional component at the core, the wing profiles of these Balmer lines were fitted. The self-consistent mass-loss rate, log M_sc = −6.015 agrees surprisingly well with the value provided by Leitherer (1988). This might partly be attributed to the reduction of the stellar radius made by us from 26 down to 18 R_⊙, which is crucial to reduce the intensity of He I lines. The broadness of He I profiles also was satisfied by increasing the rotational velocity up to 80 km s⁻¹, in contrast with the 60 km s⁻¹ reported by previous studies (Leitherer 1988; Grunhut et al. 2017; Petit et al. 2019) and especially assuming a macroturbulence of 15 km s⁻¹, far below the 94 km s⁻¹ reported by Grunhut et al. (2017).

5.5 HD 195592

HD 195592 is another runaway star, classified as O 9.7 Ia (Sota et al. 2011). It presents evidence of being a binary system, which can explain the line profile variations (LPV) observed for the spectrum (De Becker et al. 2010). The initial mass M_* = 30 M_⊙, radius R_* = 30 M_⊙, and luminosity L_* = 3.1 × 10⁵ L_⊙ were taken from De Becker et al. (2010), who derived them from Martins et al. (2005) given its spectral classification for the closest subtype to HD 195592 (O 9.5 I). We derived our initial parameters also from that calibration, while the final values are tabulated in Table 9 and the synthetic spectra are shown in Fig. 11.

The search for a self-consistent solution for this star led us to find the most extreme parameters of our sample: T_eff = 29.5 kK and log g = 3.2. To achieve this result, it was necessary to increase the number of points for our stratificated hydrodynamic solution to 999 to avoid crashes and false convergences. We call a false convergence the finding of different sets of (k, α, δ) given different initial values (k₀, α₀, δ₀). This situation violates the basic principle of uniqueness for a self-consistent solution: one and only one set of line-force parameters satisfies the e.o.m and recovers the same line-acceleration (see Fig. 2 from Gormaz-Matamala et al. 2019). As we showed in Fig. 4, the selection of the number of points produces significant variations in the resulting wind solution when effective temperature and surface gravity take too low values. Although this discrepancy was observed more over the terminal velocity than over the mass-loss rate, it is evident that the resulting hydrodynamics have an impact on the future synthetic spectra to be generated, which can create a false positive depending on the topology of the solution (Curé & Rial 2007). A deeper analysis is required for a better understanding of the reasons for this exotic behaviour, but this is beyond the scope of this study (see the discussion in Sect. 6.1). The importance of this situation is that the applicability of the m-CAK self-consistent wind prescription is less confident for effective temperatures below 30 kK, and therefore this value could be reconsidered as a new threshold for the recipe. This hypothesis is reinforced by the fit obtained for Hα, where we can infer that the self-consistent velocity profile for a wind solution with a cool temperature of T_eff = 29.5 kK is not adequate to shape the wings and core of the line. Marcolino et al. (2017) provided an even lower effective temperature (T_eff = 28 kK) and a lower surface gravity (log g = 2.9), which confirms that HD 195592 lies outside the range where m-CAK self-consistent solutions can currently be found. In any case, an emission in Hα is at least predicted, and the mass-loss rate provided by $\dot{M}\sqrt{f_{\text{cl}}}$ is relatively close to the homogeneous value provided by Marcolino etal. (2017, log M_hom = −5.14).

5.6 HD 210809

The last star of our sample is HD 210809, another blue super-giant, classified as O 9 Ia (Sota et al. 2011). According to this, the stellar parameters in the catalog of Martins et al. (2005) are 31.3 kK for the effective temperature, 3.23 for log g, and 21.76 for the stellar radius, but the survey of Petit et al. (2019) provided a value of R_* = 15.2R_⊙. Final parameters calculated for the self-consistent solution are tabulated in Table 10, whereas spectral fit is shown in Fig. 12.

The rotational velocity was assumed to be 100 km s⁻¹ (Grunhut et al. 2017). This makes HD 210809, together with HD 192639, the fastest rotators of our sample. This is also clear in the P Cygni profile observed for the Hα line in both stars. Hence the fit of Hα is focused on the emission component, as for the HD 192639 case (Fig. 7). Although we do not reproduce the absorption component of Hα, the accurate fit achieved for the rest of the Balmer series makes us confident about the solution found for HD 210809.

The self-consistent value for mass-loss rate found for HD 210809, $\dot{M}=2.55\times {10}^{-6}{M}_{\odot }{\text{yr}}^{-1}$, is ~30% lower than the value of log $\dot{M}=-5.44$ from Fullerton et al. (2006), but it agrees with the theoretical value from the usage of the Vink formulae with the stellar parameters tabulated in Table 10 $\left({\dot{M}}_{\text{Vink}}=2.5\times {10}^{-6}{M}_{\odot }{\text{yr}}^{-1}\right)$. On the other hand, KK17 provided a theoretical value of ${\dot{M}}_{\text{KK17}}=1.9\times {10}^{-7}{M}_{\odot }{\text{yr}}^{-1}$. Therefore, as shown in Table 2, the m-CAK self-consistent mass-loss closely agrees with Vink’s formula for the cases of T_eff ~ 32 kK.

Fig. 11 Best fit for HD 195592 with the self-consistent solution tabulated in Table 9.

6 Discussion

6.1 Multiple self-consistent solution and saturation

As we pointed out in Fig. 4, the self-consistent solutions for line-force parameters depend on the number of depth points, which produce variations in the resulting wind parameters to be calculated. These variations look more profound for the case of low temperatures, where our standard model of T_eff = 32 kK and log g = 3.4 showed large uncertainties in the resulting self-consistent terminal velocity.

This erratic behaviour is also expressed in the uncertainty for the self-consistent mass-loss rates (our most important wind parameter) when the effective temperature and the surface gravity decreases the threshold imposed in Paper I. These uncertainties are expressed in the error bars associated with the derived self-consistent $\dot{M}$, due to a larger contribution of the term ${\text{Δ}}_{N_{\text{points}}}^2$ in Eq. (13), which are now in the order of the errors assumed to be produced by uncertainties in the stellar parameters such as in the case of HD 218915 (Sect. 5.4, see also the considerations incorporated in Sect. 2.3).

However, as we previously pointed out for the analysis of the star HD 195592 (Sect. 5.5), a decrement in the temperature not only produced a broadening in the error bar for the mass-loss rate, but also a problem in finding a unique solution for any number of depth points implemented. These are the so-called false convergences, and they are better illustrated in Table 11, where we present different sets of self-consistent line-force parameters (k, α, δ) given different numbers of depth points. When we use T_eff < 30 kK and log g = 3.2, we observe more than one solution for the case of the lower number of points (200). This situation does not appear with a larger grid for the wind structure. Even though false solutions disappear when the number of depth points is increased, however, we still observe a huge discrepancy between the sets of (k, α, δ) calculated, which makes the error bar in these cases ~0.5 (such as for HD 195592). A more detailed analysis of the coupling of the line-acceleration and the topology of hydrodynamics is required to understand the phenomenon that is generated for temperatures <30 kK, but this is beyond the scope of this work. At this point, we can state that the discrepancies in the resulting self-consistent wind parameters cause us to set our threshold of the validity for our m-CAK self-consistent prescription at T_eff = 30 kK and log g = 3.2.

Another important case to consider is HD 210809 (Sect. 5.6). The effective temperature is above our threshold of 30 kK, but we still found huge divergences in the self-consistent solutions when we modified the number of depth points. Particularly for HD 210809, we find that these large differences appear when the rotation velocity reaches the value Ω = 0.2. This is better displayed in Table 12, where below this value, variations in the line-force parameters lie within the normal discrepancies due to the selection of different number of depth points; but above this limit, every set of self-consistent k, α, and δ changes dramatically. Moreover, for Ω > 0.2, increasing the number of the depth points always produces the same self-consistent solution where the δ parameter reaches the value of 0.204, which makes it closer to the so-called ‘delta-slow solutions’ (Curé et al. 2011; Venero et al. 2016; Araya et al. 2021). We call this a ‘saturated’ solution because it is independent of the initial value for Ω > 0.2 and because the resulting wind solution (with a mass-loss rate of $\dot{M}\simeq 1.4\times {10}^{-5}{M}_{\odot }{\text{yr}}^{-1}$) cannot generate a synthetic spectrum. Therefore, analogously to the previous case of HD 195592, we found a region in the set of initial stellar parameters in which the self-consistent m-CAK prescription presents problems to provide a unique well-converged solution, this time related with a high rotational velocity. However, given that HD 210809 is not our fastest rotator of the sample (the fastest is HD 192639 with Ω = 0.26, see Table 5), we guess that its low effective temperature (T_eff = 31.5) might act as another threshold for the cases when the rotational velocity becomes significant. In a forthcoming paper, we expect to explore the effects of the rotational velocity on the resulting self-consistent solutions with effective temperatures below 32 kK in more detail. Hence, we consider that for stars with Ω > 0.2, the range of validity of our m-CAK prescription is restricted to T_eff ≳ 31.5 kK.

Nevertheless, despite the false convergence and saturation discussed in this section, we rely on the final parameters introduced for HD 195592 and HD 210809 in Tables 9 and 10, respectively. They represent the best global spectral fit we achieved for each star, as shown in Figs. 11 and 12.

Fig. 12 Best fit for HD 210809 with the self-consistent solution tabulated in Table 10.

6.2 Limits of the self-consistent m-CAK prescription and the issue of clumping

The work introduced in Paper I has provided a grid of wind solutions to self-consistently calculate the line-force parameters for the line-acceleration with hydrodynamics. Because line-force parameters come from the m-CAK framework, these wind solutions are calculated following the formal assumptions of the theory such as the Sobolev approximation, and neglect the detailed influence of multi-scattering effects. Despite these issues, the m-CAK prescription from Paper I has been shown to provide values for M in agreement with homogeneous mass-loss rates derived from observations (Bouret et al. 2005; Markova et al. 2018), in addition to the calculation of reliable synthetic spectra to fit clumped stellar winds (see Sect. 6.2 from Gormaz-Matamala et al. 2019). Now, the inclusion of a stratified temperature for the stellar wind, based on the prescription used by Sundqvist et al. (2019) and Lucy (1971), has led to a decrement in the self-consistent values for the mass-loss rate, particularly for effective temperatures of T_eff ~ 40 kK or higher, mostly because of the decrement in the line-force parameter k. Decrements in $\dot{M}$ however, still agree with theoretical values presented by other self-consistent studies such as KK17 (Krtička & Kubát 2017), as shown in Table 2. On the one hand, this would be expectable if we were to consider that both studies used flux fields from TLUSTY (Hubeny & Lanz 1995). On the other hand, the fact that the results are close, even considering that these authors calculated their line-acceleration from CMF radiative transfer equations (whereas we calculate the line-force parameters), may suggest that some of the drawbacks commonly attributed to the m-CAK framework (such as line-blanketing and multi-scattering effects, see Puls 1987) are not relevant when self-consistent calculations are implemented (at least in the range of temperatures we study), as was previously argued in Gormaz-Matamala et al. (2019, Sect. 2).

Moreover, regardless of the differences with Paper I, the implementation of the temperature structure T(r) shows successful spectral fittings with our sample of stars, as was demonstrated in Sect. 5. This could be discussed in favour of the correctness of the self-consistent values for $\dot{M}$ presented in this work, but in this case, we would also require a proper theoretical determination of the clumping factor for each model. It has been shown that the implementation of different clumping factors leads to different line-accelerations, which in turn will lead to a different self-consistent solution for the wind (Gormaz-Matamala et al. 2021). Including these effects on the calculated gline from the line-force parameters would require a meticulous description of the inhomogeneities through the wind structure, such as the 3D Monte-Carlo calculations developed by Šurlan et al. (2012, 2013). However, these calculations also use a β-law prescription for their velocity profiles, so that a new iterative procedure would be required to couple the two calculations and then to obtain a self-consistent clumping description and υ(r). We expect to move in this direction in a forthcoming study, but this is currently beyond of the scope of this work.

Even when a fully consistent calculation of every parameter of the stellar wind is the most desirable scenario and the final goal of every theoretical prescription, in practice, this achievement cannot be reached. The spectral fits we presented and their respective self-consistent wind parameters ${\dot{M}}_{\text{sc}}$ and υ_∞,sc were made on the basis that the fitted stellar parameters are correct, especially because the error bars associated with each wind parameter do not consider the uncertainties of temperature, mass, radius, or abundances. This is the price to pay for reducing the number of free parameters for the spectral fit. Therefore, we do not claim that our predicted mass-loss rate and terminal velocity are ‘the truth’, but we state that our prescription leads to reliable theoretical wind parameters that can adequately reproduce (despite particular exceptions such as the Hα profile of HD 218915) spectral observations and agree with other studies. Based on these results, self-consistent m-CAK wind solutions show a promising perspective for the theoretical determination of mass-loss rates, and hence we are confident in the potential of this prescription in the future.

Fig. 13 Comparison of the best fits for 9 Sge (left-panel) and HD 210809 (right-panel), with the self-consistent solution introduced in Fig. 8 for 9 Sge in Fig. 12 for HD 210809 (red) and with the improved fits individualised for Hα (blue).

6.3 Accuracy of spectral fitting

In Sect. 5 we introduced a set of stellar and wind parameters to simultaneously fit a sample of optical spectra observed with the HERMES spectrograph. The fit was achieved by eye, and our main criterion was the global accuracy for the lines tabulated in Table 4. However, despite our general focus on providing accurate fits for most of the listed lines (mainly the photospheric lines), the quality the fits achieved for some of the Hα profiles in our sample might not be considered satisfactory.

Improved fits for the Hα lines of 9 Sge and HD 210809 are shown in Fig. 13. For the case of 9 Sge, we can improve the fit of Hα if we increment the turbulence velocity to 40 km s⁻¹ and reduce the rotational velocity to 50 km s⁻¹. In addition, the macroturbulence needs to be decreased to 4 km s⁻¹ in order to properly reproduce the sharpness of the Hα profile. This improvement produces a sharper emission component, even when it was not possible to reproduce the absorption portion of the P-Cygni profile of the line. For the case of HD 210809, a more accurate fit is achieved also by setting υ_turb = 40 km s⁻¹, decreasing υ_macro = 2 km s⁻¹ and increasing the rotational velocity to υ_rot = 170 km s⁻¹. In this case, the absorption and emission components of the line are accurately reproduced, but these new modifications do not fit the other spectral lines.

Therefore, the fits introduced in Sect. 5 are the best global fits provided by our method when we aim to obtain a simultaneous fit for the visible range of the spectra rather than focus only on one single line such as Hα. Studies focused on the self-consistent description of the wind either do not provide a spectral fit (Krtička & Kubát 2017; Björklund et al. 2021) or the presented spectral fits do not allow a flexible fine-tuning (Sander et al. 2017). The reason for this dichotomy is, as stated in Gormaz-Matamala et al. (2021), that a perfect fully self-consistent wind prescription must consider every feature present in the wind (e.g. the clumping prescription and even effects from X-rays), which is beyond the scope of this work. In the case of our m-CAK prescription, even when the selection of the wind parameters is based on a unique solution coming from stellar parameters (and then they are not free parameters), the subsequent spectral fit is capable to adequately reproduce the lines in the optical range. Some lines cannot be adjusted, such as the Hα emission components for HD 57682 and HD 218915, because spherically symmetric codes were used, but this is not a failure attributable to the m-CAK method. As shown by Grunhut et al. (2012), CMFGEN is also not able to fit these additional emission components. The exception is the poor fit for Hα in HD 195592, but as discussed in Sect. 5.5, this is an indicator that the limits of our prescription coincide with the lower threshold of T_eff ≃ 30 kK, as stated in Sect. 2.4. Because of these reasons, we are confident that our spectral fitting using the m-CAK prescription is reliable, and therefore the stellar and wind parameters determined by our method are trustworthy. We expect in the future to extend this analysis to UV spectra, where P Cygni profiles such as O V λ1371, C IV λ1548, and N IV λ1718 will be useful to evaluate the terminal wind velocities and clumping factors (Bouret et al. 2012; Šurlan et al. 2012).

6.4 Self-consistent WML relation

The wind-momentum-luminosity (WML) relation predicts a strong dependence of the wind momentum rate on the stellar luminosity with an exponent determined by the statistics of the strengths of the hundreds of thousands of lines driving the wind (see Puls et al. 2000). It also contains a weak dependence on stellar radius. This is due to the work against the gravitational potential that the stellar wind must overcome. This relation has been confirmed observationally using objects from the Milky Way, the LMC, SMC, and other galaxies of the local group by Bresolin et al. (2001, 2003), Urbaneja et al. (2002, 2003), Evans et al. (2005, 2007), and Mokiem et al. (2005, 2006).

According to Kudritzki & Puls (2000), D_mom is defined as $$D_{\text{mom}}=\dot{M}{\upsilon }_{\infty }\sqrt{R_{*}/R_{\otimes },}$$(15)

and the WML reads $$\log D_{\text{mom}}=\log D_0+x\log (L_{*}/L_{\odot }).$$(16)

Figure 14 shows the WML for the sample of O I stars from Kudritzki & Puls (2000) as grey circles and our results as red squares. Additionally, we included D_mom values for the sample of field stars from Bouret et al. (2012), whose values for mass-loss rate and terminal velocities are based on spectral fitting; and Markova et al. (2018), which are based on unclumped values for the mass-loss rate and hence were previously used to compare with self-consistent solutions (Gormaz-Matamala et al. 2019, Fig. 13). Although our sample is reduced, our log D_mom clearly accurately matches the values from Kudritzki & Puls (2000). To make this comparison coherent, we used unclumped values for our self-consistent mass-loss rates following $$D_{mom,\mathrm{hom}}={\dot{M}}_{\text{clump}}{f}_{\text{cl}}^{1/2}{\upsilon }_{\infty }\sqrt{R_{*}/R_{\odot }}.$$(17)

Finally, Table 13 compares the values of log D₀ and x (Eq. (16)) from Kudritzki & Puls (2000) and from Markova et al. (2018) with our results. We observe that both Bouret et al. (2012) and Markova et al. (2018) list a deviation from the original fit of Kudritzki & Puls (2000), but are still within the limits of the standard deviations of log D₀ and x. On the other hand, our homogeneous values for D_mom follow almost the same trend as in the previous studies. This is particularly clear in Fig. 14, where empty squares (clumped $\dot{M}$) lie slightly out of the main band, but filled squares (unclumped $\dot{M}$) lie inside it. Despite these encouraging results, the fits introduced in Table 13 represent a simple linear regression based on a few points, and we certainly need to expand our sample and extend our analysis to UV wavelengths to provide a more accurate conclusion.

Fig. 14 Wind momemtum, Dmom for the stars of our sample as a function of their luminosity. Filled red squares represent the momentum derived from a homogeneous mass-loss rate (see Eq. (17)), and empty red squares represent the clumped values for $\dot{M}$. In addition, we over-plot the results for O I stars from Kudritzki & Puls (2000, grey circles, from their Fig. 8), the Galactic O-supergiants from Bouret et al. (2012, green circles), and the field stars taken from Markova et al. (2018, blue circles). The respective linear regressions from Table 13, following the same colour convention, are also included.

7 Conclusions

We have exploited the potential of the m-CAK prescription performed by Gormaz-Matamala et al. (2019, Paper I) to derive self-consistent wind parameters by performing spectral fittings for a sample of six O-type stars (HD 192639, 9 Sge, HD 57682, HD 218915, HD 195592, and HD 210809) observed with the high-resolution echelle spectrograph HERMES. In addition, the m-CAK prescription has been improved to include a temperature structure throughout the wind, together with a more exhaustive quantitative implementation of the uncertainties associated with each model of self-consistent solution.

In a broader context, self-consistent synthetic spectra agree well with HERMES observations, which are no longer based on a β-law for their velocity profiles. Improvements could be made using this method to fit Hα in particular, but this would undermine the quality of the fit for the remaining lines. This success reinforces the validity of the set of new stellar and wind parameters introduced in this work for our sample of stars, especially when there are significant differences with previous literature. These differences can be attributed to the previous usage of a β-law, in contrast with our self-consistent hydrodynamics calculations. Nevertheless, our method treats clumping as a free parameter that is only modified when the line fits are fine-tuned (Sect. 4.2), even when it is known that clumps affect the line-acceleration (Gormaz-Matamala et al. 2021) because inhomogeneities are currently not included when the line-force parameters are calculated (Sect. 2). A full implementation of clumping would require a thorough analysis to describe the wind inhomogeneities (Šurlan et al. 2012, 2013). In addition, the inclusion of UV spectral analysis will complement the constraint of mass-loss rate and clumping. However, this is beyond the scope of this work, and here we just introduce the first attempt to derive stellar and wind parameters for O-type stars through a spectral analysis considering the self-consistent wind hydrodynamics. We expect that the values introduced here are helpful for future studies of all the stars constituting this sample.

The qualitative diagnostics and the detailed analysis of the self-consistent wind solutions obtained for some stars establish that the validity of our m-CAK prescription applies to O-type stars with T_eff ≥ 30 kK and log g ≥ 3.2 (although the threshold for effective temperature should increase if we consider rotational velocity Ω > 0.2). We remark that this range of validity is where the line-force parameters (k, α, δ)_sc can be fitted from the force multiplier $\mathcal{M}\left(t\right)$ assuming them to be constants instead as a function of radius, because their variance is below the common uncertainties associated with stellar parameters. Subsequently, these thresholds indicate that self-consistent solutions can describe the winds of massive stars from their birth on the main sequence and cover then a large part of their lifetime in this stage before their evolution into later-type stars. This aspect is incorporated in Gormaz-Matamala et al. (in prep.).

Acknowledgements

We sincerely thank the anonymous referee for the valuable feedback and comments provided to us. We also thank J. Puls for helpful discussions that improved this work and for having put at our disposal his code FASTWIND. This project has received funding from the European Unions Framework Programme for Research and Innovation Horizon 2020 (2014–2020) under the Marie Skłodowska-Curie grant Agreement N° 823734. A.C.G.M. has been financially supported by the PhD Scholarship folio N° 2116 1426 from National Commission for Scientific and Technological Research of Chile (CON-ICYT). A.C.G.M. and M.C. acknowledge support from Centro de Astrofísica de Valparaíso. M.C., I.A. and C.A. thank the support from FONDECYT project 1190485. A.L. acknowledges in part funding from a ESA/PRODEX Belgian Federal Science Policy Office (BELSPO) contract related to the Gaia Data Processing and Analysis Consortium. J.C. thanks the supports from FONDECYT project 1211429. J.A.P. has been funded by grants from Universidad Nacional de La Plata, Argentina (Proyecto I+D 11/G162). I.A. thanks the support from FONDECYT project 11190147. C.A. thanks to FONDECYT project N. 11190945.

Appendix A Tau

The relation between the radiative acceleration and $$g_{\text{rad}}(r)=\frac{{\kappa }_{\text{F}}{L}_{*}}{4\pi c{r}^2}=\frac{{\kappa }_{\text{F}}F(r)}{c},$$(A.1)

with $${\kappa }_{\text{F}}F(r)={\displaystyle {\int }_0^{\infty }{\kappa }_v{F}_{v}dv.}$$(A.2)

$$\begin{array}{c}{\tau }_{\text{F}}(r)={\displaystyle \int {\kappa }_{\text{F}}\rho (r){\left(\frac{R_{*}}{r}\right)}^{2}dr,}\\ ={\displaystyle \int \frac{c{g}_{\text{red}}}{F(r)}\rho (r){\left(\frac{R_{*}}{r}\right)}^{2}dr,}\\ ={\displaystyle \int g_{\text{red}}\frac{4\pi c{r}^2}{L_{*}}\rho (r){\left(\frac{R_{*}}{r}\right)}^{2}dr,}\\ =\frac{4\pi c{R}_{*}^2}{L_{*}}{\displaystyle \int g_{\text{red}}\rho (r)dr,}\\ =\frac{4\pi c{R}_{*}^2}{L_{*}}{\displaystyle \int (g_{\text{line}}+g_{\text{es}})\rho (r)dr,}\\ =\frac{4\pi c{R}_{*}^2}{L_{*}}{\displaystyle \int (g_{\text{line}}+g_{\text{grav}}{\Gamma }_{\text{es}})\rho (r)dr.}\end{array}$$(A.3)

Therefore, for the temperature field, $$\begin{array}{c}T(r)=T_{\text{eff}}{\left[W(r)+\frac{3}{4}{\tau }_{\text{F}}\right]}^{1/4}\\ =T_{\text{eff}}{\left[W+\frac{3\pi c{R}_{*}^2}{L_{*}}{\displaystyle \int (g_{\text{line}}+g_{\text{grav}}{\Gamma }_{\text{es}})\rho dr}\right]}^{1/4}\end{array}$$(A.4)

To numerically solve the integral, see Appendix B.

Appendix B Numerical integral

We have the integral $${\displaystyle {\int }_r^{\infty }(g_{\text{line}}+g_{\text{grav}}{\Gamma }_{\text{es}})\rho (r^{\prime })d{r}^{\prime }.}$$(B.1)

Being x(r) = g_rad ρ(r) and X = ∫ xdr, the integral can be rewritten in a discrete way as $$X_i={\displaystyle \sum _i^N\frac{x_i+x_{i+1}}{2}}\times (r_{i+1}-r_i).$$(B.2)

The problem with this formulation is that for the last item of the grid (i.e. when i = N), ∆r has no definite value because the limit is the infinite. In this case, it is better to perform a change of variables over the integral. Using the auxiliar variable u = −1/r and therefore $${\displaystyle {\int }_{-1/u}^{\infty }(g_{\text{line}}+g_{\text{grav}}{\Gamma }_{\text{es}})\rho (u^{\prime })\frac{d{r}^{\prime }}{d{u}^{\prime }}d{u}^{\prime }.}$$(B.3)

Discretisation is then $$X_i={\displaystyle \sum _i^N\frac{x_i+x_{i+1}}{2}\frac{u_{i+1}-u_i}{u_i^2}.}$$(B.4)

References

All Tables

All Figures

	Fig. 1 Velocity profile generated from a self-consistent wind solution using the m-CAK prescription, compared with some v(r) following β-law. The bottom panel corresponds to a zoom into the base of the wind.
In the text

	Fig. 2 Comparison of self-consistent fitting of the force multiplier M(t) between the implementation of the temperature structure (blue) and the fixed temperature (red). Calculated line-force parameters (k, α, δ) are displayed in the same colour code. The shaded area represents the t optical depth range where the force multiplier is defined within, from the sonic point outwards.
In the text

	Fig. 3 Temperature stratification for HD 57682 (leftpanel, Teff = 35 kK) and HD 210809 (rightpanel, Teff = 31.5 kK), either obtained from Eq. (6) (red) or calculated by FASTWIND (blue).
In the text

	Fig. 4 Comparison of the obtained self-consistent wind parameters v∞ and M, employing different numbers of depth points N to numerically solve the equation of motion. Values are expressed as departures from the self-consistent terminal velocity and mass-loss rate when the maximum number of points is allowed for HYDWIND (999).
In the text

	Fig. 5 Theoretical mass-loss rates tabulated in Table 2 for the values predicted by Vink’s formula (orange symbols), by KK17 (green symbols), and by this study (blue symbols). Main-sequence stars of the sample are represented by starss, giants by triangles, and supergiants by circles.
In the text

	Fig. 6 Scheme of the self-consistent m-CAK procedure combined with FASTWIND, showing the stage at which the different initial parameters (in rectangles) are set. Ellipses represent the processes or codes to be run, and the diamonds represent the outputs.
In the text

	Fig. 7 Best fit for HD 192639 with the self-consistent solution tabulated in Table 5.
In the text

	Fig. 8 Best fit for 9 Sge with the self-consistent solution tabulated in Table 6.
In the text

	Fig. 9 Best fit for HD 57682 with the self-consistent solution tabulated in Table 7.
In the text

	Fig. 10 Best fit for HD 218915 with the self-consistent solution tabulated in Table 8.
In the text

	Fig. 11 Best fit for HD 195592 with the self-consistent solution tabulated in Table 9.
In the text

	Fig. 12 Best fit for HD 210809 with the self-consistent solution tabulated in Table 10.
In the text

	Fig. 13 Comparison of the best fits for 9 Sge (left-panel) and HD 210809 (right-panel), with the self-consistent solution introduced in Fig. 8 for 9 Sge in Fig. 12 for HD 210809 (red) and with the improved fits individualised for Hα (blue).
In the text

Fig. 14 Wind momemtum, Dmom for the stars of our sample as a function of their luminosity. Filled red squares represent the momentum derived from a homogeneous mass-loss rate (see Eq. (17)), and empty red squares represent the clumped values for $\dot{M}$. In addition, we over-plot the results for O I stars from Kudritzki & Puls (2000, grey circles, from their Fig. 8), the Galactic O-supergiants from Bouret et al. (2012, green circles), and the field stars taken from Markova et al. (2018, blue circles). The respective linear regressions from Table 13, following the same colour convention, are also included.

In the text