© The Authors 2023

1 Introduction

The lives of massive stars are largely controlled by their strong radiation-driven winds. The resulting mass-loss has an impact on, for example, the stellar luminosity, its lifetime, and apparent temperature and consequently also on the ionising radiation and UV luminosity from these massive stars (Meynet et al. 1994; Smith 2014). These winds also shape the interstellar medium (ISM) providing chemical and mechanical feedback to their surroundings. The winds sweep up the ISM creating giant bubbles in the star’s surroundings and the more massive stars even insert more energy into the ISM than supernovae (SNe; Fierlinger et al. 2016). Additionally, the amount of mass-loss largely determines which type of SN the star explodes as, as well as the final mass of the remnant after that explosion (Fryer & Kalogera 2001; Belczynski et al. 2020).

This last effect plays a crucial role in explaining the masses of black holes as detected by gravitational waves (GWs) that arise from the merging of a binary black holes. The GW detections by the Advanced-LIGO interferometers provided this first observational evidence of black hole mergers (Abbott et al. 2016) with black hole masses of 36 M_⊙ and 29 M_⊙ prior to merging. The surprisingly high observed masses of these black holes as well as those from more recent observations point to the need for weak stellar winds. Indeed, mass-loss is highlighted as a main uncertainty in massive-star progenitor models of such GW sources (Abbott et al. 2016).

Accurate predictions of the mass-loss of massive stars are thus necessary for predicting their evolutionary paths (Smith 2014). Such predictions are included in stellar-evolution codes (such as the Geneva stellar-evolution code GENEC, Eggenberger et al. 2008; STERN, Heger et al. 2000, Brott et al. 2011; and MESA, Paxton et al. 2019¹) through simple scaling recipes that depend on fundamental stellar parameters. The prescriptions used in most massive-star evolution models are typically derived from either theoretical models or empirical studies that focus on certain parts of the Hertzsprung-Russel diagram (HRD), where a distinction is usually made between winds from hot stars (Vink et al. 2001; Nugis & Lamers 2000; Björklund et al. 2021) and cool stars (Reimers 1975; de Jager et al. 1988; Bloecker 1995; Kee et al. 2021).

In the first two papers of this series, we developed new atmosphere and wind models in order to provide mass-loss predictions and predicted wind structure (Sundqvist et al. 2019, hereafter Paper I) and applied them to a sample of O stars (Björklund et al. 2021, hereafter Paper II). The wind models solve the spherically symmetric, steady-state equation-of-motion (eom), relying on non-local thermodynamic equilibrium (NLTE) radiative transfer in the co-moving frame (CMF) to calculate the radiative acceleration throughout the atmosphere and wind outflow (Puls et al. 2020). In the O-star regime, the results of Paper II suggest an overall reduction in the mass-loss rate of O stars by a factor of ~2–3 compared to the rates normally applied in stellar-evolution calculations (Vink et al. 2000, 2001), in good agreement with a range of observational results (Najarro et al. 2011; Sundqvist et al. 2011; Bouret et al. 2012; Šurlan et al. 2013; Cohen et al. 2014; Shenar et al. 2015; Hawcroft et al. 2021; Rubio-Díez et al. 2022).

In this paper, we further study the dependence of the mass-loss rate on important stellar parameters, now also including the stellar mass and effective temperature on top of the stellar luminosity and metallicity. This is achieved by constructing a large grid of atmosphere and wind models that covers the O-star main-sequence and early post-main-sequence evolution of O-type stars and now also cooler B-type stars. This allows us to study the behaviour of the mass-loss over the parameter ranges and to derive a new mass-loss recipe. This includes a comparison to previous predictions as well as a study of the effect of the bi-stability region (a region in effective temperature where important chemical elements that drive the wind, such as iron, change their ionisation stage) around 22 kK (e.g. Vink et al. 2000, 2001; Driessen et al. 2019; Krtička et al. 2021). The effect of the new mass-loss rates on massive-star evolution is analysed with a focus on the creation of Wolf-Rayet (WR) stars and high-mass black holes.

2 Wind models

The wind models used here rely on the interplay between the radiative transfer and the atmospheric and wind hydrodynamics. A converged model is hydrodynamically consistent which means that the one-dimensional eom in the spherically symmetric, steady-state case is self-consistently solved. In the case of radiation-driven stellar winds, the most important contribution to the eom comes from the radiative acceleration g_rad(r), (mostly due to lines), computed for each radial point from detailed radiative transfer calculations relying on the NLTE radiative transfer of FAST WIND using a CMF technique (Paper I and II, Puls et al. 2020).

The initial guess of the velocity structure is taken to be a ß-velocity law, where the steepness of the acceleration towards the terminal wind speed is set by the exponent ß, which is stitched to a quasi-hydrostatic atmosphere in the inner region (see Eq. (3) of Paper II). As for Ṁ, we start from a first guess, for which we use a third of that predicted by the Vink et al. (2001) recipe as justified by the results from Paper II. From test calculations, we have verified that the exact initial conditions do not significantly affect our finally predicted model-structures; however, if these initial conditions are (very) poorly chosen this can affect the convergence behaviour of the simulation, slowing it down or even leading to non-convergence. During the radiative transfer the radiative acceleration is converged in an iterative process updating the occupation numbers of the included atoms. Afterwards the eom is solved by integrating in- and outwards from the sonic point from which a new velocity-structure is obtained. From the mismatch in the force balance at the sonic point, we compute a new mass-loss rate that would counter this mismatch following a basic relation $\dot{M}\propto \frac{1}{g_{\text{rad}}}$ (see Paper I). Then from υ(r) and Ṁ a new density structure is computed using the continuity equation. Additionally, the temperature structure is updated in the same way as in Paper II (see their Eqs. (6)–(7)), following Lucy (1971). As also discussed in Paper II, adopting this temperature structure leads to much better and faster model-convergence, and does not introduce significant differences in the dynamics as compared to models with a temperature structure derived from flux conservation and thermal balance of electrons (see also Fig. 5 in Paper I). Throughout the iteration procedure the stellar luminosity L_* and mass M_* are kept fixed. The radius is updated such that the final converged model has R_* ≡ r(r_F = 2/3), where R_* is the input stellar radius and τ_F the spherically modified flux-weighted optical depth (see Paper I). All models are calculated from deep optically thick layers (lower boundary typically at a column mass of 80) to a large wind radii of about R_max ~ 100 R_*, where the winds have reached their terminal wind speed υ_∞ ≡ v(r = R_max).

2.1 Convergence

The radiative transfer and the eom influence each other, such that we need to iterate between the two steps until a converged solution is reached. As already mentioned above (see also Papers I and II), in each hydrodynamic iteration cycle the NLTE occupation numbers and the radiation field are updated until g_rad(r) is converged for the given structure. The resulting radiative accelerations and flux-weighted opacities are then used to update the velocity and temperature structures, and a new density obtained by following the procedure above for updating the mass-loss rate. This new hydrodynamic structure (ρ(r), v(r), T(r)) is then used to again converge the radiative acceleration (g_rad(r)) by means of a new solution of the NLTE equations and the CMF radiative transfer. This procedure is iterated until convergence, which then also results in final convergence of the temperature structure, ionisation balance and NLTE occupation numbers (the last of which typically to within a maximum error of a few percent). In our hydrodynamic models final convergence is defined through the maximum error in the eom (see Papers I and II). This error compares Γ = g_rad/g to the other terms in the eom according to $$f_{\text{err}}\left(r\right)=1-\frac{\text{Λ}\left(r\right)}{\text{Γ}\left(r\right)},$$(1)

