11institutetext: KU Leuven, Instituut voor Sterrenkunde, Celestijnenlaan 200D,
3001 Leuven, Belgium
11email<EMAIL_ADDRESS>22institutetext: LMU München,
Universitätssternwarte, Scheinerstr. 1, 81679 München, Germany33institutetext:
Centro de Astrobiologia, Instituto Nacional de Tecnica Aerospacial, 28850
Torrejon de Ardoz, Madrid, Spain

# New predictions for radiation-driven, steady-state mass-loss and wind-
momentum from hot, massive stars

III. Updated mass-loss rates for stellar evolution

R. Björklund 11 J.O. Sundqvist 11 S.M. Singh 11 J. Puls 22 F. Najarro 33

(Received 3 August 2021; accepted 10 January 2022)

###### Abstract

Context. 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.

Aims. 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 studied.

Methods. Dynamically-consistent models are computed by solving the spherically
symmetric, steady-state equation-of-motion for a large grid of hot, massive
stars with different metallicities. The radiative acceleration is derived from
non-local thermodynamic equilibrium radiative transfer in the co-moving frame,
and all models cover a large spatial range from deep sub-sonic atmospheric
layers into the radiation-driven and highly supersonic wind outflow. The
resulting mass-loss rates are used to derive a simple scaling recipe with
stellar parameters (luminosity, mass, effective temperature, metallicity), and
the new recipe is used to evaluate some first impacts upon massive-star
evolution tracks.

Results. We provide a new prescription for steady-state, radiation-driven
mass-loss from hot, massive stars depending on their fundamental parameters.
In accordance with our previous work, the rates for O-stars are lower by about
a factor $\sim$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,
up to one or even two orders of magnitude. 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 ($T_{\rm eff}\geq 15$ kK).
Due to the lower mass-loss rates we find that massive-star envelopes are not
easily stripped by means of standard steady-state winds, making it difficult
to create classical Wolf-Rayet stars via this channel. Moreover, since the
stars retain more mass right before they die as supernovae, our new rates make
it possible to create black holes of higher masses than in previous models,
even at Galactic metallicity. 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.

###### Key Words.:

Stars: atmospheres – Stars: early-type – Stars: massive – Stars: mass-loss –
Stars: winds, outflows – Stars: evolution

##  1 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 detections 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)111These references only refer to the latest
version of the code, the original work for these appear in older references.)
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
(Sundqvist et al., 2019, 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
$\sim 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., 2021).

