Issue 
A&A
Volume 648, April 2021



Article Number  A36  
Number of page(s)  16  
Section  Stellar atmospheres  
DOI  https://doi.org/10.1051/00046361/202038384  
Published online  08 April 2021 
New predictions for radiationdriven, steadystate massloss and windmomentum from hot, massive stars
II. A grid of Otype stars in the Galaxy and the Magellanic Clouds
^{1}
KU Leuven, Instituut voor Sterrenkunde,
Celestijnenlaan 200D, 3001 Leuven, Belgium
email: robin.bjorklund@kuleuven.be
^{2}
LMU München, Universitätssternwarte,
Scheinerstr. 1,
81679 München, Germany
^{3}
Centro de Astrobiologia, Instituto Nacional de Tecnica Aerospacial,
28850
Torrejon de Ardoz,
Madrid, Spain
Received:
8
May
2020
Accepted:
12
August
2020
Context. Reliable predictions of massloss rates are important for massivestar evolution computations.
Aims. We aim to provide predictions for massloss rates and windmomentum rates of Otype stars, while carefully studying the behaviour of these winds as functions of stellar parameters, such as luminosity and metallicity.
Methods. We used newly developed steadystate models of radiationdriven winds to compute the global properties of a grid of Ostars. The selfconsistent models were calculated by means of an iterative solution to the equation of motion using full nonlocal thermodynamic equilibrium radiative transfer in the comoving frame to compute the radiative acceleration. In order to study winds in different galactic environments, the grid covers mainsequence stars, giants, and supergiants in the Galaxy and both Magellanic Clouds.
Results. We find a strong dependence of massloss on both luminosity and metallicity. Mean values across the grid are Ṁ~L_{*}^{2.2} and Ṁ~L_{*}^{0.95}; however, we also find a somewhat stronger dependence on metallicity for lower luminosities. Similarly, the mass lossluminosity relation is somewhat steeper for the Small Magellanic Cloud (SMC) than for the Galaxy. In addition, the computed rates are systematically lower (by a factor 2 and more) than those commonly used in stellarevolution calculations. Overall, our results are in good agreement with observations in the Galaxy that properly account for windclumping, with empirical Ṁ versus Z_{*} scaling relations and with observations of Odwarfs in the SMC.
Conclusions. Our results provide simple fit relations for massloss rates and wind momenta of massive Ostars stars as functions of luminosity and metallicity, which are valid in the range T_{eff} = 28 000–45 000 K. Due to the systematically lower values for Ṁ, our new models suggest that new rates might be needed in evolution simulations of massive stars.
Key words: stars: atmospheres / stars: earlytype / stars: massive / stars: massloss / stars: winds, outflows / Magellanic Clouds
© ESO 2021
1 Introduction
Hot, massive stars with masses ≿9 M_{⊙} of spectral type O and B lose a significant amount of mass due to their radiationdriven stellar winds (Castor et al. 1975). This mass loss has a dominant influence on the life cycles of massive stars as well as in determining the properties of the remnants left behind when these stars die (e.g. Smith 2014). The rates at which these stars lose mass, which is of the order of Ṁ ~ 10^{−5…−9} M_{⊙} yr^{−1}, comprise a key uncertainty in current models of stellar evolution (even on the main sequence where the stars are typically “well behaved”, e.g. Keszthelyi et al. 2017), simply because the mass of a star is the most important parameter when determining its evolution. In addition to the loss of mass, angular momentum is also lost through stellar winds affecting the surface rotation speeds of these stars. Moreover, uncertainties related to mass loss have consequences on galactic scales beyond stellar physics, as massive stars provide strong mechanical and radiative feedback to their host environment (Bresolin et al. 2008).
It is therefore important to have reliable quantitative predictions of massloss rates and windmomenta of massive stars. In the first paper of this series (Sundqvist et al. 2019, from here on Paper I), we developed a new method to provide massloss predictions based on steadystate wind models using radiative transfer in the comoving frame (CMF) to compute the radiative acceleration g_{rad}. Building on that, this paper now presents the results from a full grid of models computed for Ostars in the Galaxy, as well as the Small and the Large Magellanic Cloud, analysing the general dependence on important stellar quantities, such as luminosity and metallicity. Paper I shows that simulations using CMF radiative transfer suggest reduced massloss rates as compared to the predictions normally included in models of massive star evolution (Vink et al. 2000, 2001). Although this paper focuses on the presentation of our new rates for Ostars in different galactic environments, a key aim for future publications within this series will be to directly implement the results of the new wind models into calculations of massivestar evolution.
Since most of the important spectral lines driving hotstar winds are metallic, a strong dependence on metallicity Z_{*} is expected for the massloss rate (Kudritzki et al. 1987; Vink et al. 2001; Mokiem et al. 2007). To investigate this, here, we compute models tailored to our local Galactic environment, assuming a metal content similar to that in the Sun, as well as for the Magellanic Clouds. The metallicities of these external galaxies are about half the Galactic one for the Large Magellanic Cloud (LMC) and a fifth for the Small Magellanic Cloud (SMC). Using these three regimes, we aim to study the massloss rate as a function of Z_{*}. The Magellanic Clouds are interesting labs for stellar astrophysics because the distances to the stars are relatively well constrained, providing values of their luminosities and radii. Another reason we are focusing on the Magellanic Clouds is that quantitative spectroscopy of individual stars there has been performed and compiled into an observed set of scalingrelations for windmomenta (Mokiem et al. 2007). While the quantitative spectroscopy of individual, hot, massive stars is also possible nowadays in galaxies further away (e.g. Garcia et al. 2019), such studies are only in their infancy.
In order to derive global dependencies and relations for the massloss rates and windmomenta, we perform a study using a grid of Ostar models. Thanks to the fast performance of the method, as explained in detail in Paper I, this is finally possible for the hydrodynamically consistent steadystate models with a nonparamterised CMF lineforce computed without any assumptions about underlying linedistribution functions. In Sect. 2, we briefly review our method for computing massloss rates, highlighting one representative model from the grid. In Sect. 3, the results of the full grid of models are shown, first for the Galaxy and then including the Magellanic Clouds, in terms of computed windmomenta and massloss rates. From these results we derive simple fit relations for the dependence on luminosity and metallicity. Section 4 provides a discussion of the results, highlighting the general trends and comparing to other existing models and to observations. Additionally we address the implication for stellar evolution and current issues such as the socalled weak wind problem. Section 5 contains the conclusions and future prospects.
2 Methods
A crucial part about the radiationdriven steadystate wind models used in our research is that they are hydrodynamically consistent. This means that the equation of motion (e.o.m.) in the spherically symmetric, steadystate case is solved as described in Paper I. This e.o.m. reads (1)
Here, v(r) is the velocity, a(r) the isothermal sound speed, g_{rad}(r) the radiative acceleration, g(r) = GM_{*}∕r^{2} the gravity, with gravitation constant G, M_{*} the stellar mass used in the model, and r the radius coordinate. The temperature structure T(r) enters the equation through the isothermal sound speed (2)
with k_{B} Boltzmann’s constant, μ(r) the mean molecular weight, and m_{H} the mass ofa hydrogen atom. Equation (1) has a sonic point where v(r) = a(r). Because in the above formulation g_{rad} is only an explicit function of radius, and not of the velocity gradient, the corresponding critical point in the e.o.m. is this sonic point (also see the discussion in Paper I). The radiative acceleration also depends on velocity and mass loss, of course, but in our method these dependencies are implicit; they are accounted for through the iterative updates of the velocity and density structure and do not affect the critical point condition in an explicit way.
For given stellar parameters luminosity L_{*}, mass M_{*}, radius R_{*} (to be defined below), and metallicity Z_{*}, the e.o.m. (1) is solved to obtain v(r) for the subsonic photosphere and supersonic radiationdriven wind. For a steadystate massloss rate Ṁ the mass conservation equation Ṁ = 4πr^{2}ρv gives the density structure ρ(r). The wind modelsfurther rely on the NLTE (nonlocal thermodynamic equilibrium) radiative transfer in FASTWIND (see Paper I and Puls 2017) for the computation of g_{rad}, by means ofa comoving frame (CMF) solution and without using any parametrised distribution functions. The atomic data are taken from the WMBASIC data base (Pauldrach et al. 2001). This compilation consists of more than a million spectral lines, including all major metallic elements up to Zn and all ionisation stages relevant for Ostars. The data base is the same as the one utilised in Paper I, as well as in previous versions of the FASTWIND code (see, e.g. Puls et al. 2005). Also the hydrogen and helium model atoms are identical to those used in our previous work. We note that since WMBASIC is calibrated for diagnostic usage in the UV regime, its principal database should be ideally suited for the radiative force calculations in focus here. In the NLTE and g_{rad} calculation, we account for pure Dopplerbroadening alone^{1} as depth and mass dependent profiles, including also a fixed microturbulent velocity v_{turb} (see further Paper I, and Sects. 2.3 and 4.5 of this paper). More specific details of the NLTE and radiative transfer in the new CMF FASTWIND v11 are laid out in detail in Puls et al. (2020, see also Puls 2017 and Paper I).
The steadystate Ṁ and v(r) are converged in the model computation starting from a first guess. For all simulations presented in this paper, the startvalue for Ṁ was taken as the massloss rate predicted by the Vink et al. (2001) recipe using the stellar input parameters of the model. The initial velocity structure is obtained by assuming that a quasihydrostatic atmosphere connects at v_{tr} ≈ 0.1a(T = T_{eff}) to a socalled βvelocity law (3)
with v_{∞} the terminal wind speed, R_{*} the stellar radius, β a positive exponent, and b a constant derived from the transition velocity v_{tr}. We further define the stellar radius (4)
where is the spherically modified fluxweighted optical depth (5)
for the fluxweighted opacity κ_{F} (cm^{2} g^{−1}). The fluxweighted opacity is related to the radiative acceleration as g_{rad} = κ_{F}L_{*}∕(4πcr^{2}). After each update of the hydrodynamical structure (see below), an NLTE/radiative transfer loop is carried out, to converge the occupation numbers and g_{rad}. The velocity gradient at the wind critical point is then computed by applying l’Hôpital’s rule, after which the momentum Eq. (1) is solved with a Runge–Kutta method to obtain the velocity structure v(r) above and below the critical point. This integration is performed by shooting both outwards from the sonic point to a radius of about 100R_{*} to 140R_{*} and inwards towards the star, stopping at r = r_{min} when a column mass g cm^{−2} is reached.
In addition to the velocity structure, the temperature structure is also updated every hydrodynamic iteration. We use a simplified method similar to Lucy (1971) to speed up the convergence; the temperature throughout the radial grid is calculated as (6)
where W(r) is the dilution factor given by (7)
The effective temperature T_{eff} in Eq. (6) is defined as , with σ the Stefan–Boltzman constant and R_{*} as defined in Eq. (4). Additionally there is a floor temperature of T ≈ 0.4T_{eff} such as in previous versions of FASTWIND. This temperature structure is held fixed during the following NLTE iteration, meaning that, formally, perfect radiative equilibrium is not achieved; however, the effects from this on the wind dynamics are typically negligible (see Paper I).
As described in Paper I, the singularity and regularity conditions applied in CAKtheory cannot be used to update the massloss rate because in the approach considered here g_{rad} does not explicitly depend on density, velocity, or the velocity gradient. Instead, here at iteration i the massloss rate used in iteration i + 1 is updated to counter the current mismatch in the force balance at the critical sonic point (see also Sander et al. 2017). For a hydrodynamically consistent solution, the quantity (8)
should be equal to Γ = g_{rad}∕g at the sonic point. In order to fulfill the e.o.m. (1) the current mismatch is countered by updating the massloss rate according to , following the basic theory of linedriven winds where g_{rad} ∝ 1∕Ṁ^{b} (Castor et al. 1975). Our models take a value of b = 1 in the iteration loop, providing a stable way to converge the steadystate massloss rate. From this Ṁ and the computed velocity field, a new density is obtained from the mass conservation equation.
2.1 Convergence
The new wind structure (v(r), ρ(r), T(r), Ṁ) is used in the next NLTE iteration loop to converge the radiative acceleration once more in the radiativetransfer scheme. This method is then iterated until the error in the momentum equation is small enough to consider the model as converged and thus hydrodynamically consistent. By rewriting the e.o.m. (1) the quantity describing the current error is (9)
for each radial position r. For a hydrodynamically consistent model f_{err} is zero everywhere and Γ should thus be equal to Λ. As such, in the models we need the maximum error in the radial grid (11)
to be close enough to zero; in this paper we require a threshold 0.01, meaning the converged model is dynamically consistent to within a percent. Additional convergence criteria apply to the massloss rate and the velocity structure. These quantities are not allowed to vary by more than 2% for the former and 3% for the latter between the last two hydrodynamic iteration steps reaching convergence.
2.2 Generic model outcome
As a first illustration, some generic outcomes of one characteristic simulation are now highlighted; the model parameters are listed in Table 1. The top panel of Fig. 1 shows the evolution of Γ for several iterations in the scheme starting from an initial βlaw structure. The characteristic pattern of a steep wind acceleration starting around the sonic point, where Γ ≈1, is clearly visible throughout the iteration loop. The model starts off being far from consistency (yellow), but as the error in the hydrodynamical structure becomes smaller the solution eventually relaxes to a final converged velocity structure (dark blue). The innermost points deepest in the photosphere remain quite constant, since the deep photospheric layers relax relatively quickly. In the bottom panel of Fig. 1 a colour plot of the iterative evolution of the model error f_{err} throughout the wind can be seen. The point of maximum error at each iteration is marked with a plus and the dashdotted line shows where the velocity equals the sound speed. The dashed lines further show the boundaries within which is computed, where we note that the part at very low velocity is excluded because here the opacity is parametrised (see below). In addition a few of the outermost points are formally excluded in the calculation of (due to resolution considerations). Since the calculation of excludes the innermost region and a few outermost points, these points do not contribute to the condition of convergence based on the error in f_{err}. Nonetheless the models do provide reliable terminal wind speeds, as they additionally require the complete velocity structure (including v_{∞}) to be converged to better than 3% between the final two iteration steps (see above). The figure illustrates explicitly how both the overall and maximum errors generally decrease throughout the iteration cycle of the simulation. We note that, after some initial relaxation, for this particular model the position of maximum error always lies in the supersonic region, often quite close to the critical sonic point.
At the end of the sequence, the model is dynamically consistent and Γ matches the other terms in the equationof motion Λ. This is illustrated in Fig. 2. This figure compares Γ and Λ at each radial point of the converged model, showing a clear match between the quantities. Only below a velocity v ≲ 0.1 km s^{−1} there is some discrepancy; this arises because in these quasistatic layers the fluxweighted opacity is approximated by a Kramerlike parametrisation (see Paper I), which is useful to stabilise the base in the deep subsonic atmosphere. It is important to point out that this parametrisation is applied only at low velocities and high optical depths and so does not affect the structure of the wind or the derived global parameters.
The behaviour of the massloss rate Ṁ is important to understand for the purposes of this paper. The top panel of Fig. 3 shows of the model for all iterations versus the massloss rate computed for that iteration. The general trend is that decreases quite consistently during the iteration cycle. The massloss rate can be seen to converge to one value as the structure gets closer to dynamical consistence. Indeed, in the last couple of iterations the value of Ṁ only changes minimally from its former value. In the bottom panel of Fig. 3, the same plot is shown for the iterative evolution of the terminal wind speed. Also this quantity displays a stable convergence behaviour towards one final value. The quantities Ṁ and v_{∞} are from Fig. 3 seen to have an anticorrelation, as for this model the total windmomentum rate Ṁ v_{∞} does not vary much after the first few iterations. Finally, Fig. 4 shows the converged velocity structure versus the modified radius coordinate , where r_{min} is the innermost radial point of the grid. This figure illustrates the very steep acceleration around the sonic point that is characteristic for our Ostar models, followed by a supersonic βlawlike behaviour typical of a radiationdriven stellar wind. In the figure we add a fit to the velocity structure using a double βlaw defined by Eq. (22) and elaborate further on this comparison in Sect. 4.5.
Fig. 1 Top panel: value of Γ versus scaled radiuscoordinate for 7 (nonconsecutive) hydrodynamic iterations over the complete run. The starting structure (yellow) relaxes to the final converged structure (dark blue). Bottom panel: colour map of log (f_{err}) for all hydrodynamic iterations; on the abscissa is hydrodynamic iteration number and on the ordinate scaled wind velocity. The pluses indicate the location of for each iteration, the dashed lines the limits between which is computedand the dashdotted line the location of the sonic point. 
Fig. 2 Top panel: final converged structure of the characteristic model showing Γ in black squares and Λ (see text for definition) as a green line. The black dashed lines show the location of the sonic point approximately at Γ = 1 (but not exactly because of the additional pressure terms). Bottom panel: same as in the top panel, however versus velocity (which resolves the inner wind more). 
2.3 Grid setup
In order to study the massloss rates of O stars in a quantitative way, a modelgrid was constructed by varying the fundamental input stellarparameters. For the Galactic stars, we used the calibrated stellar parameters obtained from a theoretical T_{eff} scale of a set of O stars from Martins et al. (2005; adopting the values from their Table 1). The Martins et al. parameters are retrieved by means of a grid of nonLTE, lineblanketed synthetic spectra models using the code CMFGEN (Hillier & Miller 1998). Here, a selfconsistent wind model as described above was calculated for the stellar parameters of each star in the grid, resulting in predictions of the wind structure, terminal velocity, and massloss rate. The wind models were computed with the microturbulent velocity kept at a standard value for Ostars of v_{turb} = 10 km s^{−1} (see also Paper I) and the helium number abundance was kept fixed at Y_{He} = n_{He}∕n_{H} = 0.1. Moreover, the simulations were performed without any inclusion of clumping or highenergy Xrays. As discussed in detail in Paper I, while such wind clumping is both theoretically expected (Owocki et al. 1988; Sundqvist et al. 2018; Driessen et al. 2019)and observationally established (see Sect. 4.2), it is still uncertain what effect this might have on theoretically derived global massloss rates (indeed, in the simple tests performed in Paper I the effect was only marginal). In addition, any inclusion of wind clumping into steadystate models will inevitably be of adhoc nature (see also discussion in Paper I). A study by Muijres et al. (2011), for example, shows that introducing clumps can sometimes change their predicted massloss rate by as much as an order of magnitude. The models used in their study, however, use a global energy constraint to derive the massloss rate for an assumed fixed β velocity law. As such, they are not locally consistent and thus might not necessarily fulfill the force requirements around the sonic point. By contrast,the models presented here are (by design) both locally and globally consistent, with a massloss rate that is primarily sensitive to the conditions around the critical sonic point. It might only be influenced if these corresponding regions where Ṁ is determined, are strongly clumped. This is a key difference between the massloss rates derived here and those in Muijres et al. (2011) and a prime reason that, contrary to their findings, Ṁ in our models does not seem to change significantly when introducing clumping in the supersonic parts. However, as also discussed in Paper I, the terminal velocities are typically affected by adding such clumping. Namely, since g_{rad} is altered in the supersonic regions due to the presence of the clumps, this can lead to modified values of v_{∞}. So even though Ṁ is barely influenced when including typical windclumping, it remains an uncertainty of the current models and future work should aim for a more systematic study of possible feedback effects from clumping upon also the steadystate e.o.m.
In any case, the models presented here contain 12 spectroscopic dwarfs, 12 giants and 12 supergiants for each value of metallicity Z_{*}. For simplicity, the same stellar parameters (L_{*}, M_{*}, R_{*}) as for the Galactic grid were assumed to create models for the LMC and SMC, changing only their metallicity. This setup has the advantage of enabling a rather direct comparison for the modeldependence on metallicity, independent of the other input parameters. The used metallicities in the grid are Z_{Galaxy} = Z_{⊙}, Z_{LMC} = 0.5 Z_{⊙}, and Z_{SMC} = 0.2 Z_{⊙}, respectively. The value of the Solar metallicity was here taken to be Z_{⊙} = 0.013 (Asplund et al. 2009). In total this gives 108 models, with input stellar parameters as listed in Table A.1.
Fig. 3 Top panel: iterative behaviour of the massloss rate as decreases towards a value below 1%. The colour signifies the iteration number starting from Ṁ as predicted by the Vink et al recipe in light green. Bottom panel: iterative behaviour of the terminal velocity v_{∞} towards convergence. 
Fig. 4 Converged velocity structure for the characteristic model of Sect. 2.2, showing velocity over terminal wind speed versus the scaled radiuscoordinate in black. The green line shows a fit using a double βlaw following Eq. (22) (see text for details). 
3 Results
The results for all 108 models are added to Table A.1, containing the derived values for Ṁ and v_{∞}. The runs typically took about 50 iterative updates of the hydrodynamical structure to converge, where for most parts the corresponding calculation of g_{rad} (aside from the first ones) takes 10–15 NLTE radiative transfer steps per hydrodynamic update. Using the criteria as presented in Sect. 2.1, all models presented in this paper are formally converged. The following subsections highlight the results for the Galaxy, as well as the Small and the Large Magellanic Clouds.
3.1 The Galaxy
When studying the overall behaviour of linedriven winds, it is useful to look at a modified windmomentum rate as function of the stellar luminosity: (12)
where the left hand side is the socalled modified windmomentum rate (Kudritzki et al. 1995; Puls et al. 1996), which is proportional to the luminosity to some power x.
The key advantage when using this modified windmomentum is that basic line driven wind theory predicts the dependence on M_{*} to scale out (or at least to become of only second order impact). Namely, from (modified) CAK theory the following relations can be found (e.g. Puls et al. 2008): (13)
Here, the effective escape speed from the stellar surface is for an effective stellar mass M_{eff} = M_{*}(1 − Γ_{e}), reduced by electron scattering (14)
with an opacity κ_{e} (cm^{2} g^{−1}). Equation (13) above further introduces α_{eff} = α − δ, where α, describing the power law distribution of the line strength of contributing spectral lines in CAK theory, takes values between 0 and 1 and the parameter δ accounts for ionising effects in the wind. For a simple α_{eff} = 2∕3 we thus have^{2} (15)
as in the windmomentumluminosity relation (WLR) of Eq. (12) for x = 1∕α_{eff}. However, even if α_{eff} is not exactly 2/3, the dependence on M_{*} will still be relatively weak. The validity of the WLR relation is an overall key success of linedriven wind theory and the basic concept has been observationally confirmed by a multitude of studies (see Puls et al. 2008, for an overview).
The Galactic WLR for the radiationhydrodynamic wind models presented here is shown in the top panel of Fig. 5. The figure indeed shows a quite tight relationship between the modified windmomentum rate and luminosity; fitting the models according to Eq. (12) above, a slope x = 2.07 ± 0.32 is derived (with the error mentioned being the 1σ standard deviation). Interpreted in terms of the (modified) CAK theory above, this would imply a α_{eff}≈ 0.5 for our models in the Galaxy, in rather good agreement with the typical Ostar values α ≈0.6 and δ ≈ 0.1 (Puls et al. 2008).
Next, the massloss rate versus luminosity is plotted in the bottom panel of Fig. 5. From this figure, we infer a rather steep dependence of Ṁ on L_{*}; fitting a simple powerlaw here gives y = 2.16 ± 0.34. Within our grid, we further do not find any strong systematic trends of massloss rate (or windmomentum rate) with respect to spectral type (in the considered temperature range) and luminosity class.
Figure 6 shows the terminal wind speeds for the Galactic models. We obtain a mean value v_{∞} ∕v_{esc} = 3.7 for the complete Galactic sample, however with a relatively large scatter (1σ standard deviation of 0.8). There is a systematic trend of increasing v_{∞}∕v_{esc} ratios for lower luminosities (see also Paper I), but also here the corresponding scatter is significant. Overall, although these Galactic v_{∞} ∕v_{esc} values are quite high for Ostars, the significant scatter we find is generally consistent with observational studies. Section 4.3 further addresses this, including a discussion about if a reduction in v_{∞} might also affect the prediction of Ṁ.
Fig. 5 Top panel: modified windmomentum rate versus luminosity for all Galactic models. The solid black line is a linear fit through the points and the dashed line shows the theoretical relation by Vink et al. (2000). Bottom panel: massloss rates versusluminosity for all Galactic models with a linear fit as a solid black line. The dashed line is a fit through the massloss rates computed using the Vink et al. recipe, the dashdotted line is the relation derived by Krtička & Kubát (2017) and the dotted line is the relation computed from the results from Lucy (2010). 
3.2 All models
In the top panel of Fig. 7 the modified windmomenta for the Magellanic Cloud simulations are added to those of the Galaxy, alsoincluding powerlaw fits to the models in the LMC and SMC. The bottom panel of Fig. 7 shows the same plot for Ṁ versus L_{*}. On both figures, there is a clear systematic trend that the massloss and windmomentum rates are always lower for the LMC than for the Galaxy and lower still for the SMC. This is as expected for radiationdriven Ostar winds since the majority of driving is done by metallic spectral lines. Inspection of the slope of the WLR at the LMC reveals a similar trend as for the Galaxy and we derive x = 2.12 ± 0.34 from a fit accordingto Eq. (12). On the other hand, already from simple visual inspection it is clear that for lowluminosity dwarfs in the SMC the overall slope changes significantly; indeed, an overall fit to Eq. (12) for all SMC stars here results in a higher x = 2.56 ± 0.44.
Another feature visible for the SMC supergiant models is a bump in D_{mom} at . As discussed in Sect. 4.1, this feature most likely arises because of the different effective temperatures of these models, which affects the ionisation balance of important driving elements in the wind.
The derived slopes of the massloss rate versus luminosity, , for the LMC and SMC are y = 2.17 ± 0.34 and y = 2.37 ± 0.40, respectively.We note that while the dependence for the LMC is virtually unchanged with respect to the Galaxy, the SMC models again display a steeper dependence, driven by the very low massloss rates found for the stars with the lowest luminosities. As discussed further in Sect. 4, these findings are generally consistent with the linestatistics predictions by Puls et al. (2000) that α_{eff} becomes lower both for lower density winds and for winds of lower metallicity.
Fig. 6 Terminal wind speed over photospheric effective escape speed (corrected with 1Γ_{e}, see text) versus luminosity. Values for Galactic metallicity for all three luminosity classes are shown. 
Fig. 7 Top panel: modified windmomentum rate of all models versus luminosity. The dashed lines show linear fits through each of the three sets of models. The markers show the different luminosity classes, consistent with previous plots. Bottom panel: same as the top panel, but for the massloss rate. 
3.3 Function of metallicity
Examining trends of the modified windmomentum rate for the Galaxy, the Large, and the Small Magellanic Cloud, a dependence on metallicity is next derived. As discussed above, Fig. 7 shows that the windmomenta of all models with the same metallicity follow a quite tight correlation with L_{*}. As such, to investigate the metallicity dependence we consider the three models with identical stellar parameters, but with different Z_{*}, assuming for each triplet a simple dependence (16)
The derived values of n are plotted in Fig. 8, where the distribution gives a mean value n = 0.85 with a 1σ standard deviation of 0.29; this significant scatter is not surprising since we consider only three different metallicities in the fits. However, inspection of Fig. 8 also reveals a trend of decreasing n with increasing luminosity. Approximating it here by a simple linear fit with respect to log (L_{*}∕10^{6} L_{⊙}), we find n(L_{*}) = −0.73log(L_{*}∕10^{6} L_{⊙}) + 0.46, providing an analytic approximation for the dependence of D_{mom} on Z_{*} in function of L_{*}.
The same analysis is performed also for Ṁ, assuming a dependence (17)
The distribution of the exponent m here gives amean value 0.95. The scatter around this mean is also significant, with a 1σ standard deviation of 0.21. If we again approximate the dependence of the factor m with L_{*}, we find m(L_{*}) = −0.32log(L_{*}∕10^{6} L_{⊙}) + 0.79.
Building on the combined results above, we can now construct final relations for both the modified windmomentum rate and massloss rate of the form: (18)
To obtain the fitting coefficients (which will be different for the windmomenta and massloss relations), we simply combine the scalings found in Sect. 3.1 for the Galaxy ( and ) with those found above for the metallicity dependence ( and ). Moreover, we also performed a multilinear regression where log(Ṁ) depends on log(Z_{*}), log (L_{*}), and log (Z_{*}) ⋅ log(L_{*}); the fitting coefficients found using these two alternative methods are indeed identical to the second digit. Thus the final relations are (19)
with the wind momentum rate in units of M_{⊙} yr^{−1} km s^{−1} R_{⊙}^{0.5} and the mass − loss rate in units of M_{⊙} yr^{−1}. When the complete modelgrid sample is considered, these fitted relations give mean values that agree with the original simulations to within 10%, with standard deviations 0.33 and 0.39 for the windmomentum rate and the massloss rate, respectively. For none of the models is the ratio between the actual simulation and the fit larger than a factor 2.1 or smaller than 0.36.
Finally, we can set up the same analysis to derive the dependence of the terminal wind speed on metallicity, assuming v_{∞} ∝ Z^{p(L*)}. From this we find that the exponent p also varies approximately linearly with the log of luminosity as p(L_{*}) = − 0.41 log(L_{*}∕10^{6} L_{⊙}) − 0.32. This is consistent with the results above and indeed can be alternatively obtained by simply combining the previously derived Ṁ and D_{mom} relations (because p = n − m). The linear behaviour of p here causes low luminosity stars to have a positive exponent with Z, which flattens out and gets negative for higher luminosities. So as a general trend, stars with log (L_{*}∕10^{6} L_{⊙}) above roughly −0.78 tend to have slightly decreasing v_{∞} with increasing metallicity, while for stars below −0.78 it is the other way around.
Fig. 8 Value of the exponent n showing the metallicity dependence of D_{mom} together with a linear fit, showing the linear dependence of this exponent with log (L_{*} ∕10^{6} L_{⊙}). See text for details. 
4 Discussion
4.1 General trends
The WLR results presented in the previous section show that our results would be overall consistent with standard CAK linedriven theory for α_{eff} ≈ 0.5, at least for the Galactic and LMC cases. As already mentioned, this is in reasonably good agreement with various CAK and linestatistics results (see overview in Puls et al. 2008). For the Ostars in the SMC, on the other hand, we find a steeper relation with L_{*}, implying an overall α_{eff} ≈ 0.42 if interpreted by means of such basic CAK theory. However, Fig. 7 shows that the WLR here exhibits significant curvature with a steeper dependence for lower luminosities, making interpretation in terms of a single slope somewhat problematic. The modelgrid indicates that α_{eff} decreases both with decreasing metallicity and with decreasing luminosity, suggesting that α_{eff} generally becomes lower for lower wind densities. This is consistent with the linestatistics calculations by Puls et al. (2000) and may (at least qualitatively) be understood via the physical interpretation of α as the ratio of the lineforce due to optically thick lines to the total lineforce (lower wind densities should generally mean a lower proportion of optically thick contributing lines).
We further find a quite steep dependence of massloss and windmomentum rate on metallicity. Powerlaw fits give average values and , however the fitting also reveals a somewhat steeper dependence for the lowerluminosity stars in our sample. The final fit relations presented in the previous section (Eqs. (19)–(20)) take this into account and give for the stars in our sample with luminosities above the mean (log(L_{*}∕10^{6} L_{⊙}) > −0.7) and for the ones below this mean (log(L_{*}∕10^{6} L_{⊙}) < −0.7). Again this is generally consistent with Puls et al. (2000), who derive a scaling relation . Inserting in this relation our Galactic value α_{eff} ≈ 0.48 and a typical Ostar δ ≈ 0.1 (see previous sections) gives , whereas using the lower α_{eff} ≈ 0.42 derived for the SMC yields . These values are in rather good agreement with the slopes we find from the scaling relations derived directly from the modelgrid results. This tentatively suggests that simplified linestatistics calculations such as those in Puls et al. (2000) might perhaps be used towards further calibration (and physical understanding) of the scaling relations derived in this paper, provided accurate values of α_{eff} (and , which is the line strength normalisation factor due to Gayley 1995) can be extracted from full hydrodynamic models.
As for the dependence of the terminal wind speed on metallicity, the exponent seems to vary across the grid, being positive for lowluminosity stars and negative for higher luminosity stars. This might be a manifestation of two (or more) competing processes. We already found that Ṁ always increases when increasing metallicity, which means that winds of high Z_{*} tend to be denser and thus harder to accelerate to high speeds. On the other hand, a higher metallicity also means higher abundances of important driving lines, providing a stronger acceleration that should increase the speed. For the low luminosity dwarfs, the second effect might be more dominant because these stars already have a low massloss rate. As a mean value we find that v_{∞} depends on the metallicity as . Previous studies such as Leitherer et al. (1992) and Krtička (2006) find this dependence to be and respectively.While these exponents are positive, all dependencies are very shallow.
As mentioned in previous sections, we do not find any strong trends with spectral type and luminosity class within our modelgrid. The massloss rates and windmomenta follow almost the same relations for spectroscopic dwarfs, giants, and supergiants. The dwarfs do show some deviation from this general trend, in particular for the SMC models, but this is mainly due to the fact that they span a much larger stellarparameter range than the other spectral classes, thus reaching lower luminosities and so also the low massloss rates where the SMC WLR starts to exhibit significant curvature (see above).
The dominant dependence of Ṁ on just L_{*} and Z_{*} within our Ostar grid may seem a bit surprising, especially in view of CAK theory which predicts an additional dependence (see Eq. (13)). However, including an additional massdependence in our powerlaw fits to the grid does not significantly improve the fitquality. Moreover, a multiregression fit to for the Galactic sample shows that if the massluminosity relation within our grid is also accounted for, the derived individual exponents return the previously found (which means that we get the same result as when previously neglecting any massdependence). Figure 9 shows this massluminosity relation, where a simple linear fit yields for our model grid. We note that these results are generally consistent with the alternative CMF models by Krtička & Kubát (2018), who also find no clear dependence on mass or spectral type within their Ostar grid.
On the other hand, when running some additional testmodels outside the range of our Ostar grid, by assuming a fixed luminosity but changing the mass, we do find that including an additional massdependence sometimes can improve the fits. Specifically, it is clear that reducing the mass for such constantluminosity models generally tends to increase Ṁ. This is in qualitative agreement with the results expected from CAK theory (see above), however, we note that these testmodels have stellar parameters that no longer fall within the Ostar regime that is the focus of this paper. In any case, when extending our grid towards a more full coverage of massive stars in various phases the question regarding an explicit massdependency will need to be revisited.
A notable feature for the SMC models is the bump that occurs at for the supergiants. Figure 10 shows a zoomin on the simulations that comprise this bump. For each of these models we examined the dominant ionisation stage at the sonic point of a selected list of important winddriving elements (Fe, C, N, O, and Ar). As we note above, the massloss rate in our simulations is most sensitive to the conditions around the sonic point. If the ionisation stage of an important element changes there, the opacity and thus also the radiative acceleration may also be modified. In Fig. 10 the models have increasing effective temperature and luminosity from left to right. The simple scaling relation of Eq. (20) would then predict increasing massloss rates from left to right, in contrast to what is observed in some of these models. As illustrated in the figure, this irregularity coincides with the temperatures at which several key drivingelements change their ionisation stages, which seem to produce a highly nonlinear effect over a restricted range in the current modelgrid. Although we have not attempted to explicitly account for such nonlinear temperature/ionisationeffects in the fittingrelations derived in this paper, this will certainly be important to consider when extending our grids towards lower temperatures (in particular the regime where Fe IV recombines to Fe III, where this effect may become very important, e.g. Vink et al. 1999).
Fig. 9 Mass versus luminosity (from Martins et al. 2005) of all models, both on logarithmic scale to present a linear relation. 
Fig. 10 Zoom in on the ‘bump’ that is notably present for the supergiants at SMC metallicity plotted as massloss rate versus luminosity on the lower axis and effective temperature on the upper axis. The dashed line shows the reference slope of 2.37 derived from all SMC models. Each model shows the dominant ionisation stage at the critical point of five elements. They are connected when a transition occurs between models. 
4.2 Comparison to other models
A key result from the overall grid analysis is that the computed massloss rates are significantly and systematically lower than those predicted by Vink et al. (2000, 2001), which are the ones most commonly used in applications such as stellar evolution and feedback. Figure 5 compares the Galactic modified windmomenta and massloss rates of our models directly to these Vink et al. predictions. Their theoretical WLR presented in the top panel is constructed from a fit to the objects of their sample of observed stars. This sample is used to ensure realistic parameters for v_{∞} and R_{*}, which together with the theoretical mass loss obtained from the recipe for Ṁ make up the modified windmomentum rate (see Vink et al. 2000, their Sect. 6.2). The original predictions are calculated with a higher value for metallicity of Z_{⊙} = 0.019, but this is scaled down here to our value of Z_{⊙} to remove the metallicity effect in this comparison. We note, however, that such a simple scaling actually overestimates the effect somewhat, since the abundance of important driving elements such as iron has not changed much (see Paper I). Although both set of models follow rather tight powerlaws with similar exponents (x = 2.06 vs. x = 1.83), there is a very clear offset of about 0.5 dex between the two modelrelations over the entire luminosity range. Indeed, every single one of our models lies significantly below the corresponding one by Vink et al. This is further illustrated in Fig. 11, showing a direct comparison between the massloss rates from the full grid (that is to say for all metallicities) and those computed by means of the Vink et al. recipe. In addition to a systematic offset, this figure also highlights the very low massloss rates we find for lowluminosity stars, in particular in the SMC (reflected in the final fitrelations, Eqs. (19)–20, through the steeper metallicity dependence at low luminosities). Regarding the metallicity scaling, the obtained byVink et al. (2001) (neglecting then any additional dependence on v_{∞} ∕v_{esc}) is in quite good agreement with the overall values derived here. As discussed in Paper I, the systematic discrepancies between our models and those by Vink et al. are most likely related to the fact that their massloss predictions are obtained on the basis of a global energy balance argument using Sobolev theory (and a somewhat different NLTE solution technique) to compute the radiative acceleration for a prespecified β velocity law. By contrast, the models presented here use CMF transfer to solve the full equation of motion and to obtain a locally consistent hydrodynamic structure (which may strongly deviate from a βlaw in those regions where Ṁ is initiated).
The bottom panel in Fig. 5 further compares our galactic results to those by Krtička & Kubát (2017, 2018). These authors also make use of CMF radiative transfer in the calculation of g_{rad} and they also find lower massloss rates as compared to Vink et al. However, the dependence on luminosity derived by Krtička & Kubát is weaker than that obtained here; this occurs primarily because they predict higher massloss rates for stars with lower luminosities, which tends to flatten their overall slope (see Fig. 5 and also Paper I). Although there are some important differences between the modelling techniques in Krtička & Kubát (2017) and those applied here^{3}, the overall systematic reductions of Ṁ found both here and by them may point towards the need for a reconsideration of the massloss rates normally applied in simulations of the evolution of massive Ostars (as further discussed in Sect. 4.6 below).
A third comparison of our results, now with those from Lucy (2010), is shown in the bottom panel of Fig. 5. Using a theory of moving reverse layers (MRL, Lucy & Solomon 1970), they compute the mass flux J = ρv for a grid of Ostars with different T_{eff} and log g. The MRL method assumes a plane parallel atmosphere and therefore does not yield a spherical massloss rate Ṁ directly. As such, we obtained Ṁ for the Lucy (2010) simulations by simply assuming the stellar parameters of the models in our own grid, extracting the massloss rate from . The relation shown in the bottom panel of Fig. 5 is then computed by performing a linear fit through the resulting massloss rates. This relation derived from the Lucy (2010) models also falls well below the Vink et al. curve, even predicting somewhat lower rates than us at low luminosities. We note that in the MRL method again no Sobolev theory is used for the MonteCarlo calculations determining g_{rad} and the mass flux.
Fig. 11 Direct comparison of the massloss rates calculated in this work versus those predicted by the Vink et al. recipe. The dashed line denotes the onetoone correspondence. The different markers show the luminosity classes consistent with previous plots. 
4.3 Comparison to observations
Observational studies aim to obtain empirical massloss rates based on spectral diagnostics in a variety of wavebands, ranging from the radio domain, over the IR, optical, UV, and all the way to highenergy Xrays. As discussed in Paper I (see also Puls et al. 2008; Sundqvist et al. 2011, for reviews), a key uncertainty in such diagnostics regards the effects of a clumped stellar wind on the inferred massloss rate. Indeed, if neglected in the analysis, such wind clumping may lead to empirical massloss rates that differ by very large factors for the same star, depending on which diagnostic is used to estimate this Ṁ (Fullerton et al. 2006). More specifically, Ṁ inferred from diagnostics depending on the square of density (for instance radio emission and H_{α}) are typically overestimated if clumping is neglected in the analysis. On the other hand, if porosity in velocity space (see Sundqvist & Puls 2018) is neglected, this may cause underestimations of rates obtained from UV line diagnostics (Oskinova et al. 2007; Sundqvist et al. 2010; Šurlan et al. 2013). Finally, regarding Ṁ determinations based on absorption of highenergy Xrays (Cohen et al. 2014), these have been shown to be relatively insensitive to the effects of clumping and porosity for Galactic Ostars (Leutenegger et al. 2013; Hervé et al. 2013; Sundqvist & Puls 2018).
4.3.1 The Galaxy
Considering the above issues, Fig. 12 compares the predictions from this paper with a selected sample of observational studies of Galactic Ostar winds. The selected studies are based on Xray diagnostics (Cohen et al. 2014), UV+optical analyses accounting for the effects of velocityspace porosity (Sundqvist et al. 2011; Šurlan et al. 2013; Shenar et al. 2015), and UV+optical+IR (Najarro et al. 2011) and UV+optical (Bouret et al. 2012) studies accounting for optically thin clumping (but neglecting the effects of velocityspace porosity^{4}). The figure shows the modified windmomentum rate from the selected observations (using different markers) together with our relation Eq. (19) (the solid line). Simple visual inspection of this figure clearly shows that our newly calculated models seem to match these observations quite well. More quantitatively, the dashed line is a fit through the observed windmomenta, excluding the two stars with L_{*} ∕L_{⊙} < 10^{5} but otherwise placing equal weights for all observational data points (also for those stars that are included in several of the chosen studies). For the observed stars with luminosities L_{*}∕L_{⊙} > 10^{5}, there is an excellent agreement with the theoretically derived relation Eq. (19), although there is also a significant scatter present in the data (quite naturally considering the difficulties in obtaining empirical Ṁ, see above).
Regarding the two excluded lowluminosity stars, these seem to indicate a weak wind effect also at Galactic metallicity. Also observational studies of Galactic dwarfs (Martins et al. 2005; Marcolino et al. 2009) and giants (de Almeida et al. 2019) find that the onset of the weak wind problem seems to lie around log (L_{*}∕L_{⊙}) ≈ 5.2, where such an onset was already indicated in the data provided by Puls et al. (1996). Although our massloss rates in this regime are significantly lower than those of Vink et al. (2000), they are still not as low as those derived from these observational studies. Further work would be required extending our current Milky Way grid to even lower luminosities, in order to examine if such an extension would yield a similar downward curvature in the WLR (albeit at a lower onset luminosity) as suggested by the observational data. Based on our results for SMC stars, we do expect that such a gridextension might eventually start to display a significant curvature in the WLR also for Galactic objects. Nonetheless, the mismatch in the onsetluminosity between theory and observations would then still need to be explained. On the other hand, further work is also needed to confirm the very low empirical massloss rates at these low luminosities, as none of the UVbased analyses cited above account for velocityporosity.
Fig. 12 Observed windmomenta for stars in the Galaxy from the studies discussed in the text are shown by different markers. The solid black line is our derived relation (19) and the dashed line is a fit through the observations excluding the data points with L_{*} ∕L_{⊙} < 10^{5} (see text). 
Fig. 13 Observations of SMC stars from Bouret et al. (2013) are shown with black squares. Both the relation as predicted by Vink et al. (2001) and from Eq. (20) are shown as comparison together with a linear fit through the observational data. 
4.3.2 Low metallicity
The comparison above focused exclusively on Galactic Ostars, since similar analyses including adequate corrections for clumping are, to this date, scarce for Magellanic Cloud stars. However, a few such studies do exist, for example the one performed by Bouret et al. (2013) for O dwarfs in the SMC. Again, porosity in velocity space has not been accounted for in deriving Ṁ from the observationsin this study. Nonetheless, a comparison of our predicted rates to the empirical ones derived by Bouret et al. is shown in Fig. 13. The figure displays our fit relation found for the SMC stars together with the results from Bouret et al. (2013). Additionally the figure shows a fit to this observed sample of stars (dashdotted line) and again the relation for the samestars as would be predicted by the Vink et al. recipe. It is directly clear from this figure that, again, our relation matches these observations of SMC O dwarfs significantly better than previous theoretical predictions. Indeed, since no big systematic discrepancies are found here between our predictions and the observations, the previously discussed mismatch for lowluminosity Galactic Odwarfs (“weak wind problem”) does not seem to be present here for SMC conditions. This effect is also reflected through the significantly higher overall WLR slope we find for SMC models as compared to Galactic ones, x_{SMC} = 2.56 ± 0.44 versus x_{Gal} = 2.07 ± 0.32. In terms of linestatistics, a quite natural explanation for this behaviour is that α_{eff} decreases with wind density (see previous sections) which produces a steeper slope both for decreasing metallicity and decreasing luminosity. Another consequence of the found agreement could be that velocity porosity effects are negligible in the observations of SMC stars.
Moreover, we may also (in a relative manner) compare our predicted scaling relations to observations in different galactic environments, as long as we assume that clumping properties do not vary significantly between the considered galaxies (and so does not significantly affect the relative observational results). In a large compilationstudy of Ostars in the Galaxy, LMC, and SMC, Mokiem et al. (2007) obtain observationally inferred values of D_{mom}. Specifically, Mokiem et al. (2007) used the H_{α} emission line to derive the massloss rate. This diagnostic is highly dependent on the clumping factor, but comparison with our scaling results will still be reasonable if the clumping in the H_{α} forming region does not vary too much between the different metallicity environments. Neglecting such potential clumping effects, Mokiem et al. derived the empirical WLR slopes x_{Gal} = 1.86 ± 0.2, x_{LMC} = 1.87 ± 0.19, and x_{SMC} = 2.00 ± 0.27. Within the 1σ errors, this is in agreement with the theoretical values x_{Gal} = 2.07 ± 0.32, x_{LMC} = 2.12 ± 0.34, and x_{SMC} = 2.56 ± 0.44 obtained here (although the general trend is that we find somewhat steeper relations). In this respect, we note also that the tight 1σ limits in the empirical relations obtained by Mokiem et al. (2007) only reflect their fiterrors and not additional systematic errors due to uncertainties in, such as stellar luminosities and adopted metallicities. As such, the overall agreement between their empirical relations and the theoretical WLR slopes presented in this paper is encouraging. Moreover, by considering their results at a fixed luminosity log L_{*}∕L_{⊙} = 5.75, Mokiem et al. (2007) also derived an empirical massloss vs. metallicity relation ; again this is in rather good agreement with the theoretical obtained at such a logL_{*}∕L_{⊙} = 5.75 from Eq. (20) of this paper.
4.3.3 Terminal wind speed
As can be seen in Fig. 6, v_{∞}∕v_{esc} ratios range between approximately 2.5 to 5.5 for our Galactic models. Similar values are found for the Magellanic Cloud simulations, but the scatter is significant and it is difficult to draw general conclusions from the model data. A significant scatter in v_{∞}∕v_{esc} is however quite consistent with observational compilations such as the one presented by Garcia et al. (2014; see in particular their Fig. 9) and Lamers et al. (1995). Overall though, our mean values are somewhat higher than observations generally suggest (see also Kudritzki & Puls 2000), in particular for lowluminosity objects in the Galaxy. This issue was also addressed in Paper I, where it was argued that the high v_{∞} for such lowdensity winds might be naturally reduced as a consequence of inefficient cooling of shocks in the supersonic wind (see also Lucy 2012). Namely, if a significant portion of the supersonic wind were shockheated, this would lower g_{rad} in these layers and so also potentially reduce v_{∞}. Moreover, if the wind is shockheated, empirical estimates of v_{∞} might be misinterpreted as it may no longer be possible to accurately derive this quantity from observations of UV wind lines. Tailored radiationhydrodynamical simulations of linedriven instability (LDI) induced shocks in lowdensity winds are underway in order to examine this in detail and will be presented in an upcoming paper (Lagae et al., in prep.). For now we study how a reduced wind driving in supersonic layers may affect v_{∞} and Ṁ in our steadystate simulations by running a few additional models where g_{rad} in layers with speeds above 1000 km s^{−1} simply is reduced by an adhoc factor of two. While these models then indeed converged to a significantly lower v_{∞} than previously, Ṁ was only marginally affected. This suggests that (at least within the steadystate simulation setup used here) a reduction of the outer wind driving does not necessarily lead to an increased massloss rate, because in our simulations Ṁ is mostly sensitive to the conditions at the sonic point and is not as much affected by the reduction of g_{rad} in the highly supersonic regions. Although shockheating in the outer winds in lowluminosity objects might explain the high predicted values of the low luminosity objects, it is unclear to what degree the increasing trend in v_{∞} /v_{esc} for decreasing L_{*} might still persist in such models. Similarly, it is also somewhat unclear how this matches with observations, as Lamers et al. (1995) do not find such a trend while Garcia et al. (2014) do show a similar trend, at least for dwarf stars. All these observations show significant scatter though, such that a clear prediction is difficult to obtain from the present data. To a major part, the scatter in v_{∞} is understandable though, since in the outer wind only a few dozen strong resonance lines are responsible for the acceleration (mostly from C, N, O, Ne, and Ar). Thus, with only a few lines, a certain difference in the abundances (for example because of different ages and mixing) can have a significant effect on v_{∞}, leading to the observed scatter.
4.4 Influence of v_{turb}
As discussed in Paper I, the turbulent velocity v_{turb} can have a significant effect on g_{rad} around the sonic point and so also affect the predicted Ṁ. In general, increasing v_{turb} tends to decrease Ṁ (Lucy 2010; Krtička & Kubát 2017; Paper I), mostly because then higher velocities are required to Doppler shift line profiles out of their own absorption shadows, reducing g_{rad} in the critical layers. To study the behaviour of our grid with v_{turb}, we ran a set of additional simulations where, for each of our previous models, we increased v_{turb} to 12.5 km s^{−1} and decreased it to 7.5 km s^{−1}. This way a slope of logṀ with logv_{turb} could be found for each model in the grid. The mean value obtained is (21)
where the error is derived from the 1σ spread of the slope for all models. Also Lucy (2010) studied the behaviour of the mass flux J with v_{turb}. Similar to this study, he finds an inverse dependence, but with a somewhat steeper slope −1.46. We note, however, that this slope was derived from a single model at T_{eff} = 40 000 K and log g = 3.75 whereas we here consider an average across our full grid. Indeed, inspection of a single model similar to the one used in Lucy (2010), with T_{eff} = 40 062 K and log g = 3.92, reveals a slope −1.49 that is in good agreement with their result. Furthermore, considering the large scatter on the slope that we find, with values between a maximum of 0.07 and a minimum of −1.88, the result ofEq. (21) is also still in agreement with that of Lucy (2010). But although the scatter around our derived meanslope thus is significant, there are no clear trends within the grid. As such, to obtain a firstorder approximation accounting for a v_{turb} that deviates from the standard value 10 km s^{−1}, the Ṁ relations in Sect. 3.3 may simply be scaled according to Eq. (21).
4.5 Comparison to β velocity law
Most wind models of hot, massive stars used for spectroscopic studies actually do not solve the e.o.m. for the outflow. Instead, these models assume an empirical ‘βtype’ wind velocitylaw that connects to a quasihydrostatic photosphere. This is the case also for the “standard” version of FASTWIND, used here as a starting condition for our selfconsistent simulations (see Sect. 2). For the prototypical simulation presented in Sect. 2.2, Fig. 4 illustrates a fit to the selfconsistent calculated velocity using such a βlaw. More specifically, a “double” βlaw similar to that presented in Sect. 5.2 of Paper I is used, however in order to obtain a better fit to the very steep acceleration of the transition region we have also included a modificationterm in dependence on the photospheric scaleheight. The full expression used to match the velocity structure is (22)
In this relation, the parameters v_{∞}, β_{1}, β_{2}, v_{exp}, and H are obtained by fitting to the numerically derived velocity structure above the transition radius r_{tr}, defined according to v(r_{tr}) = v_{tr}. Introduced here is the parameter v_{exp}, which is roughly the velocity at which v(r) has its biggest curvature and thus controls how far the exponential behaviour of the velocity holds in the inner wind. Also introduced is H, which is thescale height setting the stratification in the photosphere, controlled by the density structure of the model.
With this, the model displayed in Fig. 4 can be fitted down to a transition velocity v_{tr} ≈ 0.1a, for a primary betafactor β_{1} = 0.8. We note though, that in particular the more lowdensity winds in our sample typically require much higher transition velocities in order to be well fit, often values well above v_{tr} ≈ 0.5a are found (see also discussion in Paper I). Moreover, most of the bestfit β_{1} values lie between 0.5–0.8, indicating a steep acceleration of the inner wind across our grid. A more detailed study aiming to develop an improved velocityparametrisation for spectroscopic studies is underway and will be presented in an upcoming paper. There, we also plan to examine in detail what effects the predicted steep acceleration in the transition region may have upon the formation of various strategic spectral lines used for diagnostic work.
4.6 Implications for stellar evolution
The stellar mass is the most important parameter defining the evolution of a star and accurate massloss rates are crucial to determine the corresponding evolutionary pathways. Codes for stellar structure and evolution use prescribed recipes calculating the change in stellar mass in between time steps. Many codes, such as MESA (Paxton et al. 2019, and references therein), offer the Vink et al. recipe as an option to calculate Ṁ for hot, hydrogenrich stars (excluding then the classical WolfRayet (WR) stars which require different prescriptions for Ṁ). As clearly illustrated by Fig. 11, the results presented here suggest that the Ostar massloss rates should be significantly lower. In addition, this is also supported by the comparison to observations in Sect. 4.3. Even though the calculations here are performed for massive stars in their early phases of evolution (mostly on the mainsequence), the lower rates will not only affect the stars there, but also impact the properties of postmainsequence stages. One consequence is that the luminosity at which the stars end their mainsequence evolution and cross the Hertzprung gap is changed, as seen directly in our evolution calculations of a 60 M_{⊙} star using MESA (where we have simply reduced the amount of mass loss on the main sequence to be in accordance with the models presented here). Another effect is that a lower mass loss means that angular momentum is lost less rapidly so that the star keeps a higher surface rotation. Keszthelyi et al. (2017) computed models with reduced massloss rates (factor 2 to 3) and find the surface rotation speeds at the end of the mainsequence to remain rather high, possibly requiring an additional source of angular momentum loss to reconcile with observed values.
Moreover, Belczynski et al. (2010) studied the masses of compact objects originating from single stars. Here, the adopted mass loss during the life of the progenitor star is crucial in determining the maximum black hole mass that can be achieved. The determinations of the black hole masses in Belczynski et al. (2010) were done assuming the Vink et al. rates; adopting lower rates would reduce the amount of mass lost and thus possibly increase the resulting black hole mass. Getting direct measurements of these black hole masses is not straightforward. However, by taking advantage of gravitational wave astronomy, detections such as GW150914 (Abbott et al. 2016) provide observational constraints. Derived blackhole masses turn out to be relatively high, ≿ 25 M_{⊙}, which, with current wind prescriptions, can only be created in low metallicity environments. With lower values of the massloss rate during the evolution of the star, such as proposed here in this paper (or as in magnetic massive stars, Petit et al. 2017), the ‘heavy’ black holes as detected by gravitational waves might in principle also be produced in highmetallicity environments (Belczynski et al. 2020). So far though, in the Galaxy, no observationsof such heavy mass black holes (from single stars) exist as all of them have a mass less than 15 M_{⊙} (as found in studies such as Shaposhnikov & Titarchuk 2007; Torres et al. 2020). Depending on the initial mass of the progenitor stars, these values can be explained both by our proposed massloss rates as well as with those predicted by Vink et al. (2001).
Significant mass loss is further necessary to create the naked Helium core that is a classical WRstar through windstripping. In low metallicity environments such as the SMC, this is generally considered to be difficult because of the strong metallicity dependence resulting in low massloss rates. Considering our results here of lower Ostar massloss rates, such (steady) windstripping would be more difficult to obtain also in the Galaxy, potentially leading to an increase of the lower limit for the initial mass of WRstars created by this channel. In order to explain the observed number of WRstars, alternative pathways might thus be necessary. In this respect, a straightforward option is binarity where the outer layers of stars can be removed through binary interaction such as Rochelobe overflow (e.g. Götberg et al. 2018). Recent studies have shown that a large majority of massive stars reside in such binary systems (Sana et al. 2014). On the other hand, a second possible channel is eruptive mass loss in the luminous blue variable stage (LBV; Smith 2014). Indeed, significant fractions of the stellar mass can be removed in such eruptive events; the LBV ηCarina, for example, lost 10 M_{⊙} in just 10 yr in the 19th century. Also considering that the (so far confirmed) binary fraction of WRstars in the SMC is lower than that of the Galaxy, or at least similar (Foellmi et al. 2003), this latter pathway might prove to be of increased importance.
5 Summary and future prospects
We calculated a grid of steadystate wind models of Ostars by varying fundamental stellar parameters in three metallicity regimes corresponding to the Galaxy, the Small, and the Large Magellanic Clouds. The models provide predictions of global wind parameters such as massloss rate and windmomentum rate, allowing us to analyse how these quantities depend on fundamental stellar parameters such as luminosity and metallicity.
From our grid, we find steep dependencies of the massloss rate with both luminosity and metallicity, with mean values and . The metallicity dependence is further found to vary across the luminosity range. Accounting for this, results in the final fit relations for the windmomentum rate and massloss rate presented in Sect. 3.3. Additionally, a clear change in the slope for the predicted WLR for dwarfs in the SMC is found, pointing towards the occurrence of weak winds in the models.
Our computed massloss rates are significantly lower for all models than those predicted by Vink et al. (2000, 2001), which are the ones usually implemented in evolution calculations of massive stars. Such lower Ostar windmomenta and massloss rates are also in general accordance with observational studies in the Galaxy that properly account for the effects of clumping upon the diagnostics used to infer the empirical massloss rates. Regarding the metallicity dependence, our scaling predictions are (within the errors) in agreement with the larger empirical study by Mokiem et al. (2007).
The systematically reduced massloss rates for all models strengthens the claim that new rates might be needed in evolution simulations of massive stars. Namely, adopting different rates can significantly affect the evolution of the massive star, for example by changing its spindown time and altering the initial mass needed in order to produce a windstripped WolfRayet star. As such, a key followup work to this study will be to now extend the grid presented here to include massive stars outside the Ostar domain, to incorporate our new models into simulations of massivestar evolution, and to analyse in detail the corresponding effects.
Acknowledgements
R.B. and J.O.S. acknowledge support from the Odysseus programme of the Belgian Research Foundation Flanders (FWO) under grant G0H9218N. J.O.S. also acknowledges support from the KU Leuven C1 grant MAESTRO C16/17/007. F.N. acknowledges financial support through Spanish grants ESP201786582C41R and PID2019105552RBC41 (MINECO/MCIU/AEI/FEDER) and from the Spanish State Research Agency (AEI) through the Unidad de Excelencia “María de Maeztu”Centro de Astrobiología (CSICINTA) project No. MDM20170737. We would also like to thank the referee for useful comments that led to additional improvements of the paper.
Appendix A Model parameters
Input parameters of all models in the grid together with the resulting massloss rate and terminal wind speed.
References
 Abbott, B. P., Abbott, R., Abbott, T. D., et al. 2016, ApJ, 818, L22 [Google Scholar]
 Asplund, M., Grevesse, N., Sauval, A. J., & Scott, P. 2009, ARA&A, 47, 481 [NASA ADS] [CrossRef] [Google Scholar]
 Belczynski, K., Bulik, T., Fryer, C. L., et al. 2010, ApJ, 714, 1217 [Google Scholar]
 Belczynski, K., Hirschi, R., Kaiser, E. A., et al. 2020, ApJ, 890, 113 [Google Scholar]
 Bouret, J.C., Hillier, D. J., Lanz, T., & Fullerton, A. W. 2012, A&A, 544, A67 [Google Scholar]
 Bouret, J. C., Lanz, T., Martins, F., et al. 2013, A&A, 555, A1 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Bresolin, F., Crowther, P. A., & Puls, J., eds. 2008, IAU Symposium, 250, Massive Stars as Cosmic Engines [Google Scholar]
 Castor, J. I., Abbott, D. C., & Klein, R. I. 1975, ApJ, 195, 157 [NASA ADS] [CrossRef] [Google Scholar]
 Cohen, D. H., Wollman, E. E., Leutenegger, M. A., et al. 2014, MNRAS, 439, 908 [Google Scholar]
 de Almeida, E. S. G., Marcolino, W. L. F., Bouret, J. C., & Pereira, C. B. 2019, A&A, 628, A36 [Google Scholar]
 Driessen, F. A., Sundqvist, J. O., & Kee, N. D. 2019, A&A, 631, A172 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Foellmi, C., Moffat, A. F. J., & Guerrero, M. A. 2003, MNRAS, 338, 360 [Google Scholar]
 Fullerton, A. W., Massa, D. L., & Prinja, R. K. 2006, ApJ, 637, 1025 [NASA ADS] [CrossRef] [Google Scholar]
 Garcia, M., Herrero, A., Najarro, F., Lennon, D. J., & Alejandro Urbaneja, M. 2014, ApJ, 788, 64 [Google Scholar]
 Garcia, M., Herrero, A., Najarro, F., Camacho, I., & Lorenzo, M. 2019, MNRAS, 484, 422 [Google Scholar]
 Gayley, K. G. 1995, ApJ, 454, 410 [NASA ADS] [CrossRef] [Google Scholar]
 Götberg, Y., de Mink, S. E., Groh, J. H., et al. 2018, A&A, 615, A78 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Hamann, W.R. 1981, A&A, 93, 353 [NASA ADS] [Google Scholar]
 Hervé, A., Rauw, G., & Nazé, Y. 2013, A&A, 551, A83 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Hillier, D. J., & Miller, D. L. 1998, ApJ, 496, 407 [NASA ADS] [CrossRef] [Google Scholar]
 Keszthelyi, Z., Puls, J., & Wade, G. A. 2017, A&A, 598, A4 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Krtička J. 2006, MNRAS, 367, 1282 [NASA ADS] [CrossRef] [Google Scholar]
 Krtička, J., & Kubát, J. 2017, A&A, 606, A31 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Krtička, J., & Kubát, J. 2018, A&A, 612, A20 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Kudritzki, R.P., & Puls, J. 2000, ARA&A, 38, 613 [NASA ADS] [CrossRef] [Google Scholar]
 Kudritzki, R. P., Pauldrach, A., & Puls, J. 1987, A&A, 173, 293 [NASA ADS] [Google Scholar]
 Kudritzki, R.P., Lennon, D. J., & Puls, J. 1995, in Science with the VLT, eds. J. R. Walsh, & I. J. Danziger, 246 [Google Scholar]
 Lamers, H. J. G. L. M., CerrutiSola, M., & Perinotto, M. 1987, ApJ, 314, 726 [Google Scholar]
 Lamers, H. J. G. L. M., Snow, T. P., & Lindholm, D. M. 1995, ApJ, 455, 269 [NASA ADS] [CrossRef] [Google Scholar]
 Leitherer, C., Robert, C., & Drissen, L. 1992, ApJ, 401, 596 [NASA ADS] [CrossRef] [Google Scholar]
 Leutenegger, M. A., Cohen, D. H., Sundqvist, J. O., & Owocki, S. P. 2013, ApJ, 770, 80 [Google Scholar]
 Lucy, L. B. 1971, ApJ, 163, 95 [NASA ADS] [CrossRef] [Google Scholar]
 Lucy, L. B. 2010, A&A, 524, A41 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Lucy, L. B. 2012, A&A, 544, A120 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Lucy, L. B., & Solomon, P. M. 1970, ApJ, 159, 879 [NASA ADS] [CrossRef] [Google Scholar]
 Marcolino, W. L. F., Bouret, J., Martins, F., et al. 2009, A&A, 498, 837 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Martins, F., Schaerer, D., & Hillier, D. J. 2005, A&A, 436, 1049 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Mokiem, M. R., de Koter, A., Vink, J. S., et al. 2007, A&A, 473, 603 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Muijres, L. E., de Koter, A., Vink, J. S., et al. 2011, A&A, 526, A32 [Google Scholar]
 Najarro, F., Hanson, M. M., & Puls, J. 2011, A&A, 535, A32 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Oskinova, L. M., Hamann, W.R., & Feldmeier, A. 2007, A&A, 476, 1331 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Owocki, S. P., Castor, J. I., & Rybicki, G. B. 1988, ApJ, 335, 914 [Google Scholar]
 Pauldrach, A. W. A., Hoffmann, T. L., & Lennon, M. 2001, A&A, 375, 161 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Paxton, B., Smolec, R., Schwab, J., et al. 2019, ApJS, 243, 10 [Google Scholar]
 Petit, V., Keszthelyi, Z., MacInnis, R., et al. 2017, MNRAS, 466, 1052 [Google Scholar]
 Puls, J. 2017, in IAU Symposium, 329, The Lives and DeathThroes of Massive Stars, eds. J. J. Eldridge, J. C. Bray, L. A. S. McClelland, & L. Xiao, 435 [Google Scholar]
 Puls, J., Kudritzki, R.P., Herrero, A., et al. 1996, A&A, 305, 171 [NASA ADS] [Google Scholar]
 Puls, J., Springmann, U., & Lennon, M. 2000, A&AS, 141, 23 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Puls, J., Urbaneja, M. A., Venero, R., et al. 2005, A&A, 435, 669 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Puls, J., Vink, J. S., & Najarro, F. 2008, A&ARv,, 16, 209 [Google Scholar]
 Puls, J., Najaro, F., Sundqvist, J., & Sen, K. 2020, A&A, 642, A172 [EDP Sciences] [Google Scholar]
 Sana, H., Le Bouquin, J. B., Lacour, S., et al. 2014, ApJS, 215, 15 [Google Scholar]
 Sander, A. A. C., Hamann, W.R., Todt, H., Hainich, R., & Shenar, T. 2017, A&A, 603, A86 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Shaposhnikov,N., & Titarchuk, L. 2007, ApJ, 663, 445 [Google Scholar]
 Shenar, T., Oskinova, L., Hamann, W. R., et al. 2015, ApJ, 809, 135 [Google Scholar]
 Smith, N. 2014, ARA&A, 52, 487 [Google Scholar]
 Sundqvist, J. O., & Puls, J. 2018, A&A, 619, A59 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Sundqvist, J. O., Puls, J., & Feldmeier, A. 2010, A&A, 510, 11 [Google Scholar]
 Sundqvist, J. O., Puls, J., Feldmeier, A., & Owocki, S. P. 2011, A&A, 528, 64 [Google Scholar]
 Sundqvist, J. O., Owocki, S. P., & Puls, J. 2018, A&A, 611, A17 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Sundqvist, J. O., Björklund, R., Puls, J., & Najarro, F. 2019, A&A, 632, A126 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Šurlan, B., Hamann, W.R., Aret, A., et al. 2013, A&A, 559, A130 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Torres, M. A. P., Casares, J., JiménezIbarra, F., et al. 2020, ApJ, 893, L37 [Google Scholar]
 Vink, J. S., de Koter, A., & Lamers, H. J. G. L. M. 1999, A&A, 350, 181 [Google Scholar]
 Vink, J. S., de Koter, A., & Lamers, H. J. G. L. M. 2000, A&A, 362, 295 [Google Scholar]
 Vink, J. S., de Koter, A., & Lamers, H. J. G. L. M. 2001, A&A, 369, 574 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Relaxing this assumption would have an only marginal effect on the occupation numbers (e.g. Hamann 1981; Lamers et al. 1987). Also for the calculation of g_{rad}, pure Dopplerbroadening is sufficient, since (i) only few strong lines have significant natural and – in the photosphere – collisionally linebroadened wings that could contibute to g_{rad}; and (ii), because of lineoverlap effects, these wings are typically dominated by the (Dopplercore) line opacity from other transitions, which then dominate the acceleration.
