New predictions for radiation-driven, steady-state mass-loss and wind-momentum from hot, massive stars III. Updated mass-loss rates for stellar evolution

Massive stars lose a large fraction of their mass to radiation-driven winds throughout their entire life. These outflows impact both the life and death of these stars and their surroundings. Theoretical mass-loss rates of hot, massive stars are derived to be used in applications such as stellar evolution. The behaviour of these rates in the OB-star regime is analysed, and their effects on massive-star evolution predictions is studied. Dynamically-consistent models are computed by solving the spherically symmetric, steady-state equation-of-motion for a large grid of hot, massive stars. The radiative acceleration is derived from non-local thermodynamic equilibrium radiative transfer in the co-moving frame. The resulting mass-loss rates are used to derive a simple scaling recipe with stellar parameters and to evaluate some first impacts upon evolution tracks. We provide a new prescription for steady-state, radiation-driven mass-loss from hot, massive stars depending on their fundamental parameters. The rates for O-stars are lower by about a factor ~3 than the rates typically used in previous stellar-evolution calculations, where differences generally decrease with increasing luminosity and temperature. For cooler B-giants/supergiants we find larger discrepancies. This arises because we do not find any systematic increase in mass-loss rates below the so-called bi-stability region; indeed, our results do not show any sign of a significant bi-stability jump within the parameter range covered by the grid. Due to the lower mass-loss rates we find that envelopes are not easily stripped by means of standard steady-state winds, making it difficult to create classical Wolf-Rayet stars via this channel. However, a remaining key uncertainty regarding these predictions concerns unsteady mass loss for very high-luminosity stars close to the Eddington limit as well as the impact of non-line-driven winds.


Introduction
The lives of massive stars are largely controlled by their strong radiation-driven winds. The resulting mass loss has an impact on the stellar luminosity, its lifetime, apparent temperature, etc., 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 (Fierlinger et al. 2016). Additionally, the amount of mass loss largely determines which type of supernova (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 (GW) which arise from the merging of a binary black holes. The GW de-tections by the Advanced-LIGO interferometers provided this first observational evidence of black hole mergers (Abbott et al. 2016a) 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 of quite weak stellar winds. Indeed, mass loss is highlighted as a main uncertainty in massive-star progenitor models of such gravitational-wave sources (Abbott et al. 2016b).
Accurate predictions of the mass loss of massive stars are thus necessary to predict 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) 1 ) through simple scaling recipes depending on fundamental stellar parameters. The prescriptions used in most massive-star evolution models are typically derived from either theoretical models or empirical studies focusing on certain parts of the Hertzprung-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 the series we developed new atmosphere and wind models in order to provide mass-loss predictions and predicted wind-structure , from now on Paper I) and applied them to a sample of O-stars (Björklund et al. 2021, from now on Paper II). The wind models solve the spherically symmetric, steady-state equation-of-motion, relying on non-local thermodynamic (NLTE) radiative transfer in the co-moving frame 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(Vink et al. , 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. 2021).
In this paper we further study the dependence of the massloss rate on important stellar parameters, including now also 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 Ostar 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. 2000Vink et al. , 2001Driessen 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 stars and high-mass black holes.

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 equation of motion (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 (non-local thermodynamic equilibrium) radiative transfer of fastwind using a co-moving frame (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) ill-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Ṁ ∝ 1 g 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 eqns. 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(τ 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 ∼ 100R * , where the winds have almost reached their terminal wind speed ∞ ≡ (r = R max ).

Convergence
The radiative transfer and the equation of motion 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 Paper 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), (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, ionization balance and NLTE occupation numbers (the latter typically to within a maximum error of a few percent). In our hydrodynamic models final convergence is defined through the maximum error in the equation of motion (see Papers I and II). This error compares Γ = g rad /g to the other terms in the eom according to where Our models are formally defined as converged when the maximum error radially in f err 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 allow for somewhat larger deviations (see discussion below). Additionally, we require both the mass-loss rate and the velocity structure not to vary more than 2% and 3% respectively between the last two hydrodynamic steps.

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 discussion below) with an onset beyond the sonic point typically does not significantly affect the predicted massloss rates, because they are mostly sensitive to the conditions locally around the sonic point. This is different than, for example, models deriving a mass-loss rate only from a global energy balance constraint and not from a locally consistent solution to the equation of motion (Muijres et al. 2011). Also the recent simulations by Krtička et al. (2021) 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 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, let us also point out here that current inclusions of wind clumping into 1D line-driven steady-state hydrodynamic simulations (such as those presented here) are of a quite adhoc 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 ionization balance and thus also have an implicit effect on g rad . But in a steadystate 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 equation of motion. Indeed, to our knowl-edge all non-parametrised 2 steady-state line-driven wind models attempting to include 'clumping', have computed the ionization for the clumps themselves, but then (inconsistently) solved the equation of motion 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 if/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 discussion in Sect. 5.2.2).