where $$\text{Λ\hspace{0.17em}}=\frac{1}{g}\left(\upsilon \frac{\text{d}\upsilon }{\text{d}r}\left(1-\frac{a^2}{{\upsilon }^2}\right)+g-\frac{2a^2}{r}+\frac{\text{d}{a}^2}{\text{d}r}\right).$$(2)

Our models are formally defined as converged when the maximum error radially in f_err $$f_{\text{err}}^{\max }=\max \left(\text{abc}\left(f_{\text{err}}\right)\right)$$(3)

is less than a certain threshold throughout the complete atmosphere. Typically, we set this limit to a very stringent 2%; however, in certain regions of the HRD we allowed for somewhat larger deviations (see the discussion below). Additionally, we required both the mass-loss rate and the velocity structure not vary by more than 2% and 3%, respectively, between the last two hydrodynamic steps.

2.2 Clumping

We have not included clumping or effects from X-rays in the models. As discussed in Paper I and Paper II, the inclusion of clumping (albeit in a highly simplified and somewhat inconsistent way, see the discussion below) with an onset beyond the sonic point typically does not significantly affect the predicted mass-loss rates, because they are mostly sensitive to the conditions locally around the sonic point. This is different than, for example, models that derive a mass-loss rate only from a global energy balance constraint and not from a locally consistent solution to the eom (Muijres et al. 2011). The recent simulations by Krtička et al. (2021) also seem to find a general dependence on clumping, but again these models are somewhat different than ours since their Sobolev-scaled line force shifts their wind critical point downstream from the sonic point to the supersonic parts of the outflow. By contrast, in our models the main ‘choke’ of the wind, which effectively controls the rate of mass-loss, is in principle always located in the near-sonic regions (see the discussions in Papers I and II). Regarding potential inhomogeneities in sub-sonic optically thick layers, this cannot really be treated using any of the current methods for modelling clumping in hot-star atmospheres with winds. As such, potential effects of sub-surface structures upon line-driven mass-loss are currently very difficult to assess.

Moreover, we will also point out here that current inclusions of wind clumping into one-dimensional line-driven steady-state hydrodynamic simulations (such as those presented here) are of a quite ad hoc and somewhat inconsistent nature. For example, the standard way of introducing such wind clumping in corresponding models used for spectroscopic analysis is to assume that all wind mass is contained within optically thin clumps (see Puls et al. 2008 for a review). This can then shift the wind ionisation balance and thus also have an implicit effect on g_rad. But in a steady-state wind model, g_rad is rather computed for an assumed smooth (or average) wind, and not for the clumps themselves. Indeed, due to the strong inverse dependence of line-driving with density (essentially g_rad ~ 1/ρ^a where 0 < a < 1, Castor et al. 1975), test calculations reveal that the clumps themselves actually do not experience significant line-driving. As such, in order to take clumping into account in hydrodynamic models one must really follow the dynamics of the (time-dependent and possibly multi-dimensional) flow in a way that accounts properly for the very different line-driving experienced by high- and low-density gas parcels. But with the elaborate techniques used to compute g_rad in this series of papers, this is not computationally possible. Instead, time-dependent line-driven wind models attempting to compute predictions for clumping (e.g., Driessen et al. 2019) typically use a parametrised line force calibrated to reproduce a certain average mass-loss. That is, these models are not able to predict mass-loss rates but rather use simplified radiative transfer in order to make time-dependent calculations feasible, and by means of these then predict the clumpy wind structure. Vice versa, steady-state models (such as those presented here) use very detailed techniques for computing g_rad and Ṁ, on the expensive of having to assume a spherically symmetric and smooth steady-state eom. Indeed, to our knowledge all non-parametrised² steady-state line-driven wind models attempting to include ‘clumping’, have computed the ionisation for the clumps themselves, but then (inconsistently) solved the eom for the smooth density of a steady wind (e.g., our own tests in Paper I; Sander et al. 2020; Krtička et al. 2021). In our opinion, it is questionable if such an approach is able to give any reasonable answers to the question of if or how clumping may affect also the final (average) mass-loss rate of the star. Rather theoretical models targeting simultaneous investigations of clumping and mass-loss rates need to be carried out using time-dependent hydrodynamics, for which one then would need to develop more elaborate (but still sufficiently fast) techniques for computing g_rad than those assumed thus far in such simulations.

Finally, we note that the above (of course) applies to theoretical line-driven wind simulations aiming to predict mass-loss rates. The situation is very different for diagnostic models aiming to obtain empirical mass-loss rates from direct comparison to observations, where it is crucial to properly account for the effects of wind clumping upon the diagnostic under consideration (see also the discussion in Sect. 5.2.2).

3 The new grid

3.1 Grid setup

We set up a grid of radiation-driven wind models that cover the hot region of the HRD where massive stars evolve most of their lives. We excluded the cool temperatures (<15 000 K), where we find, for example, red supergiants that presumably are driven by other processes (e.g. Kee et al. 2021), and the extremely hot temperatures of stripped WR stars (Sander et al. 2020; Poniatowski et al. 2021). For a given metallicity Z_*, in this grid we vary the luminosity of the star, its effective temperature, and its mass. It is not a regular ‘cubic’ grid, because we preferentially include combinations of L_*, T_eff, and M_* that are accessible from typical predictions of stellar evolution. To set this up, we computed evolutionary tracks of stars between 15 and 80 M_⊙ using MESA (Modules for Experiments in Stellar Astrophysics, Paxton et al. 2011, 2013, 2015, 2018, 2019). MESA is a computational tool that allows one-dimensional stellar structure and evolution models to be created for a wide variety of stellar parameters. The ranges covered by the grid for each parameter are 10^4.5–10⁶ L_⊙ for L_*, 15 000–50 000 K for T_eff, and 15–80 M_⊙ for M_*. Figure 1 shows the extent of the model grid on the HRD. Each point represents a few models where the mass is varied according to the position on the HRD. Each model is calculated with a solar metallicity of Z_⊙ = 0.014 and lower metallicities of 0.5 Z_⊙ and 0.2 Z_⊙, which correspond to those of the Large Magellanic Cloud (LMC) and the Small Magellanic Cloud (SMC). This leads to a total of 516 wind models. This setup of the model grid allows us to investigate the dependences on the different stellar parameters such as luminosity and metallicity, mass, and effective temperature. The last two, as discussed in Paper II, needed a wider parameter coverage than the O-star range considered there in order to capture the proper dependence. The models all have a fixed helium number abundance $Y_{\text{He}}=\frac{N_{\text{He}}}{N_{\text{H}}}=0.1$ with also the micro-turbulent velocity kept fixed at υ_turb = 10 km s⁻¹.

Fig. 1 HRD showing the coverage of the grid of models. Each point in the L*-Teff plane represents several models with varying mass, shown by the overlaid squares. Also shown are the stellar tracks including an initial rotation of 0.4 times the critical rotation for 15, 20, 30, 40 and 60 M⊙ stars computed using MESA.

3.2 Regions of the HRD

