© 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 radiation-driven 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⁻¹, 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 mass-loss rates and wind-momenta 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 mass-loss predictions based on steady-state wind models using radiative transfer in the co-moving 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 O-stars 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 mass-loss 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 O-stars 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 massive-star evolution.

Since most of the important spectral lines driving hot-star winds are metallic, a strong dependence on metallicity Z_* is expected for the mass-loss 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 mass-loss 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 scaling-relations for wind-momenta (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 mass-loss rates and wind-momenta, we perform a study using a grid of O-star models. Thanks to the fast performance of the method, as explained in detail in Paper I, this is finally possible for the hydrodynamically consistent steady-state models with a non-paramterised CMF line-force computed without any assumptions about underlying line-distribution functions. In Sect. 2, we briefly review our method for computing mass-loss 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 wind-momenta and mass-loss 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 so-called weak wind problem. Section 5 contains the conclusions and future prospects.

2 Methods

A crucial part about the radiation-driven steady-state 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, steady-state 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² 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 radiation-driven wind. For a steady-state mass-loss rate Ṁ the mass conservation equation Ṁ = 4πr²ρv gives the density structure ρ(r). The wind modelsfurther rely on the NLTE (non-local thermodynamic equilibrium) radiative transfer in FASTWIND (see Paper I and Puls 2017) for the computation of g_rad, by means ofa co-moving frame (CMF) solution and without using any parametrised distribution functions. The atomic data are taken from the WM-BASIC 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 O-stars. 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 WM-BASIC 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 Doppler-broadening alone¹ 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 steady-state Ṁ and v(r) are converged in the model computation starting from a first guess. For all simulations presented in this paper, the start-value for Ṁ was taken as the mass-loss 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 quasi-hydrostatic atmosphere connects at v_tr ≈ 0.1a(T = T_eff) to a so-called β-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 flux-weighted optical depth (5)

for the flux-weighted opacity κ_F (cm² g⁻¹). The flux-weighted opacity is related to the radiative acceleration as g_rad = κ_FL_*∕(4πcr²). 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⁻² 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 CAK-theory cannot be used to update the mass-loss 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 mass-loss 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 mass-loss rate according to , following the basic theory of line-driven 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 steady-state mass-loss 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 radiative-transfer 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)

with (10)

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 mass-loss 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 dash-dotted 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⁻¹ there is some discrepancy; this arises because in these quasi-static layers the flux-weighted opacity is approximated by a Kramer-like 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 mass-loss 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 mass-loss rate computed for that iteration. The general trend is that decreases quite consistently during the iteration cycle. The mass-loss 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 anti-correlation, as for this model the total wind-momentum 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 inner-most radial point of the grid. This figure illustrates the very steep acceleration around the sonic point that is characteristic for our O-star models, followed by a supersonic β-law-like behaviour typical of a radiation-driven 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 radius-coordinate for 7 (non-consecutive) hydrodynamic iterations over the complete run. The starting structure (yellow) relaxes to the final converged structure (dark blue). Bottom panel: colour map of log (ferr) 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 dash-dotted 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 mass-loss rates of O stars in a quantitative way, a model-grid 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 non-LTE, line-blanketed synthetic spectra models using the code CMFGEN (Hillier & Miller 1998). Here, a self-consistent 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 mass-loss rate. The wind models were computed with the microturbulent velocity kept at a standard value for O-stars of v_turb = 10 km s⁻¹ (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 high-energy X-rays. 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 mass-loss rates (indeed, in the simple tests performed in Paper I the effect was only marginal). In addition, any inclusion of wind clumping into steady-state models will inevitably be of ad-hoc 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 mass-loss rate by as much as an order of magnitude. The models used in their study, however, use a global energy constraint to derive the mass-loss 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 mass-loss 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 mass-loss 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 wind-clumping, 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 steady-state 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 set-up has the advantage of enabling a rather direct comparison for the model-dependence 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 mass-loss 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 radius-coordinate 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 line-driven winds, it is useful to look at a modified wind-momentum rate as function of the stellar luminosity: (12)

where the left hand side is the so-called modified wind-momentum 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 wind-momentum 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² g⁻¹). 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² (15)

as in the wind-momentum-luminosity 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 line-driven 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 radiation-hydrodynamic wind models presented here is shown in the top panel of Fig. 5. The figure indeed shows a quite tight relationship between the modified wind-momentum 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 O-star values α ≈0.6 and δ ≈ 0.1 (Puls et al. 2008).

Next, the mass-loss 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 power-law here gives y = 2.16 ± 0.34. Within our grid, we further do not find any strong systematic trends of mass-loss rate (or wind-momentum 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 O-stars, 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 wind-momentum 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: mass-loss rates versusluminosity for all Galactic models with a linear fit as a solid black line. The dashed line is a fit through the mass-loss rates computed using the Vink et al. recipe, the dash-dotted 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 wind-momenta for the Magellanic Cloud simulations are added to those of the Galaxy, alsoincluding power-law 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 mass-loss and wind-momentum rates are always lower for the LMC than for the Galaxy and lower still for the SMC. This is as expected for radiation-driven O-star 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 low-luminosity 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 mass-loss 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 mass-loss rates found for the stars with the lowest luminosities. As discussed further in Sect. 4, these findings are generally consistent with the line-statistics 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 wind-momentum 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 mass-loss rate.

3.3 Function of metallicity

Examining trends of the modified wind-momentum 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 wind-momenta 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⁶ L_⊙), we find n(L_*) = −0.73log(L_*∕10⁶ 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⁶ L_⊙) + 0.79.