Grid setup
We set up a grid of radiation-driven wind models covering the hot region of the HR-diagram where massive stars evolve most of their lives, excluding only the cool temperatures (<15000 K), where we find, for example, red supergiants which presumably are driven by other processes (e.g. Kee et al. 2021), and the extremely hot temperatures of stripped Wolf-Rayet stars 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. 2011Paxton et al. , 2013Paxton et al. , 2015Paxton et al. , 2018Paxton et al. , 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 6 L for L * , 15000 − 50000 K for T eff , and 15 − 80 M for M * . Fig. 1 shows the extent of the model grid on the Hertzprung-Russel diagram. 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.5Z and 0.2Z corresponding to those of the Large and Small Magellanic Cloud (LMC and SMC) respectively. This leads to a total of 516 wind models. This setup of the model grid allows us to investigate the dependencies on the different stellar parameters such as luminosity and metallicity, mass, and effective temperature. These latter 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 He = N He N H = 0.1 with also the micro-turbulent velocity kept fixed at turb = 10 km/s.

Regions of the HR-diagram
The grid is chosen such that we can compute steady-state wind models that experience radiation-driven outflows. Hot, mainsequence (and post-main sequence super-giants 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 wellbehaved and converge in a relatively straightforward way. There are, however, some areas of the HR-diagram 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, maximum 2 % error for the equation-of-motion over the complete radial grid, but rather show hydrodynamic solutions displaying somewhat larger maximum errors. It is important to realize, 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 equation-ofmotion in these types of steady-state line-driven wind models (see also 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 max err < 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 (actually, our 2% basecriterion is very stringent, e.g. 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 which do not reach this strict criterion then, we calculate the final mass-loss rate as an average of allṀ in the hydrodynamic iteration steps with f max err <10%, i.e. we compute Ṁ <10% . To justify this averaging-approach, we compare 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 |Ṁ <2% − Ṁ <10% Ṁ <2% | 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 well modelled with our current technique, both because our atomic data base 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 Wolf Rayet stars (WNh, sometimes called 'slash stars' or 'O-stars on steroids') and luminous blue variables (LBV).

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 a form: that is where M eff = M * (1 − Γ e ) is the effective stellar mass, reduced by the electron-scattering Eddington parameter Γ e = κ e L * 4πcGM * applying a constant κ e in the model-fit. The value for κ e is calculated by 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 N He = 0.1, this gives the typical value of 0.34 cm 2 /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 of 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 Ostar 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 dependencies on other parameters and the limited number of metallicity points used in our model grid. Finally then, we construct a re- lation for the mass-loss rate from a multi-linear fit including the Z * (T eff ) dependence as follows logṀ = − 5.52 This recipe is valid within the ranges 4.5 ≤ log L * /L ≤ 6.0, 15 ≤ M * /L ≤ 80, 15000K ≤ T eff ≤ 50000K, 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 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 to 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 ra-tioṀ 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 emphasize, 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 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 turb , which can be applied to the predicted mass-loss rates here as well in case turb should deviate significantly from 10 km/s.

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 =25000 K (some later models seem to suggest a somewhat lower T eff =21500 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. Also these models by Krticka et al., however, are a bit different than the ones presented here, since they are based on a solution to the equation-of-motion 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 1D steady-state models of line-driven flows are really able to properly capture such effects. In our models, we do not find an increase ofṀ when crossing to the cool side of the bi-stability region. To further investigate this absence of a bi-stability jump, we compute additional models at 15000 K (i.e. below the previously predicted bi-stability region). We select 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.3 esc and a β = 1.0 (we have also tested both a lower and higher β, 0.6 and 3, with essentially the same result). Fig. 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 170 higher than the resulting self-consistentṀ of the converged model. For completeness, the corresponding results for the highluminosity 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 Γ( = 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 discussions in Paper I and II), which makes it impossible to drive an outflow of such a high mass loss through the sonic point using the steady-state modelling technique applied here. These results are qualitatively the same for the highluminosity 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 Paper 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 which 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 drastiċ M 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 , using the set-up 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 massloss 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 3 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.

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 equation (7) to those of Eq.
(20) in Paper II for the set of O-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 higheṙ M values in the new recipe, displaying a standard deviation of 13%. When including O-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 massloss 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.

Comparison to alternative theoretical mass-loss recipes
Next we compare our new prescription to other current mass-loss predictions for hot, massive stars. Fig. 5 shows this comparison, using the recipe from Vink et al. (2000Vink et al. ( , 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 lower panel of Fig. 5). On the cool side of the bi-stability region the differencefactor 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 to 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) 4 as compared to ours. There is also disagreement between prescriptions at lower and higher effective temperatures. At higher effective temperatures, the Krticka 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 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 B type stars can vary significantly depending on whether or not this jump is included in the massloss 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. 2000Vink et al. , 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).

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/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. (2021) who also include an analysis of Herschel/PACS continuum emission data. As these mass-loss rates are obtained from radio-diagnostics assuming an unclumped 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 f cl , where f cl = ρ 2 ρ 2 is the clumping factor in the radio-emitting region. While Benaglia et al. (2007) find a modestly increased wind efficiency (Ṁ ∞ /L * c) around T eff =21.5 kK, both studies still find a continuation of the decrease of 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 fastwind, respectively. Crowther et al. (2006) find no significant increase in the mass-loss rate below T eff = 24000K even without the inclusion of clumping (which would further de-crease 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 which are not generally found in the observational literature. Across the grid, we have mean values of about 3300 km/s and ∞ / esc,eff ≈ , which corresponds to a value of ∞ / 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 ultraviolet 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 hightemperature 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 steadystate models.

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 in 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 M , we perform 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  Choi et al. (2016) for an extended list of the prescriptions.

Ingredient
Adopted physics Details Mass loss This paper T eff > 11kK, X surf > 0.4 de Jager et al.
T eff < 10kK Nugis & Lamers T eff > 10kK, X surf < 0.4 Overshoot f ov,core = 0.016 exponential overshoot f ov,envelope = 0.0174 Convection Ledoux + MLT Includes superadiabatic convection (MLT++) Rotation init = 350 km/s Ω/Ω crit = 0.4 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 Tab. 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 > 11kK), hydrogen-rich phases we apply our mass-loss formula (7). This implementation is different to previous stellar evolution calculation, such as those from MIST, where instead the mass-loss rates from Vink et al. (2000Vink et al. ( , 2001 are applied. When the star evolves to the cooler side of the HR-diagram (T eff < 10kK), we switch to the empirical rates by de Jager et al. (1988). Even though our grid only extends to T eff = 15kK, 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 > 10kK), stripped (with a surface hydrogen mass fraction X surf < 0.4) stars we apply the WR-star mass loss 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 superadiabaticity, 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 have used this setup to calculate tracks of stars with different initial masses, comparing to models using an equiva-  lent setup however applying the (previously standard) Vink et al. recipe instead. Fig. 6 shows the resulting Hertzprung-Russel (HR) diagram 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 massloss rates.

No WR stars from standard line-driven wind stripping?
Classical Wolf-Rayet (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 discussion in 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 time-scale 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 metalicities. Indeed, from further test-models using the same set-up, we find that we must increase the initial mass to > 100M in order to form an evolved, hot WR helium star with our new mass-loss rates 5 . 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 coreassumption of a steady-state outflow becomes highly questionable in this regime. As such, our recipe may well underestimate the mass loss in this regime 6 .

Massive black holes
The final fate of many massive stars is to die in a supernova (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 that 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 Hecore mass greater than 64 M to become PI supernova 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 supernovae 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 the PPI. Due to this, the final helium core masses that manage to avoid PPI, or that have come down as a result of PPI falls in the range 35 − 45M (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. Fig.  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 massloss recipe. For the template model displayed here we find this final mass to be between 45-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.  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.

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 and the Large and Small Magellanic Clouds. From the dynamically-consistent, steady-state wind models we obtain theoretical mass-loss rates and investigate their dependencies 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 dependencies on M eff and T eff in order to properly capture the behaviour throughout our larger grid. We additionally found that we could best represent the results of this grid when including also 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, when comparing to observations using radio-diagnostics we find that empirical upper limits of OBstar mass-loss rates are higher than our predictions with factors around 3.5, corresponding to reasonable clumping factors of 12.
When comparing the new recipe to existing predictions of Vink et al. (2000Vink et al. ( , 2001, we again find a downwards shift by ap-proximately 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 also agrees with empirical studies of Hα (Markova & Puls 2008;Crowther et al. 2006) and radio emission (Benaglia et al. 2007;Rubio-Díez et al. 2021).
Evolution models of massive stars using the new predictions show significant differences as compared to equivalent models employing 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 to turn into evolved, classical Wolf-Rayet 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 to strip its hydrogen-envelope and still create a WR star despite the lower steady-state mass-loss rates found here.   Fig. 3 of a selfconsistent 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.

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 windstructure (both with respect to mass-loss rate and velocity law) compared to previous predictions (e.g., from Vink et al. 2000Vink et al. , 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 latter code has been used, e.g., also in 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 super-/hypergiants, assuming a prototypical β = 1 velocity law as also used in the calculations by Vink et al. (2000Vink et al. ( , 2001. We neglected windclumping 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 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. (2000Vink et al. ( , 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 photospheric accelerations (even in those cases where Γ is close to unity and thus has a significant effect on the structure) agree quite well, and that also the predicted wind acceleration (at least until 0.5 ∞ ) 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 which might have influenced our results from the previous sections. Particularly in the most important region around the sonic point (wherė M is determined) until few hundreds of km s −1 (a region which 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 data bases, the outcome of our comparison is ensuring, and we are confident that our calcula-tions 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 −1 or more), whereas, in self-consistent windmodels, this should happen very close to the sonic point (here: 15 . . . 18 km s −1 , 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 ∝ (d /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.