The grid is chosen such that we can compute steady-state wind models that experience radiation-driven outflows. Hot, main-sequence (and post-main-sequence supergiants that have not evolved too far) O stars fall nicely within this description, and as in Paper II we find that these models normally are quite well behaved and converge in a relatively straightforward way. There are, however, some areas of the HRD where it is not so straightforward to converge the models. Two such regions are discussed below.

Several models with effective temperatures around 20 kK do not reach the strict convergence criterion discussed above, the maximum 2% error for the eom over the complete radial grid, but rather show hydrodynamic solutions displaying somewhat larger maximum errors. It is important to realise, however, that even though we do not always reach this strict criterion, this does not mean that the models in general also fluctuate significantly in their predicted mass-loss rates; rather this reflects the general difficulty of converging the full eom in these types of steady-state line-driven wind models (see also the discussions in Sander et al. 2017; Krtička & Kubát 2010; Gormaz-Matamala et al. 2021). Indeed, in many of these iteration steps the mass-loss rate changes only marginally and $f_{\text{err}}^{\max }<10\%$, with most of the atmosphere showing significantly lower errors. As such, we are still able to use these results for the wind structure and mass-loss rate. In fact, our 2% base-criterion is very stringent; for example, the alternative CMF O-star model by Sander et al. (2017) applied a 5% criterion and also allowed for larger deviations close to the inner and outer boundaries. In models that do not reach this strict criterion, we calculated the final mass-loss rate as an average of all Ṁ in the hydrodynamic iteration steps with $f_{\text{err}}^{\max }<10\%$ (i.e. we compute 〈Ṁ_<10%〉). To justify this averaging approach, we compared the predicted mass-loss rate of all formally converged models (i.e. those which reach Ṁ_<2%) to the Ṁ that would have been derived if we instead had applied the averaging method 〈Ṁ_<10%〉 for also these model stars. From this, we find that the difference $\left|\frac{{\dot{M}}_{<2\%}-\langle {\dot{M}}_{<10\%}\rangle }{{\dot{M}}_{<2\%}}\right|$ is generally very small, <5%; indeed, in our complete grid we only find three outliers, for which the difference is then still a modest ~20%. This comparison supports the use of 〈Ṁ_<10%〉 as the mass-loss rate of models that do not formally reach Ṁ_<2%.

We limit the parameter range of the grid to a maximum luminosity log L_*/L_⊙ = 6. Above this limit, the strong radiation field often induces a radiation force that exceeds gravity already in the dense, sub-surface regions of the star. Such deep-seated outflows cannot be modelled well with our current technique, both because our atomic database does not contain the high ionisation stages that are needed in these hot atmospheric layers, and because a steady-state approach there becomes highly questionable (e.g. Jiang et al. 2018; Schultz et al. 2020). For now we therefore do not include very high-luminosity regions in our model grid, presumably involving stars that would spectroscopically be classified as hydrogen-rich WR stars (WNh, sometimes called ‘slash stars’ or ‘O stars on steroids’) and luminous blue variables (LBVs).

4 The new recipe

The newly computed atmosphere and wind models allow us to analyse trends of the mass-loss rate with the stellar parameters that are varied in the grid. To this end, we aim to provide a simple theoretical prescription of Ṁ that depends on the stellar luminosity, mass, effective temperature and metallicity. For this we assume the form: $$\dot{M}\propto L_{*}^m{M}_{\text{eff}}^n{T}_{\text{eff}}^p{Z}_{*}^q,$$(4)

that is $$\log \dot{M}=C+m\log L_{*}+n\log M_{\text{eff}}+p\log T_{\text{eff}}+q\log Z_{*},$$(5)

where M_eff = M_*(1 – Γ_e) is the effective stellar mass, reduced by the electron-scattering Eddington parameter ${\text{Γ}}_{\text{e}}=\frac{{\kappa }_{\text{e}}{L}_{*}}{4\pi cG{M}_{*}}$ applying a constant K_e in the model-fit. The value for K_e is calculated by $${\kappa }_{\text{e}}=\frac{0.4\left(1+i_{\text{He}}{Y}_{\text{He}}\right)}{1+2Y_{\text{He}}},$$(6)

with i_He the ionisation stage of helium (with i_He = 0 for neutral helium). For O stars, with doubly-ionised helium and our assumed Y_He = 0.1, this gives the typical value of 0.34 cm² g⁻¹. Using the mass-loss rates obtained in our grid we can fit Eq. (5) and derive the best matching, m, n, p, q and C coefficients. As was also found in Paper II, though, one value for the exponent q (metallicity dependence) does not properly represent the grid results, because the parameter itself depends quite strongly on the other stellar parameters. In Paper II we showed that a linear dependence with log L_* can capture the variation in the metallicity dependence across the grid. Having access now to a wider range of parameters for the wind models, this variation is characterised more properly by using log T_eff instead of log L_*. In the previous grid of O stars, this was however essentially the same because of the strong relationship of L_* and T_eff for models within that O-star parameter range. Indeed, when replacing the luminosity by the effective temperature in the exponent q when re-deriving the fitting formula, the predicted Ṁ only changes by at most ~4% to those predicted by Eq. (20) in Paper II. Now, when we turn to the current grid of models we find that q shows the tightest relation with T_eff, which is shown in Fig. 2. The scatter in this figure likely arises due to a combination of less significant dependences on other parameters and the limited number of metallicity points used in our model grid. Finally then, we construct a relation for the mass-loss rate from a multi-linear fit including the Z_*(T_eff) dependence as follows log Ṁ = − 5.52 $$\begin{array}{ll}\log \dot{M}=\hfill & -5.52\hfill \\ \hfill & +2.39\log \left(\frac{L_{*}}{{10}^6{L}_{\odot }}\right)\hfill \\ \hfill & -1.48\log \left(\frac{M_{\text{eff}}}{45M_{\odot }}\right)\hfill \\ \hfill & +2.12\log \left(\frac{T_{\text{eff}}}{45000\text{\hspace{0.17em}K}}\right)\hfill \\ \hfill & +\left(0.75-1.87\log \left(\frac{T_{\text{eff}}}{45000\text{\hspace{0.17em}K}}\right)\right)\log \left(\frac{Z_{*}}{Z_{\odot }}\right).\hfill \end{array}$$(7)

This recipe is valid within the ranges 4.5 ≤ log L_*/L_⊙ ≤ 6.0, 15 ≤ M_*/L_⊙ ≤ 80, 15 000 K ≤ T_eff ≤ 50 000 K, and 0.2 ≤ Z_*/Z_⊙ ≤ 1.0. Despite its simple form, it is still a reasonably good representation of the results of the grid where most values of the mass-loss rates are represented by the recipe to within a factor of 2. More quantitatively, taking a mean across the complete grid of the ratio between the actual calculated mass-loss rate, Ṁ_mod, and that predicted by the simple fitting formula above, Ṁ_fit, yields 1.14 (median 0.91), with almost 90% of all models having a ratio lying in the range 0.5–2. Some larger factors exist, mostly originating from the bi-stability region, where, as discussed below, an overall larger scatter on the mass-loss rate is found in the grid. This then increases also the overall scatter found between the models and our fitting formula, so that when taken across the complete grid we find a standard deviation 1.67 for the ratio Ṁ_mod/Ṁ_fit. By contrast, when considering only models with T_eff ≥ 30 000 K, this standard deviation decreases significantly to 0.34 (consistent with Paper II). We emphasise, however, that this relatively large scatter does not affect a key conclusion found from our models with lower T_eff, namely that we do not find any evidence for a systematic increase in Ṁ within the bi-stability region (see the next section). As mentioned above, the models are further all computed with a standard value for the micro-turbulent velocity of 10 km s⁻¹. In Paper II we derived a dependence of Ṁ on v_turb, which can be applied to the predicted mass-loss rates here as well in case v_turb should deviate significantly from 10 km s⁻¹.