Such a value results from considering the distribution of oscillator strengths for resonance lines within a hydrogenic ion and neglecting δ, see Puls et al. (2000).
In particular, since Krtička & Kubát (2017, 2018) scale their CMF line force to the corresponding Sobolev force, this means that their critical point is no longer the sonic point and so that the nature of their basic hydrodynamic steadystate solutions may be quite different from those presented here. See also discussion in Paper I.
Although these two last studies indeed do not account for velocityporosity, they do attempt to adjust their studies accordingly; Najarro et al. (2011) by analysis of IR lines that should be free of such effects and Bouret et al. (2012) by scaling down the phosphorus abundance, thus mimicking the effect velocity–porosity would have on the formation of the unsaturated UV PV lines. As such, we opt here to include also these two studies in our selected sample for observational comparisons.
All Tables
Input parameters of all models in the grid together with the resulting massloss rate and terminal wind speed.
All Figures
Fig. 1 Top panel: value of Γ versus scaled radiuscoordinate for 7 (nonconsecutive) hydrodynamic iterations over the complete run. The starting structure (yellow) relaxes to the final converged structure (dark blue). Bottom panel: colour map of log (f_{err}) for all hydrodynamic iterations; on the abscissa is hydrodynamic iteration number and on the ordinate scaled wind velocity. The pluses indicate the location of for each iteration, the dashed lines the limits between which is computedand the dashdotted line the location of the sonic point. 