In this paper we further study the dependence of the mass-loss 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 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 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 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_{\text{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 $\beta$-velocity
law, where the steepness of the acceleration towards the terminal wind speed
$\varv_{\infty}$ is set by the exponent $\beta$, which is stitched to a quasi-
hydrostatic atmosphere in the inner region (see Eq. 3 of Paper II). As for
$\dot{M}$, 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
$\dot{M}\propto\frac{1}{g_{\text{rad}}}$ (see Paper I). Then from $\varv(r)$
and $\dot{M}$ 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_{\ast}$ and mass
$M_{\ast}$ are kept fixed. The radius is updated such that the final converged
model has $R_{\ast}\equiv r(\tau_{F}=2/3)$, where $R_{\ast}$ is the input
stellar radius and $\tau_{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_{\rm max}\sim 100R_{\ast}$, where the winds have almost reached
their terminal wind speed $\varv_{\infty}\equiv\varv(r=R_{\rm max})$.

###  2.1 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_{\rm 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 ($\rho(r),\varv(r),T(r)$) is then used to
again converge the radiative acceleration ($g_{\rm 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
$\Gamma=g_{\text{rad}}/g$ to the other terms in the eom according to

$\displaystyle f_{\rm err}(r)=1-\frac{\Lambda(r)}{\Gamma(r)},$ (1)

where

$\displaystyle\Lambda=\frac{1}{g}\left(\varv\frac{d\varv}{dr}\left(1-\frac{a^{2}}{\varv^{2}}\right)+g-\frac{2a^{2}}{r}+\frac{da^{2}}{dr}\right).$
(2)

Our models are formally defined as converged when the maximum error radially
in $f_{\rm err}$

$\displaystyle f_{\rm err}^{\text{max}}=\text{max}(\text{abs}(f_{\rm err}))$
(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 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.

###  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 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
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 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 ionization balance and thus also have an implicit
effect on $g_{\rm rad}$. But in a steady-state wind model, $g_{\rm 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_{\rm rad}\sim 1/\rho^{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_{\rm 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_{\rm rad}$ and $\dot{M}$, on the expensive of
having to assume a spherically symmetric and smooth steady-state equation of
motion. Indeed, to our knowledge all non-parametrised222e.g. the attempts by
Sundqvist et al. (2014) to account for porosity in velocity space when
computing the line-force still used a parametrised approach. 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_{\rm 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).

##  3 The new grid

###  3.1 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 (Sander et al.,
2020; Poniatowski et al., 2021). For a given metallicity $Z_{\ast}$, 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_{\ast}$, $T_{\rm eff}$, and $M_{\ast}$ 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^{6}$ L⊙ for $L_{\ast}$, $15000-50000$ K for $T_{\rm eff}$, and
$15-80$ M⊙ for $M_{\ast}$. 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.

Figure 1: HR-diagram showing the coverage of the grid of models. Each point in
the $L_{\ast}-T_{\rm eff}$ plane represents several models with varying mass,
shown by overlaying squares in the figure. 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.

Each model is calculated with a Solar metallicity of $Z_{\odot}=0.014$ and
lower metallicities of $0.5Z_{\odot}$ and $0.2Z_{\odot}$ 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_{\rm He}=\frac{N_{\rm He}}{N_{\rm
H}}=0.1$ with also the micro-turbulent velocity kept fixed at
$\varv_{\text{turb}}$ = 10 km/s.

###  3.2 Regions of the HR-diagram

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 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 well-behaved 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-of-motion 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_{\rm
err}^{\text{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 (actually, our 2% base-criterion 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
$\dot{M}$ in the hydrodynamic iteration steps with $f_{\rm
err}^{\text{max}}$¡10%, i.e. we compute $\langle\dot{M}_{<10\%}\rangle$. To
justify this averaging-approach, we compare the predicted mass-loss rate of
all formally converged models, i.e. those which reach $\dot{M}_{<2\%}$, to the
$\dot{M}$ that would have been derived if we instead had applied the
averaging-method $\langle\dot{M}_{<10\%}\rangle$ for also these model-stars.
From this, we find that the difference
$|\frac{\dot{M}_{<2\%}-\langle\dot{M}_{\rm<10\%}\rangle}{\dot{M}_{<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 $\sim$20%. This
comparison supports the use of $\langle\dot{M}_{\rm<10\%}\rangle$ as the mass-
loss rate of models that do not formally reach $\dot{M}_{<2\%}$.

We limit the parameter-range of the grid to a maximum luminosity $\log
L_{\ast}/\text{L}_{\odot}=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).

##  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 $\dot{M}$
that depends on the stellar luminosity, mass, effective temperature and
metallicity. For this we assume a form:

$\displaystyle\dot{M}\propto L_{\ast}^{m}M_{\rm eff}^{n}T_{\rm
eff}^{p}Z_{\ast}^{q},$ (4)

that is

$\displaystyle\log\dot{M}=C+m\log L_{\ast}+n\log M_{\rm eff}+p\log T_{\rm
eff}+q\log Z_{\ast},$ (5)

where $M_{\rm eff}=M_{\ast}(1-\Gamma_{e})$ is the effective stellar mass,
reduced by the electron-scattering Eddington parameter
$\Gamma_{e}=\frac{\kappa_{e}L_{\ast}}{4\pi cGM_{\ast}}$ applying a constant
$\kappa_{e}$ in the model-fit. The value for $\kappa_{e}$ is calculated by

$\displaystyle\kappa_{e}=\frac{0.4\left(1+i_{\rm He}Y_{\rm
He}\right)}{1+4Y_{\rm He}},$ (6)

with $i_{\rm He}$ the ionisation stage of helium (with $i_{\rm He}=0$ for
neutral Helium). For O-stars, with doubly-ionised helium and our assumed
$N_{\rm He}=0.1$, this gives the typical value of $0.34\rm\ 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_{\ast}$ 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_{\rm eff}$ instead of $\log L_{\ast}$. In the previous grid
of O-stars, this was however essentially the same because of the strong
relationship of $L_{\ast}$ and $T_{\rm 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 $\dot{M}$ only changes by at most $\sim 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_{\rm 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.

Figure 2: Dependence of the exponent $q$ on $\log T_{\rm eff}$ including a
linear fit with $\log T_{\rm eff}$.

Finally then, we construct a relation for the mass-loss rate from a multi-
linear fit including the $Z_{\ast}(T_{\rm eff})$ dependence as follows

$\displaystyle\log\dot{M}=$ $\displaystyle-5.52$
$\displaystyle+2.39\log\left(\frac{L_{\ast}}{10^{6}\text{L}_{\odot}}\right)$
$\displaystyle-1.48\log\left(\frac{M_{\rm eff}}{45\text{M}_{\odot}}\right)$
(7) $\displaystyle+2.12\log\left(\frac{T_{\rm eff}}{45000K}\right)$
$\displaystyle+\left(0.75-1.87\log\left(\frac{T_{\rm
eff}}{45000K}\right)\right)\log\left(\frac{Z_{\ast}}{Z_{\odot}}\right).$

This recipe is valid within the ranges $4.5\leq\log
L_{\ast}/\text{L}_{\odot}\leq 6.0$, $15\leq M_{\ast}/\text{L}_{\odot}\leq 80$,
$15000{\rm K}\leq T_{\rm eff}\leq 50000{\rm K}$, and $0.2\leq
Z_{\ast}/Z_{\odot}\leq 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, $\dot{M}_{\rm mod}$, and that predicted
by the simple fitting formula above, $\dot{M}_{\rm 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 ratio $\dot{M}_{\rm
mod}/\dot{M}_{\rm fit}$. By contrast, when considering only models with
$T_{\rm eff}\geq 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_{\rm eff}$, namely that we do not find any evidence for a
systematic increase in $\dot{M}$ 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 $\dot{M}$ on $\varv_{\text{turb}}$, which can be
applied to the predicted mass-loss rates here as well in case
$\varv_{\text{turb}}$ should deviate significantly from 10 km/s.

##  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 $\sim$ 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_{\rm eff}$=25000 K (some later models seem to suggest a somewhat
lower $T_{\rm eff}$=21500 K, see Petrov et al., 2016). Recently, Krtička et
al. (2021) also found that $\dot{M}$ increased when lowering $T_{\rm 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 $\dot{M}$ 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 $\beta$-law and selecting $\dot{M}$ as
predicted by Vink et al. together with the recommended
$\varv_{\infty}=1.3\varv_{\rm esc}$ and a $\beta=1.0$ (we have also tested
both a lower and higher $\beta$, 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 $\dot{M}$ of the converged model. For
completeness, the corresponding results for the high-luminosity models are
shown in Fig. 7, where $\dot{M}_{\rm Vink}/\dot{M}=70$.

Figure 3: Wind structure for models with
$\log(L_{\ast}/\text{L}_{\odot})=4.5$, $M_{\ast}=15\text{M}_{\odot}$, and
$T_{\rm eff}=15000$K. The upper panels show $\Gamma$ and the other terms in
the eom throughout the wind. The lower panels show the ionisation fractions
for FeIII and FeIV. The panels on the upper left and lower left are the
results for the self-consistent models (with $\dot{M}=2.64\cdot 10^{-10}$
M⊙/yr and $\varv_{\infty}=3028$ km/s), while the panels on the upper right and
lower right are the results using a $\beta$-velocity law of $\beta=1.0$ (with
$\dot{M}=4.67\cdot 10^{-8}$ M⊙/yr and $\varv_{\infty}=588$ km/s). The abscissa
show a radius-coordinate below and the corresponding velocity on top.

The first thing to notice from these figures is that the calculated $\Gamma$
in the models using a $\beta$-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
$\Gamma(\varv=a)\approx 0.1$ in the $\beta$-law model (compared to the
required $\Gamma=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 high-luminosity model; in this case $\Gamma=0.46$ at the sonic point for
the model with the high $\dot{M}$ 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 $\beta$-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_{\rm 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 $\dot{M}$ and a
$\beta$-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
$\dot{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 Vink & Sander (2021), using the set-up of Müller
& Vink (2008), do have dynamic solutions. However, due to their parametrised
way of calculating $g_{\text{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 law333Since cmfgen does not solve the hydrodynamical equation we make
use of $\beta$-law models for the purpose of a one to one comparison. This
should produce almost identical wind and density structures, except for a
subtle difference in the definition of the $\beta$-law and how the connection
photosphere-wind is performed in each code. 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.

###  5.2 Comparison to other results

Here we analyse how the mass-loss prescription from Eq. (4) compares to other
theoretical and observational results. First, we show in Fig. 4 the
predictions from equation (4) to those of Eq. (20) in Paper II for the set of
O-stars presented in that paper (for a Galactic metallicity).

Figure 4: Ratio of $\dot{M}$ as computed from Eq. (4) and from Paper II for
the O-stars presented in that paper which show a mean ratio of 1.12 and a
standard deviation of 0.13.

The predictions agree well with each other, with the mean ratio only changing
by 12%, due to slightly higher $\dot{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 mass-loss rates of
the O-star models, because it of course captures the additional (small)
effects from the mass and $T_{\rm 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. Fig. 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.

Figure 5: Mass-loss rates predicted from Eq. (4) (black), from the Vink et al.
recipe as implemented in MESA (red), and from the Krticka et al. recipe
(blue), at constant $\log(L_{\ast}/\text{L}_{\odot})=4.5$ and
$M_{\ast}=15\text{M}_{\odot}$ for the upper panel and
$\log(L_{\ast}/\text{L}_{\odot})=5.75$ and $M_{\ast}=50\text{M}_{\odot}$ for
the lower panel.

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_{\rm
eff}$¿50 kK and $L_{\ast}/L_{\odot}=5.75$ the rates are in good agreement (see
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 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 $\dot{M}$ on stellar
luminosity for the results by Krtička & Kubát (2018)444The results from
Krtička & Kubát (2018) for O-stars have been incorporated in the more general
recipe in Krtička et al. (2021) for O and B stars. 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 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
($\sim 300$ km/s, e.g. Howarth et al., 1997) are too high (Simón-Díaz &
Herrero, 2014).

####  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. (4), 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 ($\dot{M}_{\rm
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\rangle^{2}}$ is the clumping
factor in the radio-emitting region. While Benaglia et al. (2007) find a
modestly increased wind efficiency ($\dot{M}\varv_{\infty}/L_{\ast}c$) around
$T_{\rm 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_{\rm
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 which are not generally found
in the observational literature. Across the grid, we have mean values of about
3300 km/s and $\varv_{\infty}/\varv_{\rm esc,eff}\approx 4.5$ (where
$\varv_{\infty}/\varv_{\rm
esc,eff}=\sqrt{\frac{2GM_{\ast}(1-\Gamma_{e})}{R_{\ast}}}$), which corresponds
to a value of $\varv_{\infty}/\varv_{\rm esc}\approx 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 $\varv_{\infty}$ 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 $\dot{M}$ 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_{\rm 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
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.

Table 1: Parameters and prescriptions of a massive-star model in mesa. See Tab. 1 in Choi et al. (2016) for an extended list of the prescriptions. Ingredient | Adopted physics | Details
---|---|---
Mass loss | This paper |  $T_{\rm eff}>11$kK, $X_{\rm surf}>0.4$
| de Jager et al. |  $T_{\rm eff}<10$kK
| Nugis & Lamers |  $T_{\rm eff}>10$kK, $X_{\rm surf}<0.4$
Overshoot | $f_{\rm ov,core}=0.016$ | exponential overshoot
| $f_{\rm ov,envelope}=0.0174$ |
Convection | Ledoux + MLT | Includes superadiabatic
|  | convection (MLT++)
Rotation |  $\varv_{\rm init}=350$ km/s | $\Omega/\Omega_{\rm crit}=0.4$

Concerning the mass loss, we apply three different prescriptions depending on
temperature and evolution phase. For hot ($T_{\rm eff}>11$kK), hydrogen-rich
phases we apply our mass-loss formula (4). This implementation is different to
previous stellar evolution calculation, 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 HR-diagram ($T_{\rm eff}<10$kK), we switch
to the empirical rates by de Jager et al. (1988). Even though our grid only
extends to $T_{\rm 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_{\rm eff}>10$kK), stripped (with a surface hydrogen mass fraction
$X_{\rm 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_{\rm
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
$\Omega=0.4\Omega_{\rm crit}$, with $\Omega_{\rm 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 equivalent 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_{\ast}=60$ M⊙
and $Z_{\ast}=Z_{\odot}$, and here below we now discuss two important aspects
that are seen in these stellar tracks resulting from our new mass-loss rates.

Figure 6: HR-diagram for a star with an initial mass of $M_{\ast}=60$ M⊙and an
initial surface rotation of 350 km/s. The two tracks are calculated assuming
different mass-loss prescriptions.

####  5.3.1 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 $\approx 3$ on the main
sequence and $\approx 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_{\rm
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 $>100\text{M}_{\odot}$ in
order to form an evolved, hot WR helium star with our new mass-loss
rates555These models applied an extrapolation of the mass-loss rates beyond
the grid at high mass and luminosity, where as discussed above our
prescription might not fully capture the nature of these stars’ winds..
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, 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 regime666even though some
first 3D radiation-hydrodynamical simulations of LBV envelopes by Jiang et al.
(2018) only found moderate mass-loss rates, not significantly higher than
those predicted by our steady-state recipe here..

####  5.3.2 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 He-core 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-45\text{M}_{\odot}\ $ (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 mass-loss
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. 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_{\ast}\leq 0.1Z_{\odot}$ 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 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. (4). As in the previous O-star calculations by Björklund
et al. (2021), we find strong trends of the mass-loss rate with $L_{\ast}$ and
$Z_{\ast}$, which are supplemented here with additional dependencies on
$M_{\rm eff}$ and $T_{\rm 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
OB-star 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. (2000,
2001), we again find a downwards 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 $\beta$-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$\alpha$ (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.

###### Acknowledgements.

RB and JOS acknowledge support from the Odysseus program of the Belgian
Research Foundation Flanders (FWO) under grant G0H9218N. We thank the referee
for their constructive and useful comments on the manuscript.

## References

  * Abbott et al. (2016a) Abbott, B. P., Abbott, R., Abbott, T. D., et al. 2016a, ApJ, 818, L22
  * Abbott et al. (2016b) Abbott, B. P., Abbott, R., Abbott, T. D., et al. 2016b, ApJ, 818, L22
  * Asplund et al. (2009) Asplund, M., Grevesse, N., Sauval, A. J., & Scott, P. 2009, ARA&A, 47, 481
  * Belczynski et al. (2020) Belczynski, K., Hirschi, R., Kaiser, E. A., et al. 2020, ApJ, 890, 113
  * Benaglia et al. (2007) Benaglia, P., Vink, J. S., Martí, J., et al. 2007, A&A, 467, 1265
  * Björklund et al. (2021) Björklund, R., Sundqvist, J. O., Puls, J., & Najarro, F. 2021, A&A, 648, A36
  * Bloecker (1995) Bloecker, T. 1995, A&A, 297, 727
  * Böhm-Vitense (1958) Böhm-Vitense, E. 1958, ZAp, 46, 108
  * Bond et al. (1984) Bond, J. R., Arnett, W. D., & Carr, B. J. 1984, ApJ, 280, 825
  * Bouret et al. (2012) Bouret, J.-C., Hillier, D. J., Lanz, T., & Fullerton, A. W. 2012, A&A, 544, A67
  * Brott et al. (2011) Brott, I., de Mink, S. E., Cantiello, M., et al. 2011, A&A, 530, A115
  * Castor et al. (1975) Castor, J. I., Abbott, D. C., & Klein, R. I. 1975, ApJ, 195, 157
  * Choi et al. (2016) Choi, J., Dotter, A., Conroy, C., et al. 2016, ApJ, 823, 102
  * Cohen et al. (2014) Cohen, D. H., Wollman, E. E., Leutenegger, M. A., et al. 2014, MNRAS, 439, 908
  * Crowther et al. (2006) Crowther, P. A., Lennon, D. J., & Walborn, N. R. 2006, A&A, 446, 279
  * de Jager et al. (1988) de Jager, C., Nieuwenhuijzen, H., & van der Hucht, K. A. 1988, A&AS, 72, 259
  * Dotter (2016) Dotter, A. 2016, ApJS, 222, 8
  * Driessen et al. (2019) Driessen, F. A., Sundqvist, J. O., & Kee, N. D. 2019, A&A, 631, A172
  * Eggenberger et al. (2008) Eggenberger, P., Meynet, G., Maeder, A., et al. 2008, Ap&SS, 316, 43
  * Fierlinger et al. (2016) Fierlinger, K. M., Burkert, A., Ntormousi, E., et al. 2016, MNRAS, 456, 710
  * Fryer & Kalogera (2001) Fryer, C. L. & Kalogera, V. 2001, ApJ, 554, 548
  * Glatzel et al. (1985) Glatzel, W., Fricke, K. J., & El Eid, M. F. 1985, A&A, 149, 413
  * Gormaz-Matamala et al. (2021) Gormaz-Matamala, A. C., Curé, M., Hillier, D. J., et al. 2021, ApJ, 920, 64
  * Hawcroft et al. (2021) Hawcroft, C., Sana, H., Mahy, L., et al. 2021, å, submitted
  * Heger et al. (2000) Heger, A., Langer, N., & Woosley, S. E. 2000, ApJ, 528, 368
  * Heger & Woosley (2002) Heger, A. & Woosley, S. E. 2002, ApJ, 567, 532
  * Hillier & Miller (1998) Hillier, D. J. & Miller, D. L. 1998, ApJ, 496, 407
  * Howarth et al. (1997) Howarth, I. D., Siebert, K. W., Hussain, G. A. J., & Prinja, R. K. 1997, MNRAS, 284, 265
  * Jiang et al. (2018) Jiang, Y.-F., Cantiello, M., Bildsten, L., et al. 2018, arXiv e-prints, arXiv:1809.10187
  * Kee et al. (2021) Kee, N. D., Sundqvist, J. O., Decin, L., de Koter, A., & Sana, H. 2021, A&A, 646, A180
  * Keszthelyi et al. (2017) Keszthelyi, Z., Puls, J., & Wade, G. A. 2017, A&A, 598, A4
  * Krtička & Kubát (2010) Krtička, J. & Kubát, J. 2010, A&A, 519, A50
  * Krtička & Kubát (2018) Krtička, J. & Kubát, J. 2018, A&A, 612, A20
  * Krtička et al. (2021) Krtička, J., Kubát, J., & Krtičková, I. 2021, A&A, 647, A28
  * Lagae et al. (2021) Lagae, C., Driessen, F. A., Hennicker, L., Kee, N. D., & Sundqvist, J. O. 2021, A&A, 648, A94
  * Langer (1998) Langer, N. 1998, A&A, 329, 551
  * Lucy (1971) Lucy, L. B. 1971, ApJ, 163, 95
  * Maeder & Meynet (2000) Maeder, A. & Meynet, G. 2000, ARA&A, 38, 143
  * Markova & Puls (2008) Markova, N. & Puls, J. 2008, A&A, 478, 823
  * Meynet et al. (1994) Meynet, G., Maeder, A., Schaller, G., Schaerer, D., & Charbonnel, C. 1994, A&AS, 103, 97
  * Muijres et al. (2011) Muijres, L. E., de Koter, A., Vink, J. S., et al. 2011, A&A, 526, A32
  * Müller & Vink (2008) Müller, P. E. & Vink, J. S. 2008, A&A, 492, 493
  * Najarro et al. (2011) Najarro, F., Hanson, M. M., & Puls, J. 2011, A&A, 535, A32
  * Nugis & Lamers (2000) Nugis, T. & Lamers, H. J. G. L. M. 2000, A&A, 360, 227
  * Ober et al. (1983) Ober, W. W., El Eid, M. F., & Fricke, K. J. 1983, A&A, 119, 61
  * Owocki (2015) Owocki, S. P. 2015, Instabilities in the Envelopes and Winds of Very Massive Stars, Vol. 412, 113
  * Paxton et al. (2011) Paxton, B., Bildsten, L., Dotter, A., et al. 2011, ApJS, 192, 3
  * Paxton et al. (2013) Paxton, B., Cantiello, M., Arras, P., et al. 2013, ApJS, 208, 4
  * Paxton et al. (2015) Paxton, B., Marchant, P., Schwab, J., et al. 2015, ApJS, 220, 15
  * Paxton et al. (2018) Paxton, B., Schwab, J., Bauer, E. B., et al. 2018, ApJS, 234, 34
  * Paxton et al. (2019) Paxton, B., Smolec, R., Schwab, J., et al. 2019, ApJS, 243, 10
  * Petrov et al. (2016) Petrov, B., Vink, J. S., & Gräfener, G. 2016, MNRAS, 458, 1999
  * Poniatowski et al. (2021) Poniatowski, L. G., Sundqvist, J. O., Kee, N. D., et al. 2021, A&A, 647, A151
  * Puls et al. (2020) Puls, J., Najarro, F., Sundqvist, J. O., & Sen, K. 2020, A&A, 642, A172
  * Puls et al. (2008) Puls, J., Vink, J. S., & Najarro, F. 2008, A&A Rev., 16, 209
  * Reimers (1975) Reimers, D. 1975, Memoires of the Societe Royale des Sciences de Liege, 8, 369
  * Rubio-Díez et al. (2021) Rubio-Díez, M. M., Sundqvist, J. O., Najarro, F., et al. 2021, accepted by A&A
  * Sander et al. (2017) Sander, A. A. C., Hamann, W.-R., Todt, H., Hainich, R., & Shenar, T. 2017, A&A, 603, A86
  * Sander et al. (2020) Sander, A. A. C., Vink, J. S., & Hamann, W. R. 2020, MNRAS, 491, 4406
  * Schultz et al. (2020) Schultz, W. C., Bildsten, L., & Jiang, Y.-F. 2020, ApJ, 902, 67
  * Shenar et al. (2020) Shenar, T., Gilkis, A., Vink, J. S., Sana, H., & Sander, A. A. C. 2020, A&A, 634, A79
  * Shenar et al. (2016) Shenar, T., Hainich, R., Todt, H., et al. 2016, A&A, 591, A22
  * Shenar et al. (2015) Shenar, T., Oskinova, L., Hamann, W. R., et al. 2015, ApJ, 809, 135
  * Simón-Díaz & Herrero (2014) Simón-Díaz, S. & Herrero, A. 2014, A&A, 562, A135
  * Smith (2014) Smith, N. 2014, ARA&A, 52, 487
  * Sundqvist et al. (2019) Sundqvist, J. O., Björklund, R., Puls, J., & Najarro, F. 2019, A&A, 632, A126
  * Sundqvist et al. (2011) Sundqvist, J. O., Puls, J., Feldmeier, A., & Owocki, S. P. 2011, A&A, 528, 64
  * Sundqvist et al. (2014) Sundqvist, J. O., Puls, J., & Owocki, S. P. 2014, A&A, 568, 59
  * Šurlan et al. (2013) Šurlan, B., Hamann, W.-R., Aret, A., et al. 2013, A&A, 559, A130
  * Vink et al. (1999) Vink, J. S., de Koter, A., & Lamers, H. J. G. L. M. 1999, A&A, 350, 181
  * Vink et al. (2000) Vink, J. S., de Koter, A., & Lamers, H. J. G. L. M. 2000, A&A, 362, 295
  * Vink et al. (2001) Vink, J. S., de Koter, A., & Lamers, H. J. G. L. M. 2001, A&A, 369, 574
  * Vink et al. (2021) Vink, J. S., Higgins, E. R., Sander, A. A. C., & Sabhahit, G. N. 2021, MNRAS, 504, 146
  * Vink & Sander (2021) Vink, J. S. & Sander, A. A. C. 2021, MNRAS, 504, 2051
  * Woosley (2017) Woosley, S. E. 2017, ApJ, 836, 244
  * Woosley et al. (2007) Woosley, S. E., Blinnikov, S., & Heger, A. 2007, Nature, 450, 390

##  Appendix A Additional figures

Fig. 7 shows a similar comparison as in Fig. 3 of a self-consistent model and
a model using a $\beta$-velocity law of $\beta=1$, but now at higher mass and
luminosity. The mass-loss rate of the $\beta$-model is about 70 times higher
than the dynamically-consistent value. This model then again shows a mismatch
between $\Gamma$ and $\Lambda$ with a value of $\Gamma=0.45$ at the sonic
points, which is significantly below the required value $\approx 1$.

Figure 7: Same figure as Fig. 3 but now with
$\log(L_{\ast}/\text{L}_{\odot})=5.75$, $M_{\ast}=50\text{M}_{\odot}$, and
$T_{\rm eff}=15000$K. The self-consistent model has $\dot{M}=7.67\cdot
10^{-8}$ M⊙/yr and $\varv_{\infty}=3111$ km/s, while the $\beta$-model has
$\dot{M}=5.38\cdot 10^{-6}$ M⊙/yr and $\varv_{\infty}=451$ km/s.

##  Appendix B Comparison with $\Gamma$ from CMFGEN

Table 2: Stellar and wind parameters of our model grid used to check the radiative acceleration in the B super-/hypergiant domain. All models have a luminosity of $\log L_{\ast}/\text{L}_{\odot}=5.76$, and have been calculated with a velocity field exponent $\beta=1$, a micro-turbulent velocity $\varv_{\text{turb}}=10\leavevmode\nobreak\ {\rm km}\,{\rm s}^{-1}$, an unclumped wind, no X-ray emission from wind-embedded shocks, and solar abundances from Asplund et al. (2009), in particular $Y_{\rm He}$=$N_{\rm He}/N_{\rm H}$ = 0.1. Model | Mass | $T_{\rm eff}$ | $R_{\ast}$ | $\log{g}$ | $\Gamma_{\rm e}$ | $\dot{M}$ | $\varv_{\infty}$
---|---|---|---|---|---|---|---
|  M⊙ | (K) | (R⊙) | (cgs) |  | (M⊙yr-1) | (${\rm km}\,{\rm s}^{-1}$)
B225_1 | 31.5 | 22500 | 49.3 | 2.55 | 0.43 | 8.20$\cdot 10^{-6}$ | 678
B225_2 | 31.5 | 22500 | 49.3 | 2.55 | 0.43 | 0.73$\cdot 10^{-6}$ | 678
B200_1 | 30.1 | 20000 | 62.4 | 2.33 | 0.45 | 11.20$\cdot 10^{-6}$ | 556
B200_2 | 50.5 | 20000 | 62.4 | 2.55 | 0.27 | 3.70$\cdot 10^{-6}$ | 560
B175_1 | 33.0 | 17500 | 71.3 | 2.25 | 0.31 | 14.30$\cdot 10^{-6}$ | 397
B175_2 | 33.0 | 17500 | 71.3 | 2.25 | 0.31 | 3.00$\cdot 10^{-6}$ | 430
B150_1 | 37.8 | 15000 | 83.2 | 2.17 | 0.20 | 14.30$\cdot 10^{-6}$ | 397
B150_2 | 45.0 | 15000 | 83.2 | 2.25 | 0.17 | 3.50$\cdot 10^{-6}$ | 430

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

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_{\rm eff}$$\approx$ 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_{\rm
eff}$ $\approx$ 20 kK and concentrating on B super-/hypergiants, assuming a
prototypical $\beta=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 $\beta$ law and the quasi-
hydrostatic stratification was set (in most cases) to 0.5 $\varv_{\rm 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 2.

The outcome of our modelling is displayed in Fig. 8, 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 $\Gamma$ as defined in Eq. 1). The red dashed-dotted lines indicate
the conventional Gamma-factor for electron scattering, $\Gamma_{\rm 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 $\Gamma$ 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$\varv_{\infty}$) 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 (where $\dot{M}$ is determined) until few hundreds of
${\rm km}\,{\rm 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 calculations presented in the main part of this paper
are not affected by too low $g_{\text{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 $\Gamma_{\rm tot}=1$ is reached only at
substantial velocities (100 ${\rm km}\,{\rm s}^{-1}$or more), whereas, in
self-consistent wind-models, this should happen very close to the sonic point
(here: 15 … 18 ${\rm km}\,{\rm s}^{-1}$, in dependence of $T_{\rm 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_{\text{rad}}\propto({\rm d}\varv/{\rm d}r)/\rho$ (e.g., Castor
et al. 1975), only a steeper velocity law and a lower mass-loss rate will lead
to a higher $g_{\text{rad}}$ in the transonic region, which is also the basic
outcome of our self-consistent calculations.