Fig. 2 Dependence of the exponent q on log Teff including a linear fit with log Teff.

5 Analysis

5.1 The bi-stability region

An important aspect of our new grid regards the behaviour within the so-called ‘bi-stability’ region around ~20 kK. Here, several important metals change their ionisation stage which affects the radiative acceleration throughout the wind. One notable transition is that of Fe IV to Fe III. Due to recombination of iron, Vink et al. (1999) find that the mass-loss rate ‘jumps’ by an order of magnitude when going from the hot to the cool side of the region at T_eff = 25 000 K (some later models seem to suggest a somewhat lower T_eff = 21 500 K; see Petrov et al. 2016). Recently, Krtička et al. (2021) also found that Ṁ increased when lowering T_eff across the bi-stability region, at least in the high-luminosity region and with more modest factors than those obtained by Vink et al. These models by Krtička et al. are, however, a bit different than the ones presented here, since they are based on a solution to the eom wherein the critical point of the flow is shifted downstream from the sonic point to the supersonic parts of the flow. As such, they also find an additional dependence on clumping for their mass-loss rates. But as discussed in detail in Sect. 2.2, it is questionable if current techniques trying to account for clumping within one-dimensional steady-state models of line-driven flows are really able to properly capture such effects.

In our models, we do not find an increase in Ṁ when crossing to the cool side of the bi-stability region. To further investigate this absence of a bi-stability jump, we computed additional models at 15 000 K (i.e. below the previously predicted bi-stability region). We selected two models from the grid, one with high and one with low luminosity, calculating the CMF transfer by fixing the velocity structure to a β-law and selecting Ṁ as predicted by Vink et al. together with the recommended υ_∞ = 1.3v_esc and β = 1 0 (we also tested both a lower and higher β, 0.6 and 3, with essentially the same result). Figure 3 shows the resulting force balance and iron ionisation balance for the low-luminosity model, including as a reference the results for the self-consistent and fully converged model with the same stellar parameters. The mass-loss rate predicted by the Vink et al. recipe is here a factor of 170 higher than the resulting self-consistent Ṁ of the converged model. For completeness, the corresponding results for the high-luminosity models are shown in Fig. A.1, where Ṁ_Vink/Ṁ = 70.

The first thing to notice from these figures is that the calculated Γ in the models using a ß-velocity law does not match the other terms in the eom (i.e. the model is not hydrodynamically consistent). Particularly important is the discrepancy around the sonic point a, where Γ(v = a) ≈ 0.1 in the β-law model (compared to the required Γ = 0.97 for a consistent model). This points to the significantly reduced radiation force in this critical region (see also the discussions in Papers I and II), which makes it impossible to drive an outflow of such high mass-loss through the sonic point using the steady-state modelling technique applied here. These results are qualitatively the same for the high-luminosity model; in this case Γ = 0 46 at the sonic point for the model with the high Ṁ predicted by Vink et al.

For these tests, it is clear that our detailed CMF-calculations do not provide enough radiation force to drive a flow with such a high mass-loss rate through the sonic point. As discussed in Papers I and II (see also Owocki & Puls 1999), the basic reason for this might be related to the dip in radiative acceleration in atmospheric layers leading up to the sonic point. The ‘force dip’ directly arises from source-function gradients in the resonance zone that are accounted for in the CMF line transfer. In order for the model to be dynamically consistent, this force-dip induces a very steep acceleration in near-sonic atmospheric regions and requires the mass-loss rate to be substantially reduced as compared to the β-law models displayed in Fig. 3.

Another important conclusion we can draw from the figures concerns the ionisation stage of iron. Although our code does not allow us to explicitly analyse the contributions to g_rad from individual elements, the alternative models by Vink et al. (2000) and Krtička et al. (2021) both point to the importance of iron for determining the overall scale of the total line force in this region. The models with very high Ṁ and a β velocity law generally have higher fractions of FeIII than the consistent models in the outer part of the wind. However, when inspecting the conditions around the sonic point, we can see that both models have iron mostly in FeIII. This means that even though the consistent models have access to the lower ionisation stages of iron, there is still not enough acceleration to launch a wind with a very high mass-loss rate. To compensate for this lack of force, the model then self-adjust towards a lower mass-loss rate in order to become dynamically consistent. This then suggests that the drastic Ṁ increase found in earlier models in this region might simply be an artefact of not being dynamically consistent around the sonic point, and not allowing properly for the feedback between radiative and velocity accelerations. The models of Vink & Sander (2021), using the setup of Müller & Vink (2008), do have dynamic solutions. However, due to their parametrised way of calculating g_rad (essentially Eq. (14) in Müller & Vink 2008) and deriving the mass-loss rate (using a global energy constraint), these models still do not capture the dynamical properties around the sonic point stemming from the CMF-force here.

The radiative acceleration predicted by FASTWIND has been thoroughly benchmarked to the alternative NLTE radiative transfer code CMFGEN (Hillier & Miller 1998), however so far mainly focusing on the O-star and early B-star regime (Puls et al. 2020). The results of the models with an effective temperature around and below the bi-stability region, have thus not yet been subjected to such an extensive calibration. The question then arises whether our atomic data at these lower temperatures include sufficient lines to represent a physically realistic result for the mass-loss rates, or if additional lines are required, which might then increase the radiative driving and result in higher mass-loss rates. To address this question, we calculated additional models from 15 kK to 22.5 kK both using FASTWIND with a beta-velocity law³ and CMFGEN (which uses a line list completely independent from that used in FASTWIND) using exactly the same stellar and wind parameters. The results for our additional models show overall a very good agreement for the radiative acceleration between both codes, especially below and around the sonic point (which is the crucial region setting the mass-loss rate). A detailed description of these models and comparisons is given in Appendix B.

Fig. 3 Wind structure for models with log(L*/L⊙) = 4.5, M* = 15 M⊙, and Teff = 15 000 K. The upper panels show Γ and the other terms in the eom throughout the wind. The lower panels show the ionisation fractions for FeIII and FeIV. The upper-left and lower-left panels are the results for the self-consistent models (with Ṁ = 2.64 × 10−10 M⊙ yr−1 and = 3028 km s−1), while the upper-right and lower-right panels are the results using a β-velocity law of β = 1.0 (with Ṁ = 4.67 × 10−8 M⊙ yr−1 and v∞ = 588 km s−1). The abscissae show a radius coordinate (lower) and the velocity (upper).

Fig. 4 Ratio of Ṁ as computed from Eq. (7) and from Paper II for the Ο stars presented in that paper that show a mean ratio of 1.12 and a standard deviation of 0.13.

5.2 Comparison to other results

Here, we analyse how the mass-loss prescription from Eq. (7) compares to other theoretical and observational results. First, we show in Fig. 4 the predictions from Eq. (7) to those of Eq. (20) in Paper II for the set of Ο stars presented in that paper (for a Galactic metallicity). The predictions agree well with each other, with the mean ratio only changing by 12%, due to slightly higher Ṁ values in the new recipe, displaying a standard deviation of 13%. When including Ο stars of LMC and SMC metallicity, the mean changes to 1.15 and the standard deviation to 0.20. The new recipe also matches slightly better with the calculated mass-loss rates of the O-star models, because it of course captures the additional (small) effects from the mass and T_eff variation also within this region.