In the text 
Fig. 2 Top panel: final converged structure of the characteristic model showing Γ in black squares and Λ (see text for definition) as a green line. The black dashed lines show the location of the sonic point approximately at Γ = 1 (but not exactly because of the additional pressure terms). Bottom panel: same as in the top panel, however versus velocity (which resolves the inner wind more). 

In the text 
Fig. 3 Top panel: iterative behaviour of the massloss rate as decreases towards a value below 1%. The colour signifies the iteration number starting from Ṁ as predicted by the Vink et al recipe in light green. Bottom panel: iterative behaviour of the terminal velocity v_{∞} towards convergence. 

In the text 
Fig. 4 Converged velocity structure for the characteristic model of Sect. 2.2, showing velocity over terminal wind speed versus the scaled radiuscoordinate in black. The green line shows a fit using a double βlaw following Eq. (22) (see text for details). 

In the text 
Fig. 5 Top panel: modified windmomentum rate versus luminosity for all Galactic models. The solid black line is a linear fit through the points and the dashed line shows the theoretical relation by Vink et al. (2000). Bottom panel: massloss rates versusluminosity for all Galactic models with a linear fit as a solid black line. The dashed line is a fit through the massloss rates computed using the Vink et al. recipe, the dashdotted line is the relation derived by Krtička & Kubát (2017) and the dotted line is the relation computed from the results from Lucy (2010). 