Building on the combined results above, we can now construct final relations for both the modified wind-momentum rate and mass-loss rate of the form: (18)

To obtain the fitting coefficients (which will be different for the wind-momenta and mass-loss 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 multi-linear 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)

and (20)

with the wind momentum rate in units of M_⊙ yr⁻¹ km s⁻¹ R_⊙^0.5 and the mass − loss rate in units of M_⊙ yr⁻¹. When the complete model-grid 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 wind-momentum rate and the mass-loss 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⁶ 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⁶ 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 Dmom together with a linear fit, showing the linear dependence of this exponent with log (L* ∕106 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 line-driven 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 line-statistics results (see overview in Puls et al. 2008). For the O-stars 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 model-grid 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 line-statistics calculations by Puls et al. (2000) and may (at least qualitatively) be understood via the physical interpretation of α as the ratio of the line-force due to optically thick lines to the total line-force (lower wind densities should generally mean a lower proportion of optically thick contributing lines).

We further find a quite steep dependence of mass-loss and wind-momentum rate on metallicity. Power-law fits give average values and , however the fitting also reveals a somewhat steeper dependence for the lower-luminosity 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⁶ L_⊙) > −0.7) and for the ones below this mean (log(L_*∕10⁶ 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 O-star δ ≈ 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 model-grid results. This tentatively suggests that simplified line-statistics 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 low-luminosity 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 mass-loss 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 model-grid. The mass-loss rates and wind-momenta 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 stellar-parameter range than the other spectral classes, thus reaching lower luminosities and so also the low mass-loss rates where the SMC WLR starts to exhibit significant curvature (see above).

The dominant dependence of Ṁ on just L_* and Z_* within our O-star grid may seem a bit surprising, especially in view of CAK theory which predicts an additional dependence (see Eq. (13)). However, including an additional mass-dependence in our power-law fits to the grid does not significantly improve the fit-quality. Moreover, a multi-regression fit to for the Galactic sample shows that if the mass-luminosity 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 mass-dependence). Figure 9 shows this mass-luminosity 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 O-star grid.

On the other hand, when running some additional test-models outside the range of our O-star grid, by assuming a fixed luminosity but changing the mass, we do find that including an additional mass-dependence sometimes can improve the fits. Specifically, it is clear that reducing the mass for such constant-luminosity models generally tends to increase Ṁ. This is in qualitative agreement with the results expected from CAK theory (see above), however, we note that these test-models have stellar parameters that no longer fall within the O-star 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 mass-dependency 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 zoom-in 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 wind-driving elements (Fe, C, N, O, and Ar). As we note above, the mass-loss 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 mass-loss 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 driving-elements change their ionisation stages, which seem to produce a highly non-linear effect over a restricted range in the current model-grid. Although we have not attempted to explicitly account for such non-linear temperature/ionisation-effects in the fitting-relations 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 mass-loss 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 mass-loss 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 wind-momenta and mass-loss 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 wind-momentum 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 power-laws with similar exponents (x = 2.06 vs. x = 1.83), there is a very clear offset of about 0.5 dex between the two model-relations 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 mass-loss 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 mass-loss rates we find for low-luminosity stars, in particular in the SMC (reflected in the final fit-relations, 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 mass-loss 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 pre-specified β 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 mass-loss 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 mass-loss 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³, the overall systematic reductions of Ṁ found both here and by them may point towards the need for a re-consideration of the mass-loss rates normally applied in simulations of the evolution of massive O-stars (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 O-stars with different T_eff and log g. The MRL method assumes a plane parallel atmosphere and therefore does not yield a spherical mass-loss 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 mass-loss rate from . The relation shown in the bottom panel of Fig. 5 is then computed by performing a linear fit through the resulting mass-loss 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 Monte-Carlo calculations determining g_rad and the mass flux.

Fig. 11 Direct comparison of the mass-loss rates calculated in this work versus those predicted by the Vink et al. recipe. The dashed line denotes the one-to-one correspondence. The different markers show the luminosity classes consistent with previous plots.

4.3 Comparison to observations

Observational studies aim to obtain empirical mass-loss 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 high-energy X-rays. 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 mass-loss rate. Indeed, if neglected in the analysis, such wind clumping may lead to empirical mass-loss 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 high-energy X-rays (Cohen et al. 2014), these have been shown to be relatively insensitive to the effects of clumping and porosity for Galactic O-stars (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 O-star winds. The selected studies are based on X-ray diagnostics (Cohen et al. 2014), UV+optical analyses accounting for the effects of velocity-space 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 velocity-space porosity⁴). The figure shows the modified wind-momentum 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 wind-momenta, excluding the two stars with L_* ∕L_⊙ < 10⁵ 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⁵, 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 low-luminosity 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 mass-loss 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 grid-extension might eventually start to display a significant curvature in the WLR also for Galactic objects. Nonetheless, the mismatch in the onset-luminosity 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 mass-loss rates at these low luminosities, as none of the UV-based analyses cited above account for velocity-porosity.

Fig. 12 Observed wind-momenta 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⊙ < 105 (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 O-stars, 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 (dash-dotted 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 low-luminosity Galactic O-dwarfs (“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 line-statistics, 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 compilation-study of O-stars 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 mass-loss 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 fit-errors 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 mass-loss 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 low-luminosity objects in the Galaxy. This issue was also addressed in Paper I, where it was argued that the high v_∞ for such low-density 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 shock-heated, this would lower g_rad in these layers and so also potentially reduce v_∞. Moreover, if the wind is shock-heated, 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 radiation-hydrodynamical simulations of line-driven instability (LDI) induced shocks in low-density 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 steady-state simulations by running a few additional models where g_rad in layers with speeds above 1000 km s⁻¹ simply is reduced by an ad-hoc 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 steady-state simulation set-up used here) a reduction of the outer wind driving does not necessarily lead to an increased mass-loss 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 shock-heating in the outer winds in low-luminosity 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⁻¹ and decreased it to 7.5 km s⁻¹. 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 mean-slope thus is significant, there are no clear trends within the grid. As such, to obtain a first-order approximation accounting for a v_turb that deviates from the standard value 10 km s⁻¹, 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 velocity-law that connects to a quasi-hydrostatic photosphere. This is the case also for the “standard” version of FASTWIND, used here as a starting condition for our self-consistent simulations (see Sect. 2). For the prototypical simulation presented in Sect. 2.2, Fig. 4 illustrates a fit to the self-consistent 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 modification-term in dependence on the photospheric scale-height. The full expression used to match the velocity structure is (22)

where (23)

In this relation, the parameters v_∞, β₁, β₂, 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 beta-factor β₁ = 0.8. We note though, that in particular the more low-density 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 best-fit β₁ -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 velocity-parametrisation 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 mass-loss 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, hydrogen-rich stars (excluding then the classical Wolf-Rayet (WR) stars which require different prescriptions for Ṁ). As clearly illustrated by Fig. 11, the results presented here suggest that the O-star mass-loss 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 main-sequence), the lower rates will not only affect the stars there, but also impact the properties of post-main-sequence stages. One consequence is that the luminosity at which the stars end their main-sequence 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 mass-loss rates (factor 2 to 3) and find the surface rotation speeds at the end of the main-sequence 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 black-hole 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 mass-loss 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 high-metallicity 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 mass-loss 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 WR-star through wind-stripping. In low metallicity environments such as the SMC, this is generally considered to be difficult because of the strong metallicity dependence resulting in low mass-loss rates. Considering our results here of lower O-star mass-loss rates, such (steady) wind-stripping would be more difficult to obtain also in the Galaxy, potentially leading to an increase of the lower limit for the initial mass of WR-stars created by this channel. In order to explain the observed number of WR-stars, 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 Roche-lobe 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 WR-stars 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 steady-state wind models of O-stars 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 mass-loss rate and wind-momentum 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 mass-loss 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 wind-momentum rate and mass-loss 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 mass-loss 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 O-star wind-momenta and mass-loss 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 mass-loss 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 mass-loss 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 spin-down time and altering the initial mass needed in order to produce a wind-stripped Wolf-Rayet star. As such, a key follow-up work to this study will be to now extend the grid presented here to include massive stars outside the O-star domain, to incorporate our new models into simulations of massive-star 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 ESP2017-86582-C4-1-R and PID2019-105552RB-C41 (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 (CSIC-INTA) project No. MDM-2017-0737. We would also like to thank the referee for useful comments that led to additional improvements of the paper.

Appendix A Model parameters

References

All Tables

All Figures

Fig. 1 Top panel: value of Γ versus scaled radius-coordinate for 7 (non-consecutive) hydrodynamic iterations over the complete run. The starting structure (yellow) relaxes to the final converged structure (dark blue). Bottom panel: colour map of log (ferr) 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 dash-dotted 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 mass-loss 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 radius-coordinate 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 wind-momentum 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: mass-loss rates versusluminosity for all Galactic models with a linear fit as a solid black line. The dashed line is a fit through the mass-loss rates computed using the Vink et al. recipe, the dash-dotted 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 wind-momentum 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 mass-loss rate.
In the text

	Fig. 8 Value of the exponent n showing the metallicity dependence of Dmom together with a linear fit, showing the linear dependence of this exponent with log (L* ∕106 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 mass-loss 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 mass-loss rates calculated in this work versus those predicted by the Vink et al. recipe. The dashed line denotes the one-to-one correspondence. The different markers show the luminosity classes consistent with previous plots.
In the text

	Fig. 12 Observed wind-momenta 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⊙ < 105 (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