5.2.1 Comparison to alternative theoretical mass-loss recipes

Next we compare our new prescription to other current mass-loss predictions for hot, massive stars. Figure 5 shows this comparison, using the recipe from Vink et al. (2000, 2001) as implemented in the standard ‘Dutch scheme’ in MESA, as well as from Krtička et al. (2021), for two sequences of effective temperature at a fixed high and low luminosity, respectively.

The figures show a clear reduction of the mass-loss rate as compared to the Vink et al. recipe, with slightly more pronounced differences for the low-luminosity sequence. The difference factor is about 2 to 3 on the hot side of the bi-stability region, but gets smaller when going towards higher luminosities and temperatures; for example for a hot star with T_eff > 50 kK and L_*/L_⊙ = 5.75 the rates are in good agreement (see the lower panel of Fig. 5). On the cool side of the bi-stability region the difference factor becomes on average 50, mainly because we do not predict the strong bi-stability jump that is present in this implementation of the Vink et al recipe.

In the region around 30–40 kK, there is reasonable agreement with the Krtička et al. (2021) prescription for the more luminous star, but less so when inspecting the less luminous star. This is in accordance with the results from Paper II, which show a shallower dependence of Ṁ on stellar luminosity for the results by Krtička & Kubát (2018)⁴ as compared to ours. There is also disagreement between prescriptions at lower and higher effective temperatures. At higher effective temperatures, the Krtička et al. recipe is actually out of the domain where their mass-loss rates are derived and a direct comparison with our rates therefore becomes difficult in this regime. At lower effective temperatures, their mass-loss shows a gradual increase followed by a rather steep decline. While their increase in mass-loss rate in this region is smaller than that predicted by Vink et al., the occurrence of a local maximum still means that there is a quite large disagreement with our predictions. Because of the different modelling techniques (see the discussions above), it is again difficult to judge exactly why these differences arise, but this should be further investigated in future work.

The evolution of O- and Β-type stars can vary significantly depending on whether or not this jump is included in the mass-loss recipe. Keszthelyi et al. (2017) investigate its effect on the early evolution of these stars and their rotation, because the stellar winds not only remove mass but also angular momentum from the surface layers (Langer 1998). This angular-momentum loss in turn has an effect on the evolution through altered internal mixing and transport of angular momentum (Maeder & Meynet 2000). Keszthelyi et al. (2017) find that, in order to explain the rotational velocities observed in late B supergiants, the used rates (from Vink et al. 2000, 2001) have to be reduced, either by reducing the early O-star rates or by removing the occurrence of the bi-stability jump. Both reductions could also occur simultaneously if an additional source of angular momentum loss is present or if the initial rotational velocities of O-type stars that is generally assumed (~300 km s⁻¹; e.g. Howarth et al. 1997) are too high (Simón-Díaz & Herrero 2014).

Fig. 5 Mass-loss rates predicted from Eq. (7), (black), from the Vink et al. recipe as implemented in MESA (red) and from the Krtička et al. recipe (blue), at constant log(L*/L⊙) = 4.5 and M* = 15 M⊙ for the upper panel and log(L*/L⊙) = 5.75 and M* = 50 M⊙ for the lower panel.

5.2.2 Comparison to observations

In Paper II, we presented a comparison of the theoretical predictions to a selected sample of observations for O stars, where we find good overall agreement. Considering the minor deviations of these O-star predictions with the new predictions of Eq. (7), as shown above, we still recover the same overall agreement with these observations. We can additionally look at the recent study by Hawcroft et al. (2021), who analyse a sample of Galactic O supergiants by means of a multi-diagnostics (UV+optical) genetic algorithm fitting accounting fully for potentially optically thick clumps and ‘velocity porosity’ in the line formation process (Sundqvist & Puls 2018). These empirical results again agree to within about 10% with our theoretical predictions.

With our new, more extended grid and recipe we can now also compare to empirical mass-loss rates at cooler temperatures, particularly investigating observational evidence pointing to the occurrence or absence of a bi-stability jump. In this respect, we first compare to the maximum mass-loss rates derived from radio observations such as in Benaglia et al. (2007) and Rubio-Díez et al. (2022) who also include an analysis of Herschel/PACS continuum emission data. As these mass-loss rates are obtained from radio diagnostics assuming an un-clumped outermost wind, the derived results are upper limits (Ṁ_max). Any inclusion of clumping would reduce the observationally inferred mass-loss rate by a factor of $\sqrt{f_{cl}}$, where $f_{cl}=\frac{\langle {\rho }^2\rangle }{\langle {\rho }^2\rangle }$ is the clumping factor in the radio-emitting region. While Benaglia et al. (2007) find a modestly increased wind efficiency (Ṁv_∞/L_*c) around T_eff = 21.5 kK, both studies still find a continuation of the decrease in mass-loss rates below the bi-stability jump (i.e. they do not find any evidence for a jump in mass-loss across the bi-stability region), in agreement with the overall trend of this paper. Actually, these upper-limit observed rates are even significantly lower than the rates predicted by the Vink et al. recipe. The measured upper limits are, however, also consistently above our predictions for all stars. The inclusion of clumping in the empirical mass-loss rates would then reduce the values to numbers comparable to our predictions, assuming a clumping factor of about 5–20 in the radio-emitting region. This shift could of course be altered if the clumping factor should depend on spectral type, which then would shift the different observations by different factors (see e.g. Driessen et al. 2019), which has been discussed in Markova & Puls (2008).

A comparison can also be made with analysis of stellar and wind parameters derived from the optical. Crowther et al. (2006) and Markova & Puls (2008) derive mass-loss rates from the optical spectra of B-type supergiants using CMFGEN (Hillier & Miller 1998) and FAST WIND, respectively. Crowther et al. (2006) find no significant increase in the mass-loss rate below T_eff = 24000 K even without the inclusion of clumping (which would further decrease their empirical results). Similarly Markova & Puls (2008) find that the observed mass-loss rates below the bi-stability region either decreases or increases only marginally. All observational results above thus suggest that the previously predicted bi-stability jump in mass-loss is not present (at least not for stars of not too high luminosity), but instead that rates in this region seem to follow similar scalings as for O stars (see also Keszthelyi et al. 2017).

We do, however, predict high terminal wind speeds that are not generally found in the observational literature. Across the grid, we have mean values of about 3300 km s⁻¹ and v_∞/v_esc,eff ≈ 4.5 (where ${\upsilon }_{\infty }/{\upsilon }_{\text{esc,eff}}=\sqrt{\frac{2G{M}_{*}\left(1-{\text{Γ}}_{\text{e}}\right)}{R_{*}}}$), which corresponds to a value of v_∞/v_esc ≈ 4 when using the regular escape speed, not reduced by electron scattering. This value does not significantly decrease for lower effective temperatures. These overall high terminal wind speeds may point towards either missing physics in the outer wind-regions of the models or inaccurate empirical terminal velocity estimates when the sharp edge of the UV P-Cygni lines is not clearly present in the observed spectrum. Lagae et al. (2021) computed time-dependent, line-driven instability (LDI) hydrodynamic simulations of weak-wind stars, and find that a significant portion of the weak wind is shock-heated and unable to cool down efficiently, so that the UV-line opacity becomes significantly reduced. This might then not only lead to erroneous empirical mass-loss estimates, but also to underestimated terminal speeds derived from UV-line profile fitting (see their Fig. 5 for an illustration). Additionally, as studied in Paper II, a shock-heated outer wind could also reduce the radiation force in those regions, such that if we were to account for the shock heated gas the predicted υ_∞ might be substantially reduced. On the other hand, significant portions of high-temperature gas is not really expected for winds of higher density (which can cool efficiently by radiative losses), and thus it is questionable whether this effect could also reduce terminal wind speeds in the high-luminosity parts of the HRD. Further investigations are clearly needed here, both regarding obtaining better empirical constraints on terminal wind speeds and in terms of theoretical analysis of additional physics not included in steady-state models.