In the text 
Fig. 6 Terminal wind speed over photospheric effective escape speed (corrected with 1Γ_{e}, see text) versus luminosity. Values for Galactic metallicity for all three luminosity classes are shown. 

In the text 
Fig. 7 Top panel: modified windmomentum rate of all models versus luminosity. The dashed lines show linear fits through each of the three sets of models. The markers show the different luminosity classes, consistent with previous plots. Bottom panel: same as the top panel, but for the massloss rate. 

In the text 
Fig. 8 Value of the exponent n showing the metallicity dependence of D_{mom} together with a linear fit, showing the linear dependence of this exponent with log (L_{*} ∕10^{6} L_{⊙}). See text for details. 

In the text 
Fig. 9 Mass versus luminosity (from Martins et al. 2005) of all models, both on logarithmic scale to present a linear relation. 

In the text 
Fig. 10 Zoom in on the ‘bump’ that is notably present for the supergiants at SMC metallicity plotted as massloss rate versus luminosity on the lower axis and effective temperature on the upper axis. The dashed line shows the reference slope of 2.37 derived from all SMC models. Each model shows the dominant ionisation stage at the critical point of five elements. They are connected when a transition occurs between models. 

In the text 
Fig. 11 Direct comparison of the massloss rates calculated in this work versus those predicted by the Vink et al. recipe. The dashed line denotes the onetoone correspondence. The different markers show the luminosity classes consistent with previous plots. 

In the text 
Fig. 12 Observed windmomenta for stars in the Galaxy from the studies discussed in the text are shown by different markers. The solid black line is our derived relation (19) and the dashed line is a fit through the observations excluding the data points with L_{*} ∕L_{⊙} < 10^{5} (see text). 

In the text 
Fig. 13 Observations of SMC stars from Bouret et al. (2013) are shown with black squares. Both the relation as predicted by Vink et al. (2001) and from Eq. (20) are shown as comparison together with a linear fit through the observational data. 

In the text 
Current usage metrics show cumulative count of Article Views (fulltext article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.
Data correspond to usage on the plateform after 2015. The current usage metrics is available 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.