5.3 Stellar evolution

Our recipe for Ṁ can readily be implemented into evolution codes. In this paper, we chose to work with MESA (Paxton et al. 2019). MESA, as mentioned above, is a modelling tool for the structure and evolution of astrophysical stellar objects and allows for simple adaptations of its input physics modules. The other_wind routine offers the option to calculate and apply a mass-loss rate from the stellar parameters at a certain time step. Here, we implemented the mass-loss recipe from Sect. 4 and applied it to all hot, hydrogen-rich stars. As mentioned before, the wind models are computed for a fixed Y_He = 0.1. As the star evolves, this value can change when, for example, material is mixed and brought to the surface. This effect is not incorporated into the mass-loss recipe, as to keep the amount of variable parameters of the grid reasonably small, and it is still unclear how much it would affect the mass-loss predictions.

To analyse the impact of our new predicted mass-loss on the evolution of stars with a mass between 15 and 80 Ṁ_⊙, we performed a differential study comparing to previous models. Naturally, many physical parameters can influence the evolution of the massive star, for example mixing processes due to overshooting from a convective region or from rotation. The purpose of the differential study here, however, is to isolate the effect of mass-loss. As a baseline for the evolution models, we chose to work with a similar setup as utilised in MIST (MESA Isochrones & Stellar Tracks Dotter 2016; Choi et al. 2016). Their stellar tracks are also computed using MESA and cover a wide range in mass and metallicity. Choi et al. (2016) present comparisons both with observational constraints and already existing models from literature. This allows them to compile a set of prescriptions and parameters applied throughout the various evolution stages (see their Table 1); in Table 1 we mention some of these that are particularly relevant for our massive-star evolution models. Concerning the mass-loss, we apply three different prescriptions depending on temperature and evolution phase. For hot (T_eff > 11 kK), hydrogen-rich phases we apply our mass-loss formula (7). This implementation is different to previous stellar-evolution calculations, such as those from MIST, where instead the mass-loss rates from Vink et al. (2000, 2001) are applied. When the star evolves to the cooler side of the HRD (T_eff < 10 kK), we switch to the empirical rates by de Jager et al. (1988). Even though our grid only extends to T_eff = 15 kK, we still apply our mass-loss formula until 10 kK, because an extrapolation of these rates to these cooler temperatures probably yields more realistic results than extrapolating the de Jager rates, which are primarily used to describe the winds of cool supergiants. For the hot (T_eff > 10 kK), stripped (with a surface hydrogen mass fraction X_surf < 0.4) stars we applied the WR-star mass-loss rates from Nugis & Lamers (2000). In the transition regions between two prescriptions, a linear interpolation between those recipes is used. This transition occurs between 10 kK and 11 kK between the cool and hot-star recipes and between X_surf = 0.4 and a few percent below for this surface abundance to transition to the stripped-star recipe.

Additionally, the table mentions parameters concerning internal transport of energy and matter. For example, exponential overshoot is applied extending the reach of the convective core and any convective layers in the envelope. This convection is described by the conventional mixing length theory (MLT; Böhm-Vitense 1958) and the location of convection is determined by the Ledoux criterion which compares the temperature gradient with the adiabatic and chemical gradients. Furthermore, to treat the radiation-dominated envelopes of massive stars, where standard MLT is not always applicable, convection is treated by artificially reducing the super-adiabaticity, effectively assuming energy is carried away by some other means than radiation or standard convection (MLT++, Paxton et al. 2013). A last prescription we mention is rotation, where the star has an initial angular rotation of Ω = 0.4Ω_crit, with Ω_crit the critical rotation where the centrifugal force equals gravity.

We used this setup to calculate tracks of stars with different initial masses, comparing to models using an equivalent setup however applying the (previously standard) Vink et al. recipe instead. Figure 6 shows the resulting HRD for a star with initial mass M_* = 60 M_* and Z_* = Z_⊙, and here below we now discuss two important aspects that are seen in these stellar tracks resulting from our new mass-loss rates.

Fig. 6 HRD for a star with an initial mass of M* = 60 M⊙ and an initial surface rotation of 350 km s−1. The two tracks are calculated assuming different mass-loss prescriptions.

5.3.1 No WR stars from standard line-driven wind stripping?

Classical WR stars are created when the hydrogen rich envelope of the massive star is somehow stripped, which in evolution models typically occurs through wind stripping or binary interactions. In the single-star models considered in this paper, only the first mechanism can be active.

In Fig. 6, we illustrate the evolution of a star with an initial mass 60 M⊙ and solar metallicity. It can be directly seen that when applying our new rates the core-helium burning phase ends without the star ever evolving into a hot WR star. On the other hand, when applying the previous higher rates, mass-loss through line-driven wind stripping is indeed enough to make the star evolve back to the blue parts of the HRD. The typical differences in mass-loss between these two characteristic models are factors of ~3 on the main sequence and ~30 for the cooler regions (where the previous rates experience the bi-stability jump, see the discussion in the previous section). The model where lower mass-loss rates are applied do not evolve to become a WR star because the outer envelope is not stripped before helium-burning ends after which the evolutionary timescale becomes too short to expel a lot more mass. This happens because mass is lost more slowly than when adopting the previous rates while the lifetime of the star is even shorter, as the star will live as a higher-mass star (compare for example 54 M_⊙ versus 45 M_⊙ in the B supergiant stage after the main sequence). So the combined effect of a shorter lifetime and reduced mass-loss makes it a lot more difficult for a star to evolve to the WR phase. For the 60 M_⊙ star, this conclusion remains unchanged, also when changing the condition when the Nugis & Lamers rates are applied. When in this model, these mass-loss rates are applied from X_surf < 0.6 instead of 0.4, more mass is lost (about 8 M_⊙ difference at the end of the evolution), but not enough to strip the star.

This differential comparison thus shows that with the implementation of our new rates in massive-star models, the formation of classical WR stars through standard, steady-state line-driven wind stripping becomes difficult, even at Galactic metallicities. Indeed, from further test models using the same setup, we find that we must increase the initial mass to >100 M_⊙ in order to form an evolved, hot WR helium star with our new mass-loss rates⁵. Because such very massive stars are also very rare, this suggests that the formation of classical WR stars might be dominated by other processes.

As mentioned above, one such alternative formation channel for classical WR stars regards envelope-stripping via binary interaction. However, observational campaigns do show signs of single WR stars, even in low-metallicity environments such as the SMC (Shenar et al. 2016, Shenar et al. 2020), so that there may a need for yet another channel by which stars strip their envelopes. One such possibility regards variable mass-loss from high-luminosity LBV-like stars that exceed the Eddington limit already in deep sub-surface layers (see e.g. Owocki 2015). Indeed, such stars are observed, but the evolution model does not pass through such a phase. As discussed in Sect. 3, LBV stars are not included in our current model grid, both because of technical difficulties with modelling deep-seated wind initiation and since the core-assumption of a steady-state outflow becomes highly questionable in this regime. As such, our recipe may well underestimate the mass-loss in this regime⁶.

5.3.2 Massive black holes

The final fate of many massive stars is to die in a SN explosion, leaving behind a neutron star or a black hole. The final mass of the black hole is largely determined by the core mass of the progenitor before the SN. Burning stages beyond carbon-burning are unstable for stars that have helium cores more massive than 30 M_⊙ (Woosley et al. 2007). The creation of an electron-positron pair requires energy at high temperatures, which reduces the amount of energy and momentum available to provide a pressure balance counteracting gravity. Contraction at high temperatures then becomes unstable, which is called the pair instability (PI). As stars get more massive, the mass of the helium core increases and results in a more violent implosion as a result of this PI. This is believed to cause stars with a He-core mass greater than 64 M_⊙ to become PI SN where a single pulse disrupts the entire star, leaving no remnant (Glatzel et al. 1985; Heger & Woosley 2002; Ober et al. 1983; Bond et al. 1984). The implosion typically happens in a series of pulses when nuclear burning is not sufficient to counter the PI; this mechanism is called pulsational pair-instability (PPI) and is believed to give rise to PPI SNe (Woosley 2017). During the pulses, the core contracts, initiates nuclear burning, and then expands and cools. The cycle repeats itself until the mass of the helium core is sufficiently lowered to avoid PPI. Due to this, the final helium core masses that manage to avoid PPI, or that have come down as a result of PPI, fall in the range 35–45 M_⊙ (Woosley 2017).

The mass of the remnant black hole is then dependent on both the He-core mass (to determine whether a PI or PPI may occur) and the mass of the remaining hydrogen envelope that gets added to the final black hole mass. In the case of the weaker winds that we predict in this paper, massive stars retain part of their hydrogen envelopes also at the end stages of their lives. Figure 6 shows the final mass after helium burning for a 60 M_⊙ Galactic star. We assume this, similarly as in Belczynski et al. (2020), to be a good representation of the final mass before core collapse. As directly seen in the figure, using our new rates we find that the total final mass is more than double than that predicted by the comparison model employing the previous standard mass-loss recipe. For the template model displayed here we find this final mass to be between 45 and 50 M_⊙ but for even higher initial masses (if the star can avoid being disrupted) this can be even higher. Indeed, Belczynski et al. (2020) found final masses as high as 70 M_⊙ from their evolution calculations using an ad hoc (factor of 5) reduction of the previously standard mass-loss rates. Here we thus recover qualitatively the same result, however with the key difference that our result does not depend on any ad hoc reduction of mass-loss rates, but rather is a natural consequence of our new theoretical predictions. We note again though, that the massive stars that would create such massive black holes indeed tend to evolve towards the LBV-like regime in the HRD, where (see discussion above) our current steady-state recipe might still not be sufficient to properly describe the mass-loss behaviour.

When considering a lower overall mass-loss, either by reducing the initial mass of the star (to below about 30 M_⊙) or by reducing the metallicity, the difference between using different mass-loss prescriptions in the prediction of the black-hole mass becomes less important. Vink et al. (2021) show, for example, that with the use of the Vink et al. rates, high-mass black holes are created at a metallicity of Z_* ≤ 0.1Z_⊙ from stars with an initial mass of 90…100 M_⊙. However, if the metallicity is not as low, the difference in the mass-loss recipe will affect how much of the envelope is blown away and thus how much mass is available to collapse into the black hole.

6 Summary

We have computed a large grid of atmosphere and wind models for hot and massive stars. The grid covers a large parameter range in stellar luminosity, mass, and effective temperature for the Galaxy, the LMC, and the SMC. From the dynamically consistent, steady-state wind models we obtain theoretical mass-loss rates and investigate their dependences on fundamental stellar parameters. The results are captured by the mass-loss prescription of Eq. (7). As in the previous O-star calculations by Björklund et al. (2021), we find strong trends of the mass-loss rate with L_* and Z_*, which are supplemented here with additional dependences on M_eff and T_eff in order to properly capture the behaviour throughout our larger grid. We additionally find that we can best represent the results of this grid when also including a temperature dependence inside the exponent for the metallicity.

We obtain generally good agreements between our predictions and observations of O stars in the Galaxy (for studies where clumping is taken into account when deriving empirical mass-loss rates). Similarly, via comparisons to observations using radio diagnostics we find that empirical upper limits of OB-star mass-loss rates are higher than our predictions by factors of around 3.5, corresponding to reasonable clumping factors of 12.

When comparing the new recipe to existing predictions of Vink et al. (2000, 2001), we again find a downward shift by approximately a factor of 3 for O stars (Björklund et al. 2021), with somewhat smaller values for the most luminous stars.

However, more notably, a much larger discrepancy is found for stars below 25 kK. Here we do not find any evidence of a bi-stability jump, which increases the mass-loss rate by an order of magnitude or more according to the mass-loss recipe by Vink et al. (2000). Additional models where we try to drive a very high mass-loss rate using a fixed β-velocity law (instead of iterating towards a consistent structure), show that there is not enough radiation force to drive such a large amount of mass through the wind sonic point. The absence of a significant mass-loss increase in the bi-stability region is also in agreement with empirical studies of Ha (Markova & Puls 2008; Crowther et al. 2006) and radio emission (Benaglia et al. 2007; Rubio-Diez et al. 2022).

Evolution models of massive stars that use the new predictions show significant differences as compared to equivalent models that employ the previous standard mass-loss rates. Here, we discuss two main influences of such decreased mass-loss rates in stellar evolution. First, with the lower mass-loss rates single stars have difficulties expelling their hydrogen-rich envelopes, which thus hinders them from turning into evolved, classical WR stars at the end of their lives. Second, since the mass of the black hole remnant depends on the final mass of its progenitor, the lower mass-loss rates potentially allow for the creation of high-mass black holes even at higher metallicities. One key uncertainty with (the mass-loss within) such massive-star evolution models, however, is that currently there is no adequate treatment of the potentially significant mass-loss from variable LBV-like stages (Smith 2014), which might ultimately help the massive star strip its hydrogen envelope and still create a WR star despite the lower steady-state mass-loss rates found here.

Acknowledgements

R.B. and J.O.S. acknowledge support from the Odysseus program of the Belgian Research Foundation Flanders (FWO) under grant G0H9218N. JOS also acknowledges support from the KU Leuven C1 grant MAESTRO C16/17/007. We thank the referee for their constructive and useful comments on the manuscript.

Appendix A Additional figures

Figure A.1 shows a similar comparison as in Fig. 3 of a self-consistent model and a model using a β-velocity law of β = 1, but now at higher mass and luminosity. The mass-loss rate of the β model is about 70 times higher than the dynamically consistent value. This model then again shows a mismatch between Γ and Λ with a value of Γ = 0.45 at the sonic points, which is significantly below the required value ≈ 1.

Fig. A.1 Same figure as Fig. 3 but now with log(L*/L⊙) = 5.75, M* = 50 M⊙, and Teff = 15 000 K. The self-consistent model has Ṁ = 7.67·0−8 M⊙ /yr and v∞ = 3111 km/s, while the β model has Ṁ = 5.38 · 10−6 M⊙ /yr and v∞ = 451 km/s.

Appendix B Comparison with Γ from CMFGEN

The self-consistent models presented and discussed in the main part of this paper show significant differences in their wind-structure (both with respect to mass-loss rate and velocity law) compared to previous predictions (e.g., from Vink et al. 2000, 2001), particularly in the B-star range. To underpin these results, we need to convince ourselves that the radiative acceleration as calculated by FASTWIND does not suffer from specific problems, especially from deficiencies in the used line-lists resulting in a predicted line-acceleration that is too low. To this end, we proceed as in Puls et al. (2020), who compared the radiative acceleration in the O/early B-star regime (down to T_eff ≈ 25 kK) with corresponding results from CMFGEN models. The CMFGEN code has also been used in, for example, theoretical investigations about the location of the potential bi-stability region preformed by Petrov et al. (2016).

We set up a small model grid enclosing the critical region around T_eff ≈ 20 kK and concentrating on B supergiants and hypergiants, assuming a prototypical β = 1 velocity law as also used in the calculations by Vink et al. (2000, 2001). We neglected wind-clumping and X-ray emission from wind-embedded shocks. The division between the β law and the quasi-hydrostatic stratification was set (in most cases) to 0.5 v_sound. For each model we considered two mass-loss rates, a high one following the trend predicted with the mass-loss recipes by Vink et al. (2000, 2001) and a lower one to check the influence of a less dense wind. All parameters of the grid-stars are provided in Table B.1.

The outcome of our modelling is displayed in Fig. B.1, where we compare the total radiative acceleration as calculated by FASTWIND v11 and CMFGEN as a function of velocity, measured in units of gravitational acceleration (i.e., we display Γ as defined in Eq. (1)). The red dashed-dotted lines indicate the conventional Gamma factor for electron scattering, Γ_e, as defined in the beginning of Sect. 4, also as a function of velocity.

The comparison shows clearly that in most cases the photo-spheric accelerations (even in those cases where Γ is close to unity and thus has a significant effect on the structure) agree quite well, and that the predicted wind acceleration (at least until 0.5v_∞) does not show significant differences (in all cases, below 20%). Only for the models at 20 and 17.5 kK, FASTWIND predicts a larger acceleration in the outer wind, by a factor 1.5 to 2. Though the origin of this discrepancy is currently unclear (work in progress), our modelling shows that, at least compared to CMFGEN, there is no deficiency in the line force that might have influenced our results from the previous sections. Particularly in the most important region around the sonic point (where M is determined) until a few hundreds km s⁻¹ (a region that also controls the lowermost wind-regime, due to backscattering), the agreement between the two independent codes and predictions is reasonable (often even better).

Given that the two codes use a rather different computational approach and different atomic databases, the outcome of our comparison is ensuring, and we are confident that our calculations presented in the main part of this paper are not affected by too low g_rad because of missing lines. As already noted in previous publications of this series, and also by Puls et al. (2020) for the O-star case, also here Γ_tot = 1 is reached only at substantial velocities (100 km s⁻¹or more), whereas, in self-consistent wind models, this should happen very close to the sonic point (here: 15 … 18 km s⁻¹, in dependence of T_eff). This highlights that the displayed models are far away from being hydrodynamically self-consistent. To achieve such a self-consistency, and due to the basic dependence g_rad ∝ (dv/dr)/ρ (e.g., Castor et al. 1975), only a steeper velocity law and a lower mass-loss rate will lead to a higher g_rad in the transonic region, which is also the basic outcome of our self-consistent calculations.

Fig. B.1 ${\text{Γ}}_{\text{tot}}=g_{\text{rad}}^{\text{tot}}/g_{\text{grav}}$ for all models from our grid (see Table B.1), as a function of velocity. Black represents results from FASTWIND v11 and green results from CMFGEN. To improve the visibility, all Γ values have been multiplied with factors 10i, i ∈ [0,4], from bottom to top. The dashed-dotted red lines correspond to Γe(r) for pure electron scattering (its local variation relates to changes in the ionisation structure, mostly of H and He), and the dashed lines indicate, for each model, the relation Γtot = 1. We note that all models displayed here have been calculated with a fixed ß = 1 velocity law, i.e. they are not self-consistent. See text.

References

All Tables

All Figures

	Fig. 1 HRD showing the coverage of the grid of models. Each point in the L*-Teff plane represents several models with varying mass, shown by the overlaid squares. Also shown are the stellar tracks including an initial rotation of 0.4 times the critical rotation for 15, 20, 30, 40 and 60 M⊙ stars computed using MESA.
In the text

	Fig. 2 Dependence of the exponent q on log Teff including a linear fit with log Teff.
In the text

Fig. 3 Wind structure for models with log(L*/L⊙) = 4.5, M* = 15 M⊙, and Teff = 15 000 K. The upper panels show Γ and the other terms in the eom throughout the wind. The lower panels show the ionisation fractions for FeIII and FeIV. The upper-left and lower-left panels are the results for the self-consistent models (with Ṁ = 2.64 × 10−10 M⊙ yr−1 and = 3028 km s−1), while the upper-right and lower-right panels are the results using a β-velocity law of β = 1.0 (with Ṁ = 4.67 × 10−8 M⊙ yr−1 and v∞ = 588 km s−1). The abscissae show a radius coordinate (lower) and the velocity (upper).

In the text

	Fig. 4 Ratio of Ṁ as computed from Eq. (7) and from Paper II for the Ο stars presented in that paper that show a mean ratio of 1.12 and a standard deviation of 0.13.
In the text

	Fig. 5 Mass-loss rates predicted from Eq. (7), (black), from the Vink et al. recipe as implemented in MESA (red) and from the Krtička et al. recipe (blue), at constant log(L*/L⊙) = 4.5 and M* = 15 M⊙ for the upper panel and log(L*/L⊙) = 5.75 and M* = 50 M⊙ for the lower panel.
In the text

	Fig. 6 HRD for a star with an initial mass of M* = 60 M⊙ and an initial surface rotation of 350 km s−1. The two tracks are calculated assuming different mass-loss prescriptions.
In the text

	Fig. A.1 Same figure as Fig. 3 but now with log(L*/L⊙) = 5.75, M* = 50 M⊙, and Teff = 15 000 K. The self-consistent model has Ṁ = 7.67·0−8 M⊙ /yr and v∞ = 3111 km/s, while the β model has Ṁ = 5.38 · 10−6 M⊙ /yr and v∞ = 451 km/s.
In the text

Fig. B.1 ${\text{Γ}}_{\text{tot}}=g_{\text{rad}}^{\text{tot}}/g_{\text{grav}}$ for all models from our grid (see Table B.1), as a function of velocity. Black represents results from FASTWIND v11 and green results from CMFGEN. To improve the visibility, all Γ values have been multiplied with factors 10i, i ∈ [0,4], from bottom to top. The dashed-dotted red lines correspond to Γe(r) for pure electron scattering (its local variation relates to changes in the ionisation structure, mostly of H and He), and the dashed lines indicate, for each model, the relation Γtot = 1. We note that all models displayed here have been calculated with a fixed ß = 1 velocity law, i.e. they are not self-consistent. See text.

In the text