Issue 
A&A
Volume 632, December 2019



Article Number  A126  
Number of page(s)  13  
Section  Stellar atmospheres  
DOI  https://doi.org/10.1051/00046361/201936580  
Published online  16 December 2019 
New predictions for radiationdriven, steadystate massloss and windmomentum from hot, massive stars
I. Method and first results
^{1}
KU Leuven, Instituut voor Sterrenkunde,
Celestijnenlaan 200D,
3001
Leuven,
Belgium
email: jon.sundqvist@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:
27
August
2019
Accepted:
11
October
2019
Context. Radiationdriven mass loss plays a key role in the life cycles of massive stars. However, basic predictions of such mass loss still suffer from significant quantitative uncertainties.
Aims. We develop new radiationdriven, steadystate wind models for massive stars with hot surfaces, suitable for quantitative predictions of global parameters like massloss and windmomentum rates.
Methods. The simulations presented here are based on a selfconsistent, iterative grid solution to the spherically symmetric, steadystate equation of motion, using full nonlocal thermodynamic equilibrium radiative transfer solutions in the comoving frame to derive the radiative acceleration. We do not rely on any distribution functions or parametrization for computation of the line force responsible for the wind driving. The models start deep in the subsonic and optically thick atmosphere and extend up to a large radius at which the terminal wind speed has been reached.
Results. In this first paper, we present models representing two prototypical Ostars in the Galaxy, one with a higher stellar mass M_{*}∕M_{⊙} = 59 and luminosity log_{10}L_{*}∕L_{⊙} = 5.87 (spectroscopically an early O supergiant) and one with a lower M_{*}∕M_{⊙} = 27 and log_{10}L_{*}∕L_{⊙} = 5.1 (a late O dwarf). For these simulations, basic predictions for global massloss rates, velocity laws, and wind momentum are given, and the influence from additional parameters like wind clumping and microturbulent speeds is discussed. A key result is that although our massloss rates agree rather well with alternative models using comoving frame radiative transfer, they are significantly lower than those predicted by the massloss recipes normally included in models of massivestar evolution.
Conclusions. Our results support previous suggestions that Galactic Ostar massloss rates may be overestimated in presentday stellar evolution models, and that new rates might therefore be needed. Indeed, future papers in this series will incorporate our new models into such simulations of stellar evolution, extending the very first simulations presented here toward larger grids covering a range of metallicities, B supergiants across the bistability jump, and possibly also WolfRayet stars.
Key words: radiation: dynamics / radiative transfer / hydrodynamics / stars: massive / stars: massloss / stars: winds, outflows
© ESO 2019
1 Introduction
For hot, massive stars of spectral type OB, scattering and absorption in spectral lines transfer momentum from the intense radiation field of the star to the plasma, thereby providing the force necessary to overcome gravity and drive a stellar wind outflow (Lucy & Solomon 1970; Castor et al. 1975; review by Puls et al. 2008). The linedriven winds of OB stars are very strong and fast, exhibiting massloss rates Ṁ ~ 10^{−5⋯−9} M_{⊙} yr^{−1} and terminal wind speeds v_{∞} that can reach several thousand kilometers per second. Furthermore, since the radiation acceleration g_{rad} driving these outflows is dominated by lines from metal ions like Fe, C, N, O, and so on, the winds are predicted and are shown to depend significantly on stellar metallicity Z_{⋆} (e.g., Kudritzki et al. 1987; Vink et al. 2001; Mokiem et al. 2007).
The first quantitative description of such linedriven winds was given in the seminal paper by Castor et al. (1975; referred to hereafter as CAK). Using the socalled Sobolev approximation^{1} in combination with an assumed powerlaw distribution of spectral linestrengths, CAK estimated the total line force, derived a parameterization for g_{rad}, and developed a windtheory that provided very useful scaling relations of Ṁ and v_{∞} with fundamental parameters like stellar luminosity L_{⋆} and mass M_{⋆}.
Extensions of this elegant theory (Pauldrach et al. 1986; Friend & Abbott 1986) have had considerable success in explaining many basic features of stellar winds from OB stars. However, over the years evidence has accumulated that such models may not be able to make quantitative predictions with the required accuracy. Namely, although CAKtype and Sobolevbased MonteCarlo linetransfer models typically agree rather well on predicted Ṁ (Pauldrach et al. 2001; Vink et al. 2000, 2001), more recent simulations based on line radiation transfer performed in the comoving frame (CMF), or by means of a nonSobolev MonteCarlo line force, seem to suggest overall lower rates (Lucy 2007, 2010; Krtička & Kubát 2010, 2017; Sander et al. 2017). Moreover, although the studies above, as well as this paper, focus solely on the global, steadystate wind, some timedependent models with g_{rad} computed inthe frame of the observer also indicate similar mass loss reductions as compared to CAKbased simulations (Owocki & Puls 1999).
Concerning attempts to derive Ṁ directly from observations of stellar spectra, a multitude of results are scattered throughout the literature. One key uncertainty is in regards to the effects of a clumped wind, which if neglected in the analysis may lead to rates that differ by large factors for the same star, depending on which spectral diagnostic isused to estimate Ṁ (Fullerton et al. 2006). When accounting properly for such “wind clumping”, including the lightleakage effects of porosity in physical and velocity space, a few first multiwavelength studies of Galactic Ostars (Sundqvist et al. 2011; Šurlan et al. 2013; Shenar et al. 2015) suggest massloss rates that are lower than or approximately equal to those predicted by the Vink et al. (2000) standard massloss recipe. These results also agree well with Xray (Cohen et al. 2014)and infrared (Najarro et al. 2011) studies, as well as with the analysis by Puls et al. (2006) that derived upper limits on Ṁ by considering diagnostics ranging from the optical to radio. On the other hand, some extragalactic studies (RamírezAgudelo et al. 2017; Massa et al. 2017) find upper Ṁ limits that typically are higher than the Vink et al. (2000, 2001) recipe. Clearly, more effort will be needed here in order to place better observational constraints on the massloss rates from hot, massive stars.
Reduced massloss rates could further have quite dramatic consequences for model predictions of stellar evolution and feedback (e.g., Smith 2014; Keszthelyi et al. 2017). Indeed, the presence of mass loss has a deciding impact on the lives and deaths of massive stars, affecting their luminosities, chemical surface abundances, rotational velocities, and nuclear burning lifetimes, as well as ultimately determining which type of supernova the star explodes as and which type of remnant it leaves behind (Meynet et al. 1994; Langer 2012; Smith 2014). As such, a key aim of this series of papers is – in addition to providing new quantitative predictions of mass loss – to incorporate our new models directly into simulations of massivestar evolution.
In this first paper, we outline the basic methods of our new wind simulations, which were computed using steadystate hydrodynamics and nonlocal thermodynamic equilibrium (NLTE)^{2} CMF radiative transfer calculations within the computercode network FASTWIND (SantolayaRey et al. 1997; Puls et al. 2005; Carneiro et al. 2016; Puls 2017; Sundqvist & Puls 2018). Moreover, we present and discuss some first results from two selected models of Ostars in the Galaxy. The second paper will present results from a larger grid of such Ostar simulations, covering metallicities of the Galaxy and the Magellanic Clouds (Björklund et al., in prep.). The third paper will then extend these calculations to Bstars, examining in particular massloss and windmomentum properties over the socalled bistability jump (where iron recombines from being three to two times ionized, thus providing many more driving lines, Pauldrach & Puls 1990; Vink et al. 2000). In a fourth paper we will fully incorporate our predictions into models of stellar evolution, and carefully investigate the corresponding effects (Björklund et al., in prep.). Further followup studies will then focus on, for example, timedependent effects and windclumping (see also Sect. 4.3) and potentially even an extension of our current models to the highly evolved, classical WolfRayet stars whose supersonic winds are optically thick also for optical continuum light.
2 Basic equations and assumptions
Global models (containing both the subsonic stellar photosphere and the supersonic wind outflow) are computed by a numerical, iterative gridsolution to the steadystate (timeindependent) radiationhydrodynamical conservation equations of mass and momentum in spherical symmetry. These are supplemented by also considering the energy balance, either by a method based on fluxconservation and the thermal balance of electrons (Puls et al. 2005) or by a simplified temperature structure derived from fluxweighted mean opacities in a spherical, diluted envelope (Lucy 1971).
For gravitational acceleration g(r) and (isothermal) sound speed (1)
where T(r) is the temperature, μ(r) the mean molecular weight, k_{b} is Bolzmann’s constant, and m_{h} is the hydrogen atom mass, the equation of motion (e.o.m.) for radial velocity v(r) is (2)
The radiative acceleration (3)
where stellar luminosity L_{*} is a fundamental input parameter of the model, c is the speed of light, and the fluxweighted mean opacity κ_{F} (cm^{2} g^{−1}): (4)
Here, opacities κ_{ν} and radiative fluxes F_{ν} are evaluated in the frame comoving with the fluid (the comoving frame, CMF); to order v∕c, this CMFevaluated g_{rad} can then be used directly in the inertial frame e.o.m. in Eq. (2) (e.g., Mihalas 1978, their Eqs. (15.102)–(15.113)). The gravity (5)
with gravitation constant G, is computed from the second fundamental input parameter stellar mass M_{*}. The necessary scaling radius for setting up the adaptive radial mesh is obtained as described below, giving a “stellar radius”: (6)
for the spherically modified fluxweighted optical depth (7)
The stellar “effective temperature” is then here defined as (8)
for the StefanBoltzmann constant σ. We note that this differs, for example, from the models by Sander et al. (2017), where the stellar radius and effective temperature are defined at a Rosseland (rather than fluxweighted) optical depth τ_{Ross} = 20.
Finally, massconservation provides the density structure for a steadystate massloss rate Ṁ, (9)
For a given (also prescribed) chemical composition and metallicity Z_{*}, g_{rad} (r) and a(r) are derived from the density, velocity, and temperature structure by using our radiative transfer and model atmosphere computercode package FASTWIND (SantolayaRey et al. 1997; Puls et al. 2005; Carneiro et al. 2016; Puls 2017; Sundqvist & Puls 2018); FASTWIND solves for the population numberdensity rate equations in statistical equilibrium (typically referred to as NLTE) within the spherically symmetric, extended envelope (containing both the subsonic stellar photosphere and the supersonic wind outflow), including the effects from millions of metal spectral lines on the radiation field. In the version of FASTWIND used here, chemical abundances are scaled to the solar values by Asplund et al. (2009) and a standard helium number abundance Y_{He} = n_{He}∕n_{H} = 0.1, T(r) is derived either as before (from a fluxcorrection method in the lower atmosphere and the thermal balance of electrons in the outer, Puls et al. 2005), or by the simplified approach described in the following section. Most importantly, the radiation field – in particular the g_{rad}(r) term responsible for the wind driving – is now selfconsistently computed from full solutions of the radiative transfer equations in the CMF for all included chemical elements and contributing lines (Puls 2017); the resulting g_{rad} is then obtained by direct integrations of the CMF frequencydependent fluxes and opacities.
3 Numerical calculation methods
Each model starts either from the structure obtained in a previous calculation or from a structure as computed in previous versions of FASTWIND (see Sect. 2 in SantolayaRey et al. 1997). Regarding the latter (standard) option, the key point is that there, a quasihydrostatic atmosphere below a transition velocity v ≈ 0.1a(T_{eff}) is smoothly combined with an analytic “β” wind velocity law , with b being a constant derived from the transition point velocity. This means that the equation of motion is generally not satisfied above v ≈ 0.1a and that Ṁ, v_{∞}, and β are simply provided as input by the user. By contrast, the iterative method presented here computes predictive values for Ṁ and v(r) for given fundamental parameters L_{*}, M_{*}, R_{*}, and Z_{*}.
The adaptive radial grid consists of 67 discrete points for computing the NLTE number densities and the radiative acceleration g_{rad}. Whether or not the radial grid is updated between two successive iterations depends on the current distribution of grid points in column mass, velocity, and radius. In the lowvelocity, optically thick parts of the atmosphere, the grid is adjusted according to a fixed distribution in column mass; around the sonic point and in the wind parts, the grid is adjusted according to the distribution of velocities and radial points. For the solution of the e.o.m. Eq. (2), a denser microgrid is used, which is produced by simple logarithmic interpolations from the coarser grid used in the NLTE network.
3.1 Velocity and density structure
For given values of g_{rad}(r) and a(r), the ordinary differential equation (ODE) Eq. (2) is solved by a simple RungeKutta method to obtain v(r) on the discrete radial mesh. Since the radiative acceleration here is an explicit function of only radius, Eq. (2) is critical at the sonic point v = a (see discussion in Sect. 5.4); the velocity gradient there is therefore derived by applying the rule of l’Hospital rule, and the rest of the velocity structure is then obtained by shooting up and down from the sonic point. The adaptive upper boundary typically lies between r∕R_{*} = 100 − 140 and the adaptive lower boundary is defined by computing down to a fixed column mass . Adopting such a fixed stabilizes the iteration cycle somewhat, and is achieved here by using a linear regression to a Kramerlike parameterization for the fluxweighted opacity in the deepest quasistatic layers satisfying both v≲0.1 km s^{−1} and (i.e., only in layers below R_{*} and well below the critical sonic point): (10)
where a^{2} is in cgsunits, κ_{e} is the opacity due to Thomson scattering, and k_{a} and k_{b} the corresponding fitcoefficients. This is essentially the same method as in the standard version of the FASTWIND code (SantolayaRey et al. 1997), except that we here parametrize in sound speed a instead of in temperature T. We note that although this parametrization indeed provides relatively good fits for the deep layers (see Fig. 1), it is essential to point out that it can only be applied in regions with very low velocity, where the Doppler shift effect upon the opacity is negligible^{3}.
For all layers above v ≈ 0.1 km s^{−1}, the g_{rad}(r) computed from the CMF solution to the transfer equation is used directly, and so does not explicitly depend on density, velocity, or the velocity gradient. Therefore, we must find an alternative way to update the iterative massloss rate other than for example the singularity and regularity conditions applied in CAK theory. This massloss rate update can be done in various ways, and to this end we have tried several options. For example, we preserved the total column mass down to some fixed lower boundary radius between successive iterations (Pauldrach et al. 1986; Gräfener & Hamann 2005) or, in analogy with timedependent models of linedriven Ostar winds, we tried fixing the density at some lower boundary radius r_{0} deep in the atmosphere and let the corresponding velocity v_{0} float and iteratively adjust itself to the overlying wind conditions (e.g., Owocki et al. 1988; Sundqvist & Owocki 2013). Just like Sander et al. (2017), however, we find that (at least among the methods tried thus far) the most stable way to update Ṁ is to instead consider the error in the force balance at the previous critical point r_{cr} (v = a) after a new g_{rad} has been computed for a given velocity and soundspeed structure. Defining (11)
this factor should equal Γ = g_{rad}∕g at the sonic point if the e.o.m. Eq. (2) is fulfilled. In general within the iteration cycle however, this will not be the case unless a δΓ = f_{rc} − Γ is added to the current estimate of Γ. Since it is well known (Castor et al. 1975) that for a linedriven wind g_{rad} ∝ 1∕Ṁ^{α}, with α being some power, we may thus update our massloss rate for iteration i + 1 according to . Again quite analogous to Sander et al. (2017), we find that it is typically easiest for the stability of the iterationcycle to simply assume α = 1. We stress, however, that this (of course) does not mean that really g_{rad} ∝ 1∕Ṁ, but only that this provides a good way of updating the iterative massloss rate towards convergence.
Fig. 1 Final force balance for the early Ostar model in Table 1. Lower panels of both figures: the solid lines and black squares show Γ = g_{rad} ∕g on the discrete mesh points and the red dashed line values for the rest of the terms in the e.o.m., Eq. (2), i.e., (1 − 2a^{2}∕r + da^{2}∕dr + vdv∕dr(1 − a^{2}∕v^{2}))∕g. Since the e.o.m. is perfectly fulfilled, the solid black and dashed red lines lie directly on top of each other in the figure. The dotted lines mark the sonic point v = a. Upper panels: absolute value of the error f_{err} (Eq. (13)) in the e.o.m. The abscissae in the upper figure show velocity v and in the lower scaled radius r∕r_{c} − 1, with r_{c} the lower boundary radius of the simulation. 

Open with DEXTER 
3.2 Temperature structure
The radiallydependent sound speed is obtained from a calculation of the temperature structure (and a radially dependent mean molecular weight). As mentioned above, the temperature update can be done selfconsistently by means of a fluxconservation+electron thermal balance method. In this case, the calculation is done in parallel with computing a new estimate of g_{rad}, requiring that both temperature and radiative acceleration be converged before the next update of velocity and density. However, since converging temperature and radiative acceleration in parallel is very timeconsuming, and sometimes leads to somewhat unstable iteration cycles, most models here instead update the temperature structure using a simplified method. Namely, following Lucy (1971) we may estimate the temperature structure in the spherically extended envelope as (12)
for dilution factor and where for radii r < R_{*} simply W = 1∕2. Analogous to the standard version of FASTWIND, we also here impose a floortemperature of typically T ≈0.4 T_{eff} in order to prevent excessive cooling in the outermost wind (see also Puls et al. 2005). Given current estimates for g_{rad}(r) and μ(r), Eq. (12) can now be formulated as an ODE for da^{2}∕dr and solved in parallel with the e.o.m. Eq. (2). In contrast to the flux conservation+electron thermal balance method, this temperature structure is then held fixed (along with the new v(r) and ρ(r)) within the next NLTE iteration loop for obtaining the updated g_{rad}. Although this procedure means we do not require perfect radiative equilibrium, we have found that the corresponding impact on the wind dynamics is small (see also Pauldrach et al. 1986). As such, the models in Sect. 4 all assume this structure, except for in Sect. 4.2 where a detailed comparison to a model computed from full fluxconservation+electron thermal balance is presented.
3.3 Radiative acceleration
As mentioned in the previous section, the radiative acceleration is derived from the NLTE rate equations and full solutions of the radiative transfer equation in the CMF, for all elements from H to Zn and including “all” contributing lines (i.e., all those included in the corresponding model atoms). This is different from previous versions of FASTWIND aimed at quantitive spectroscopy of optical and infrared spectral lines, where opacities from elements not included in the detailed spectroscopic investigations (the “background” elements) were added up to form the quasicontinuum quantities used to compute the radiation field. Although the CMF transfer solver is a new addition to the code (Puls 2017), the corresponding model atoms are the same as in previous versions (coming primarily from the Munich database compiled by Pauldrach et al. 2001).
Numerous tests have been performed to verify the validity of the new CMF solver for FASTWIND; for example an extensive comparison of “βlaw” models with the wellestablished CMFGEN code (Hillier & Miller 1998) shows excellent agreement between the radiative accelerations computed by the two independent programs (Puls et al., in prep., see also Fig. 1 in Puls 2017).
The radiative acceleration g_{rad} depends on opacities and fluxes, where the opacities are computed from a consistent NLTE solution that itself depends on the (CMF)radiation field. After a consistent solution for the frequencydependent opacities and fluxes have been obtained, the final g_{rad} is derived from direct numerical integrations. If the simplified temperature structure (see above) is applied, the procedure isperformed for a given velocity, density, and temperature, and the NLTE cycle iterated until g_{rad} is converged. This g_{rad} is then usedto compute a new temperature, velocity, and density, and so forth. On the other hand, if the flux conservation+electron thermal balance method is applied, the temperature is computed in parallel with g_{rad}, and the next hydrodynamic update of density and velocity is only done when both have converged. This often leads to more unstable iteration cycles and significantly slower model convergence, particularly when relatively big changes between successive hydrodynamic iterations are present.
Finally, just like in previous versions of FASTWIND (and as in basically all other 1D model atmosphere codes), we broaden all spectral line profiles with an additional isotropic “microturbulent” velocity v_{turb}, set here to a standard Ostar valueof 10 km s^{−1} (but see discussion in Sect. 4.4 regarding the impact of this parameter on the wind dynamics). We note that this v_{turb} also enters into the calculations of the NLTE number densities and opacities due to its influence on the radiative interaction rates.
3.4 Convergence criteria
The basic convergence criterion for our models is that the e.o.m. be fulfilled after a new g_{rad}(r) has been calculated from a given density and velocity structure. Rewriting the e.o.m. Eq. (2) as (13)
we note that f_{err} in Eq. (13) is zero in a dynamically consistent model. We thus require that (14)
be below a certain threshold, typically set to 0.5–1 × 10^{−2} (and neglecting the deepest quasistatic layers where the Kramerlike approximation is used for g_{rad}, see above); this means that for a converged model the e.o.m. is everywhere fulfilled to within better than 1%.
In addition, we also make sure that sufficient convergence (again typically better than 1%) is obtained for the hydrodynamic variables Ṁ, v, and T using the maximum relative change between twosuccessive iterations i − 1 and i. That is, for quantity X we require that (15)
be below ≈0.01.
A key point with the above is that the main convergence criterium Eq. (14) is applied directly on the actual e.o.m. Since the e.o.m. is by definition fulfilled for a converged model, this means that we should avoid some potential difficulties related to “false convergence”, which could (at least in principle) occur for criteria involving only relative changes between iterations of the hydrodynamic variables (such as Eq. (15)). On the other hand, there is nothing in the above that guarantees that the converged solution is also unique; this is further discussed in Sect. 5.4.
4 First results
Table 1 displays fundamental stellar parameters and the final predicted values for Ṁ and terminal wind speed v_{∞} (defined here simply as the velocity at the outermost grid point) for two models, representing a prototypical highermass and lowermass Ostar in the Galaxy; following the calibration by Martins et al. (2005) the two models would be classified as O4I and O7V stars, respectively. This section presents results from these two base models, starting with a detailed analysis of the former.
4.1 Early Ostar model
For the more luminous, highermass model, M ≈ 60 M_{⊙}, Fig. 1 displays the final, converged force balance of the simulation, demonstrating a perfect agreement between left and righthand sides in the e.o.m Eq. (2). We note further that even the lowermost part where the Kramerlike approximation (see above) is applied gives a force balance accurate to within a few percent (see upper panel) for the final fitparameters k_{a} = 10^{13.19} and k_{b} = 4.611.
The predicted massloss rate of the model is 1.51 × 10^{−6} M_{⊙} yr^{−1}, the terminal wind speed is 2480 km s^{−1}, and the scaled windmomentum rate η = Ṁv_{∞}c∕L_{*} = 0.25. Above the quasistatic deep layers the velocity field (illustrated in Fig. 2) displays a very steep acceleration around the sonic point, followed by a region of slower acceleration around ~ 100 km s^{−1}. Due to the (here implicit) dependence of g_{rad} on the velocity gradient in these highly supersonic regions, this feature also causes the corresponding radiative acceleration plateau visible in Fig. 1 around ~ 100 km s^{−1}. The radially modified fluxweighted optical depth at the sonic point is and the stellar radius (per definition at ) is located at a velocity v ≈ 0.2 km s^{−1}, which is well below the sonic point. We further also note the dip in radiative acceleration in the atmospheric regions leading up to the sonic point (see Fig. 1). This is qualitatively similar to the behavior observed in some earlier nonSobolev simulations based on simplified methods to compute g_{rad} (Owocki & Puls 1999), however further studies are required to determine potential connections to those models in a more quantitativeway. Figure 2 displays the accompanying temperature structure, confirming that but also illustrating that indeed no temperature reversal (see Sect. 4.2) is seen in this model computed by means of the simplified Eq. (12).
Figure 3 displays some characteristics of the iteration cycle towards hydrodynamic convergence. The top panel shows the maximum error in the e.o.m., (Eq. (14)), the second and third panels the maximum relative change in velocity and temperature between two iterations (Eq. (15)), and the two bottoms panels give the iterative evolution of massloss rate and terminal wind speed. The figure demonstrates how, as the error in the e.o.m. decreases, the massloss rate, velocity field, and temperature structure all stabilize toward convergence. Figure 4 displays the combination of and Ṁ for each iteration i, showing again how after an initial adjustment period the massloss rate remains fairly constant while decreases towards convergence. We note again that for each such hydrodynamic iterationupdate i, the NLTE and CMF radiative transfer computations iterate toward a new converged radiative acceleration for the given , , and . Figure 4 demonstrates this, showing Δg_{rad,j} over the NLTE iteration cycle j; in this figure, the valleys in Δg_{rad,j} correspond to the places at which g_{rad} has relaxed and a new update i of the hydrodynamic variables is performed.
For this simulation, the calculation started from a standard FASTWIND model with a “β” velocity law and consisted in total of 345 NLTE and CMF radiative transfer iterations for computing the radiative acceleration, corresponding here to 35 hydrodynamic iteration updates. The converged final model shows log_{10}ΔT_{i} = −4.3, log_{10}Δv_{i} = −2.5, log_{10}Δv_{∞,i} = −3.7, log_{10}ΔṀ_{i} = −3.1, and .
4.2 Influence of temperature structure
To test the influence of the simplified temperature structure of Eq. (12) on the wind dynamics, this section computes a model with identical input parameters as for the early Ostar in Table 1, but now using a full fluxconservation+electron thermal balance method (see above) to derive the run of the temperature. Indeed, this model can also be brought to convergence, but only after more than twice the number of NLTE iterations and a careful setup of the initial “βlaw” startmodel.Nonetheless, this early Ostar model converges to the same massloss rate Ṁ = 1.5 × 10^{−6} M_{⊙} yr^{−1} and terminal wind speed v_{∞} = 2500 km s^{−1} as previously. Figure 5 plots the temperature and velocity structures for the two converged models, demonstrating that although the temperatures as expected reveal some differences (e.g., a small reversal is now visible; see discussion in Puls et al. 2005 for more details about this feature), the resulting velocity fields are remarkably similar.
Further, we can also quite readily understand this result by: (i) noting that although g_{rad} in the optically thick, quasistatic layers is strongly dependent on the local gas temperature, in the thinner supersonic regions it becomes almost independentof it, and (ii) inspecting the force balance at the sonic point and beyond, which reveals that in those regions the “Parker terms” 2a^{2} ∕r −da^{2}∕dr are less than 1% of the corresponding values of g_{rad} (i.e., for the supersonic part the e.o.m. is completely dominated by inertia and the counteracting radiative and gravitational accelerations).
Input stellar and predicted wind parameters of models.
Fig. 2 Upper panel: velocity v (solid line) and sound speed a (dashed line) as a function of scaled radius r∕r_{c} − 1, with r_{c} the radius atthe lower boundary, for the early Ostar model in Table 1. The dotted lines mark the position at which v = a. Lower panel: temperature as a function of spherically modified fluxweighted optical depth . The dotted lines mark the position at which (and thus T = T_{eff}, see text). 

Open with DEXTER 
4.3 Influence of wind clumping and Xrays
All models presented here solve the e.o.m. for a steady outflow. However, linear perturbation analysis shows that the observer’s frame line force is subject to a very strong instability on scales smaller than the Sobolev length (Owocki & Rybicki 1984). Timedependent models (Owocki et al. 1988; Feldmeier et al. 1997; Sundqvist et al. 2018) following the nonlinear evolution of this linedeshadowing instability (LDI) show a characteristic twocomponentlike structure consisting of spatially small and dense clumps separated by large regions of very rarified material, accompanied by strong thermal shocks and a highly nonmonotonic velocity field. As shown for example in Fig. 5 of Sundqvist et al. (2018) however, the averaged density and velocity in such models still exhibit a smooth and steady behavior. As such, for computations aiming at predictions of global quantities, such as massloss rates and terminal wind speeds, it may still be reasonable to solve the hydrodynamic equations in the steady limit. Nevertheless, the overdense clumps and the highenergy radiation resulting from the wind shocks may still affect the overall wind ionization balance (e.g., Bouret et al. 2005) and g_{rad}, and so might induce feedbackeffects also upon global quantities.
Here we examine such feedback effects by computing models with identical input parameters as for the early Ostar in Table 1, but now including parametrized forms of such wind clumping and Xray emissions. Xrays are incorporated according to Carneiro et al. (2016) and clumping according to Sundqvist & Puls (2018). We note however that in this paper we restrict the analysis to clumps that are optically thin; potential dynamical effects of porosity in physical and velocity space (Muijres et al. 2011; Sundqvist et al. 2014) will be examined in a forthcoming paper. This assumption makes the basic treatment of wind clumping here very similar to that in the alternative steady models by Sander et al. (2017); except for the radial distribution of clumping factors. For the models displayed in Fig. 6, we assumed that clumping starts at a velocity v = 0.1v_{∞} and increases linearly in velocity until v = 0.2v_{∞}, outside which it remains constant at a value given by a clumping factor , for clump density ρ_{cl} and mean wind density ⟨ρ⟩ = Ṁ∕(4πvr^{2}). While both observations (Puls et al. 2006; Najarro et al. 2011) and theory (Sundqvist & Owocki 2013) indicate that clumping may be a more complicated function of radius, such a simple clumpinglaw should suffice for our testpurposes here. We further note that within the approximations used in this paper, f_{cl} is the inverse to the fractional wind volume f_{vol} occupied by clumps, for example, here (but see discussion in Sundqvist & Puls 2018). Xray emissions are assumed to start at r ≈ 1.5 R_{*} and follow the basic standard prescription outlined in Carneiro et al. (2016); the model presented in Fig. 6 has a resulting Xray luminosity L_{x}∕L_{*}≈ 0.8 × 10^{−7}.
The lower panel of Fig. 6 compares Γfactors from the simulation without clumping to (i) one including only clumping but no Xrays and (ii) one including both clumping and Xrays. As shown in the figure, for this particular model both these effects lead to a boost of the radiative acceleration in the outer parts of the wind. This in turn leads to higher wind speeds as illustrated by the upper panel of Fig. 6. Due to the steep initial wind acceleration, the assumed startvelocity for clumping v = 0.1v_{∞} corresponds to a rather low radius r∕R_{*}≈ 1.05, implying clumping also in nearstar wind regions. However, the massloss rates are only marginally different, since g_{rad} around the critical point is barely affected by any clumping or Xray feedback. The model including only clumping has a final Ṁ = 1.41 × 10^{−6} M_{⊙} yr^{−1} and v_{∞} = 3150 km s^{−1} and the model including both clumping and Xray feedback results in Ṁ = 1.41 × 10^{−6} M_{⊙} yr^{−1} and v_{∞} = 3390 km s^{−1}.
We finally note that future work should also examine potential feedback effects from models that introduce clumping in even deeper and slower atmospheric layers, which may then also affect the predicted massloss rates.
Fig. 3 Iterative evolution of quantities vs. hydrodynamic iteration number i. Uppermost panel: maximum error in the e.o.m.; second panel: relative change in velocity; third: relative change in temperature; fourth: massloss rate; and fifth: terminal wind speed. See text. 

Open with DEXTER 
Fig. 4 Upperpanel: maximum error in the e.o.m. (triangles) and the corresponding massloss rate for each hydrodynamic iteration i. Lower panel: relative change in radiative acceleration between two successive NLTE iterations j − 1 and j. The vertical dotted lines mark NLTE iterations j where hydrodynamic iteration updates i are made. Both panels show results from the early Ostar model in Table 1. 

Open with DEXTER 
Fig. 5 Comparison of converged temperature (upper panel) and velocity (lower panel) for early Ostar models computed with a simplified temperature structure (solid lines) and a fluxconservation+electron thermal balance method (dashed lines). On the abscissae are radially modified fluxweighted optical depth. 

Open with DEXTER 
Fig. 6 Comparison of converged structures for early Ostar models with no wind clumping and no Xrays (solid lines), with wind clumping but without Xrays (dashed lines), and with both wind clumping and Xrays (dasheddotted lines); see text for details. Uppermost panel: velocity vs. scaled radius; lower two panels: Γ vs. velocity (upper) and scaled radius (lower). A clumping factor f_{cl} = 10, starting at v = 0.1v_{∞}, is assumed. See text. 

Open with DEXTER 
4.4 Influence of turbulent speed
The turbulent speed v_{turb} applied in our models (see Sect. 3.3) effectively broadens the line profiles and so also affects the computation of g_{rad}. As also discussed by Poe et al. (1990) and Lucy (2007) for example, this can have a significant impact upon the line acceleration in the critical sonic point region. More specifically, since a line profile extends over a few thermal+ turbulent velocity widths, for our standard Ostar value v_{turb} = 10 km s^{−1}, in the region leading up to the sonic point a ≈ 20 km s^{−1} there is not enough velocity space available to completely Doppler shift line profiles out of their own absorption shadows. This then typically leads to a reduction of the line acceleration in these regions, as compared to Sobolev calculations where it is assumed that the line profiles are always fully deshadowed so that the corresponding (foreaft symmetric) line resonance zones do not extend into the stellar core. In turn, the reduced g_{rad} may then also reduce the predicted massloss rate as the overlying wind adjusts to the “choked” subsonic conditions.
Figure 7 illustrates this effect, displaying results from models computed with identical parameters as the early Ostar in Table 1, but with varying turbulent velocities v_{turb}. The figure shows resulting massloss rates and terminal wind speeds, normalized to the standard model with v_{turb} = 10 km s^{−1}. The clearly visible trend in the figure confirms the above discussion; higher turbulent velocities mean lower predicted massloss rates. In addition, the figure shows that this then also means higher terminal wind speeds, essentially because the supersonic wind now needs to drive less mass off the stellar surface and so can accelerate more easily. Typical numbers for this early Ostar model are a reduction in mass loss by approximately a factor of two when increasing v_{turb} from 10 to 18 km s^{−1}, accompanied by an increase in terminal wind speed by approximately 50%.
Fig. 7 Comparison of early Ostar models computed for different values of “microturbulent” velocity v_{turb}. Upper panel:compares the computed terminal wind speeds v_{∞} and lower panel: massloss rates Ṁ. Both panels have been normalized to the results of the standard model with v_{turb} = 10 km s^{−1}. 

Open with DEXTER 
4.5 Late Ostar model
For the lowermass model (M ≈ 27M_{⊙}), the two upper panels of Fig. 8 show the converged force balance analogous to Fig. 1. At first glance, the lowermass model displays similar characteristics to the highermass model, however closer inspection shows that the dip in g_{rad} in subsonic layers is now somewhat more prominent, and that once the Dopplershifted lineopacity kicks in near the sonic point, g_{rad} very quickly shoots up tovalues higher than 20 times the local gravity (compare Fig. 1 for the highermass model where this rise rather results in Γ ≈ 3). The steep acceleration gives rise to a fast, lowdensity wind with a low final massloss rate Ṁ = 2.55 × 10^{−8} M_{⊙} yr^{−1} and high terminal wind speed v_{∞} = 5320 km s^{−1}, giving a scaled windmomentum rate η = Ṁv_{∞}c∕L_{*} = 0.05.
We can readily understand the character of the fast wind by considering the e.o.m. Eq. (2) in the supersonic, zero soundspeedlimit: (16)
Noting then from Fig. 8 that Γ ≈ Const. ≈ 27 after the initial steep rise close to the stellar surface R_{*}, Eq. (16) can be analytically integrated to (17)
yielding a terminal wind speed v_{∞} = v(r →∞) ≈ 5300 km s^{−1} for surface escape speed . This agrees very well with the numerically computed value in Table 1. Moreover, Eq. (17) suggests that the supersonic wind here should be quite well described by a β = 1∕2 velocity law; this will be confirmed in Sect. 5.2, where we compare our selfconsistent simulations to models with such analytic β fields.
Further for this lowdensity wind model, only at the sonic point and the stellar surface is located at a very low v ≈ 5.0 × 10^{−3} km s^{−1}; indeed, Fig. 8 shows an overall subsonic region characterized by significantly lower velocities as compared to the early Ostar model, reaching v ≈ 10^{−4} km s^{−1} at the lower boundary .
The iterative behavior of the hydrodynamic variables is similar to that seen in the early Ostar model displayed in Fig. 3; after an initial decrease in mass loss and increase in terminal wind speed, the simulation eventually finds its way towards convergence. We note though that numerically it is somewhat more difficult to make this model converge, for example due to the sensitivity of the e.o.m. to the very steep acceleration around the sonic point. Figure 9 shows Ṁ and the maximum error in the e.o.m. for each of the hydrodynamic iteration steps, demonstrating that although the simulation indeed reaches our convergence criteria, it does also fluctuate slightly for structures quite close to this. The final model shows log_{10}ΔT_{i} = −4.3, log_{10}Δv_{i} = −2.5, log_{10}Δv_{∞,i} = −2.5, log_{10}ΔṀ_{i} = −2.4, and .
Fig. 8 Final results for the “late” Ostar model, with parameters as in Table 1. Upper two panels: Γfactors and lower two panels: velocity and temperature structures, analogous to Figs. 1 and 2, respectively,for the “early” Ostar model. 

Open with DEXTER 
Fig. 9 Maximum error in the e.o.m. (triangles) and the corresponding massloss rate for each hydrodynamic iteration i for the “late” Ostar in Table 1. 

Open with DEXTER 
5 Discussion
Having presented basic results and features above, this section now compares the two core simulations of this study to some other Ostar wind models and results that appear in the literature.
5.1 Comparison to other steadystate models
As mentioned in previous sections, the models by Sander et al. (2017) are based on similar assumptions as the ones adopted here (CMF radiative transfer, spherically symmetric steadystate hydrodynamics, gridbased g_{rad} depending explicitly only on radius). As such, we computed an additional model adopting the same stellar parameters, T_{eff} = 42 000 K, R_{*} ∕R_{⊙} = 15.9, log_{10}L∕L_{⊙} = 5.85, v_{turb} = 15 km s^{−1}, as in Sander et al. (2017). We note however that unlike Sander et al. (2017), the simulation here used a strict solar composition for the chemical abundances and also did not include clumping; nevertheless, the models are close enough to enable a quite fair and direct comparison. Indeed, our simulation converges to Ṁ = 1.58 × 10^{−6} M_{⊙} yr^{−1}, v_{∞} = 2700 km s^{−1}, and η = Ṁv_{∞}c∕L_{*} = 0.3. For the massloss rate this is in perfect agreement with Sander et al. (2017), who also find Ṁ = 1.58 × 10^{−6} M_{⊙} yr^{−1}. On the other hand, for the terminal wind speed (and thus the wind momentum) these latter authors find slightly lower values (v_{∞} = 2000 km s^{−1} and η = 0.23) than those found here. These differences may simply be an effect of the slight differences in input parameters and basic model assumptions (abundances, clumping, definition of stellar radius; see above and Sect. 2). Overall the general agreement between the two independent calculations is encouraging, however future work should investigate in more detail the impact of for example wind clumping on the wind parameters predicted by these CMF calculations (see also Sect. 4.3).
Also, Krtička & Kubát (2010, 2017) use CMF transfer to derive the radiative acceleration in their simulations. However, these models scale the CMF line force to a corresponding Sobolevbased one, which effectively means the critical point in their e.o.m. is shifted upstream from the sonic point to the supersonic point first identified by CAK (Krtička & Kubát 2017). As such, their models may potentially have quite different convergence behavior from the simulations presented here. Nonetheless, inserting the stellar luminosities of Table 1 into the massloss scaling relation Eq. (11) in Krtička & Kubát (2017) yields log_{10}Ṁ = −5.9 and log_{10}Ṁ = −7.2 for the models with higher and lower luminosities, respectively. For the former this agrees well with the rate derived here, whereas for the latter we predict a significantly lower value (see Table 1). Recalling again that there are some important differences in modeling techniques between Krtička & Kubát (2010, 2017) and this paper, further studies are required to trace the exact origin of this discrepancy.
The most widely used massloss rates for applications like stellar evolution are those compiled by Vink et al. (2000, 2001). These are based on global energy considerations using a prescribed β velocity law and a MonteCarlo line force based on approximate NLTE computations and the Sobolev approximation. An advantage with this approach is that since the vast majority of the radiative energy is transferred to the wind in the supersonic regions (where the Sobolev approximation is valid), the method does not depend sensitively on the transonic flow properties (nor on the turbulent speed there; see Sect. 4.4 and below). Moreover, since the velocity field is parametrized, v_{∞} can readily be adjusted (e.g., to an observationally inferred value) when deriving the corresponding massloss rate. Using Z∕Z_{⊙} = 1 and the recommended in the Vink et al. (2000) massloss recipe, we find log_{10}Ṁ = −5.3 and log_{10}Ṁ = −6.6 for our early and late Ostars, respectively. These massloss rates are higher than those derived in this paper by factors of ~ 3 and ~ 9 respectively. By accounting for the higher v_{∞} of our models, the agreement can be somewhat improved upon. Moreover, we note that the solar base metallicity is somewhatlower here than in Vink et al., simply because they used solar abundances by Anders & Grevesse (1989; Z_{⊙} = 0.019) whereas this paper assumes those by Asplund et al. (2009; Z_{⊙} = 0.013). However, a simple scaling using Z = 0.013∕0.019 in the Vink et al. (2001) recipe overestimates the effect (since the baseabundance of important driving ions like iron has not significantly changed). To quantify (see also Krtička & Kubát 2007), we ran an additional model with the same parameters as the early Ostar in Table 1, but now assuming solar abundances according to Anders & Grevesse (1989). This simulation converged to a massloss rate ~20% higher than our base model (as compared to ~40% if applying the simple metallicity scaling above). As such, also with these metallicity modifications a significantdiscrepancy remains between our rates and the Vink et al. scalings (a factor ~ 2.5 in the early Ostar case). Most likely these differences arise due to a combination of the use here of CMF transfer instead of the Sobolev approximation in the critical nearstar regions (see Sect. 4.4 and also, e.g., Owocki & Puls 1999) and the different NLTE solution techniques used here and in Vink et al. As mentioned above, the global Sobolev MonteCarlo method utilized by Vink et al. is not directly very sensitive to the conditions around the critical sonic point. However, since these models are parameterized to follow a β velocity law,the method will nevertheless indirectly neglect the effects the local reduction of the nearstar CMF radiation force also have upon the supersonic velocity field and the global energy balance. Regarding the nonSobolev models by Lucy (2007, 2010), these parametrize the MonteCarlo line force as purely a function of velocity, and further solve the e.o.m. only in the planeparallel limit up to a velocity ≈4a (i.e., focusing on the wind initiation only). In any case, Lucy (2010) also reported significantly lower massloss rates than those obtained by the corresponding Vink et al. recipe, in general agreement with the above.
Finally, regarding the dependence of the massloss rate on the microturbulent speed (Sect. 4.4), we note that both Lucy (2007) and Krtička & Kubát (2010) identify the same qualitative behavior as here, namely somewhat lower rates for higher v_{turb}.
Fig. 10 Comparison of the velocity fields in the selfconsistent hydrodynamic simulations (solid lines) to fits assuming “single” (dashed lines) and “double” (dasheddotted lines) βlaws. See text. The dotted background lines compare the fits to simple β = 0.5 and β = 1.0 laws. Upper and lower panels: late and early Ostar models from Table 1, respectively. 

Open with DEXTER 
5.2 Comparison to “βlaw” models
As mentioned above, in the standard version of FASTWIND (as well as in other codes used for quantitative spectroscopy of hot, massive stars), an analytic “βlaw” is typically smoothly combined with the quasistatic photosphere in order to approximate the outflowing wind velocity field. Figure 10 compares models computed with such βlaws to the selfconsistent velocity structures above. To allow for simple visual inspections, it is useful here to recast the radial coordinate in a dimensionless parameter X = 1 − r_{c}∕r, such that β in the highvelocity parts simply becomes the slope on a loglog plot showing velocity vs. radius, that is, log_{10}v∕v_{∞} = βlog_{10}X.
The figure shows fits to the selfconsistent velocity structures using a “single” βlaw of the form: (18)
with transition radius r_{tr} and derived from a transition velocity v_{tr} (see also Sect. 3, first paragraph). For comparison, we also display fits using “double” βlaws of the form: (19)
for u_{2} = 1 − r_{tr} b_{2}∕r and . Such double βlaws are sometimes also used in spectroscopic wind studies, for example in attempts to better capture the different acceleration regions in WRstars.
Using Eq. (18) (Eq. (19)), LevenbergMarquardt bestfits for the parts with v > v_{tr} were created by varying β (and β_{2}) for a range of assumed v_{tr}. For the early Ostar model, the bestfit singlelaw model shows β = 0.68 and v_{tr}∕a(T_{eff}) = 0.54. This velocity field is compared to the selfconsistent model in the lower panel of Fig. 10 (dashed line), showing good agreement in particular for the highly supersonic parts. A transonic point v ≈ 0.5a(T_{eff}) is significantly higher than the v ≈ 0.1a(T_{eff}) typically used in FASTWIND, and needed here to prevent the β velocity field from shooting up too early. This indicates that the wind acceleration sets in at relatively high subsonic velocities in the selfconsistent models. Indeed, for the late Ostar simulation the bestfit singlelaw has β = 0.57 and an even higher v_{tr}∕a(T_{eff}) = 0.8. The upper panel of Fig. 10 compares this velocity field to the hydrodynamic one, confirming the earlier discussion in Sect. 4.5 that the near constancy of Γ in the supersonic parts here implies a steeper β ≈ 0.5.
It is further interesting to note that the double βlaws (dasheddotted lines in Fig. 10) actually do not provide significantly better fits than the simpler single laws discussed above. Namely, while for this more complex parametrization β remains wellconstrained, our fits also signal that while for the early Ostar some constraints can be placed on the β_{2} parameter, for the late Ostar β_{2} is unconstrained (no visual improvements are seen to the fits either; see the dashed and dasheddotted curves in Fig. 10). This suggests that this double βlaw is not an optimal form for characterization of the present Ostar hydrodynamic velocity structures. Rather, andas can also be visually seen in Fig. 10, the main problem with the single βlaw is with regards to capturing the strong velocity rise and curvature in nearstar regions. In this respect, some first tests suggest that a promising way may be to also account for a quasistatic exponential rise v ~ e^{Δr∕h} in the lowvelocity wind regions, with some effective scale height h, in the above parameterizations. This idea will be further explored in the followup paper, where a larger grid of hydrodynamic models will be presented.
Overall,these first comparisons suggest that for early Ostars a model with a photosphere to wind transition point v ≈ 0.5a(T_{eff}) accompanied by a β ≈ 0.7 law may be a relatively good initial choice for quantitative spectroscopic studies not aiming for theoretical predictions of the global wind variables, but rather for empirical derivations of them. On the other hand, for late Ostars with lowdensity winds a steeper β ≈ 0.5−0.6 seems to be a better choice, at least for the highly supersonic parts.
5.3 Comparison to CAK line force
In CAKtheory (Castor et al. 1975), the radiation line force in a spherically symmetric wind is: (20)
where α is the CAK powerlaw index, which physically describes the ratio of the optically thick line force contribution to the total one (Puls et al. 2000), and f_{d} is the finitedisk correction factor (Pauldrach et al. 1986; Friend & Abbott 1986): (21)
for σ = d lnv∕d lnr − 1 and . In Eq. (20) is the line strength normalization due to Gayley (1995), representing the ratio between the line force and the electron scattering force g_{e} in the case where all lines are optically thin (here, α = 0). For a given α, we may use these CAK expressions to compute the “effective” Q_{eff} that the CMF line force g_{cmf} (g_{rad} reduced by the continuum acceleration) corresponds to: (22)
Figure 11 shows this Q_{eff} as a function ofvelocity and dimensionless radius X (see above),for the converged early Ostar model in Table 1 and using a prototypical α = 2∕3 (Puls et al. 2000). Because CAK theory uses the Sobolev approximation to derive the line force, only the supersonic parts are displayed. The figure shows how Q_{eff} in the CMF model increases from about 1000 at the base to about 2000 at v ≈ 100 km s^{−1}, after which it again starts to decline toward larger radii. The decline in Q_{eff} when approaching the stellar surface here reflects the very steep acceleration around the sonic point seen in our CMF models, which implies very high velocity gradients dv∕dr and so lower values of Q_{eff} for a fixed α (see Eq. (22)).
It is interesting to note that these “effective” values indeed agree quite well both with the original recommended by Gayley (1995) and with the line statistics computations by Puls et al. (2000). In addition, the comparisons provide some support to radiationhydrodynamical simulations that for early Ostars in the Galaxy typically use a radially fixed when modeling the timedependent dynamical evolution of the wind (e.g., Sundqvist et al. 2014, 2018).
We caution, however, that a similar analysis for the lowermass, latetype Ostar would result in significantly lower Q_{eff} values (reaching maximum values of only a few hundred if again assuming α = 2∕3), due primarily to the very low massloss rate and steep initial acceleration found for this star (see previous sections and Eq. (22)). In this respect, however, it is important to point out that the CMF line force here is computed quite differently from g_{CAK} in CAK theory (not using the Sobolev approximation or any assumptions about underlying linedistribution functions). In particular, g_{cmf} here has noexplicit dependence on velocity gradient or velocity, in contrast to g_{CAK} (Eq. (20)). As discussed below, this can have important consequences not only for direct comparisons of models, but also for the very nature of the steady wind solution.
Fig. 11 Effective Q_{eff} values (Eq. (22)) for the converged early Ostar model (Table 1) and α = 2∕3, as a function of velocity (upper panel) and dimensionless radius parameter X (lower panel). 

Open with DEXTER 
5.4 Stability, uniqueness, and wind topology
The two base models presented above show a quite stable iteration behavior toward convergence. However, in our calculations we find that this is not always the case. For some models, the iteration rather first proceeds as expected, toward lower and lower and toward a more and more constant Ṁ, but then never reaches the (admittedly rather strict) convergence criteria discussed above. Instead, the model may experience an almost semicyclic behavior, wherein low values of can be observed during semiregular iteration intervals, but where these minima never reach our selected criterion for convergence. To this end, it is neither clear what exactly causes this behavior, nor if it is a numerical or physical effect (see further below).
For such simulations, for example the ones with the lowermost turbulent velocities displayed in Fig. 7, we then instead use the iteration with the lowest when selecting our final model. We note that while these simulations are not formally converged, the range in the most important targeted parameter Ṁ with acceptably low is typically quite modest. Specifically for the case mentioned above, the lowest was on the order of a few percent, making it above our convergence criterion here but below that used by Sander et al. (2017) (5%) for example; the corresponding massloss rates with % lay in the range ≈ (1.5−2.0) × 10^{−6} M_{⊙} yr^{−1}.
Indeed, since the steadystate e.o.m. Eq. (2) is nonlinear, it is not a priori certain that a unique and stable solution can be found. For timedependent wind models computed using distribution functions and the observer’s frame line force, this has been analyzed by Poe et al. (1990) and Sundqvist & Owocki (2015). For their parametrized line forces, these authors carried out topology analyses of the solution behavior at the critical sonic point^{4}. It turns out that depending for example on the behavior of the source function at the sonic point (Sundqvist & Owocki 2015), the wind topology can transition from a stable Xtype with a unique solution to an unstable nodal type exhibiting degenerate solutions. When such a wind model lies on the shallow part of this nodal topology branch, Poe et al. (1990) showed that an iteration scheme like that applied here on the timeindependent e.o.m. will not converge.
In this respect, it is important to point out that since the models here compute g_{rad} as an explicit function only of radius (thus not including any explicit velocity dependence), we are effectively forcing the solution to lie on the stable Xtype branch, irrespective of what the “true” underlying solution actually is. Indeed, such forcing toward the Xtype branch is effectively done also in the alternative models by Sander et al. (2017) and Lucy (2007), where the former authors use a line acceleration depending only on radius (like here) and the latter use one that depends only on velocity (thus neglecting the radial dependence).
In any case, a further topology analysis of numerical CMF models such as those presented here would require a completely new parametrization of the line force in order to properly examine the solution behavior for different stellar ranges. This lies beyond the scope of the present paper but should be a focus for followup work targeting detailed comparisons of our new steadystatemodels with timedependent hydrodynamic calculations.
5.5 The “weakwind” problem
The low massloss rate derived for the late Ostar in Sect. 4.5 is generally consistent with the empirical finding that observationally inferred massloss rates of late Odwarfs typically are much lower than suggested by previous theoretical models (the socalled “weakwind” problem, see, e.g., Puls et al. 2008; Marcolino et al. 2009). Indeed, for our lateO model star we obtain a rate that is approximately nine times lower than that predicted by the theoretical massloss recipe by Vink et al. (2000), as already discussed in Sect. 5.1.
On the other hand, the accompanying very high terminal wind speeds we derive have not been observed. This suggests that some additional physics may act to reduce the radiation force in the outer parts of these lowdensity winds from Galactic Odwarfs (or at least make them invisible to UV observations). In this respect, one speculation is that LDIgenerated shocks (see Sect. 4.3) might not be able to efficiently cool in such lowdensity environments, meaning that also the bulk wind is left much hotter than expected. Such a shockheated bulk wind seems to be suggested by the Xray observations of lowdensity winds by Huenemoerder et al. (2012) and Doyle et al. (2017), who also suggested that the “weakwind” problem might simply be related to insufficient cooling (see also Lucy 2012). Indeed, from soft Xray ROSAT observations already Drew et al. (1994) suggested that large portions of hot, uncooled gas might exist in the outer regions of lowdensity winds. If true, such extended hot regions could then also imply a significantly reduced radiative driving in these wind parts.
We note that this scenario is quite different from that examined in Sect. 4.3 for the early Ostar wind with a higher density. Specifically, although the included Xrays there may have a significant impact on the ionization balance of some specific elements, the temperature structure of the bulk wind is basically unaffected. The possibility of a hot bulk wind should be further examined in future work, preferably by performing LDI simulations tailored to lowdensity winds and including a proper treatment of the shockheated gas and associated cooling processes.
6 Summary and future work
We have developed new radiationdriven wind models for hot, massive stars, suitable for basic predictions of global parameters like massloss and windmomentum rates. The simulations are based on an iterative grid solution to the spherically symmetric, steadystate e.o.m., using full NLTE CMF radiative transfer solutions for calculation of the radiative acceleration responsible for the wind driving.
In this first paper, we present models representing two prototypical Ostars in the Galaxy, one “early” with a higher mass M_{*} ∕M_{⊙} ≈ 59 and one “late” with a lower mass M_{*}∕M_{⊙} ≈ 27. Table 1 summarizes the results, with predicted massloss rates Ṁ = 1.51 × 10^{−6} M_{⊙} yr^{−1} and Ṁ = 2.55 × 10^{−8} M_{⊙} yr^{−1} for the two models, respectively. In Sect. 4, we further discussed in some detail the influence on the model predictions from additional parameters like wind clumping, turbulent speeds, and choice of temperature structure calculation.
Overall, our results seem to agree rather well with other steadystate wind models that do not rely on the Sobolev approximation for computation of the radiative line acceleration (Lucy 2010; Krtička & Kubát 2017; Sander et al. 2017). A key result is that the massloss rates we predict are significantly lower than those (Vink et al. 2000, 2001) normally included in models of massivestar evolution. Two likely reasons for the lower values are the use of CMF transfer versus the Sobolev approximation and also the somewhat different nature of the NLTE solution techniques used here and by Vink et al. Our results thus add to some previous notions (e.g., Sundqvist et al. 2011; Najarro et al. 2011; Šurlan et al. 2013; Cohen et al. 2014; Krtička & Kubát 2017; Keszthelyi et al. 2017) that Ostar massloss rates may be overestimated in current massivestar evolution models and that (considering the big impact of mass loss on the lives of massive stars) new rates might be needed.
In this planned paper series, we intend to work toward providing such updated radiationdriven mass loss for direct use in stellar evolution calculations (Björklund et al., in prep.). Other followups include investigations of the massloss vs. metallicity dependence in our models, comparison to timedependent simulations, and extension of the current Ostar models to Bstars (in particular supergiants across the bistability jump) and potentially also WolfRayet stars.
Acknowledgements
We thank the referee for useful comments on the manuscript. We also thank Stan Owocki for so many fruitful discussions on linedriving. J.O.S. and R.B. also thank the KU Leuven EQUATIONgroup for their valuable inputs as well as their weeklysupply of fika, and also acknowledge the Belgian Research Foundation Flanders (FWO) Odysseus program under grant number G0H9218N. J.O.S. further acknowledges previous support from the European Union Horizon 2020 research and innovation program under the MarieSklodowskaCurie grant agreement No 656725. F.N. acknowledges financial support through Spanish grants ESP201565597C41R and ESP201786582 C41R (MINECO/FEDER).
References
 Anders, E., & Grevesse, N. 1989, Geochim. Cosmochim. Acta, 53, 197 [NASA ADS] [CrossRef] [Google Scholar]
 Asplund, M., Grevesse, N., Sauval, A. J., & Scott, P. 2009, ARA&A, 47, 481 [NASA ADS] [CrossRef] [Google Scholar]
 Bouret, J.C., Lanz, T., & Hillier, D. J. 2005, A&A, 438, 301 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Carneiro, L. P., Puls, J., Sundqvist, J. O., & Hoffmann, T. L. 2016, A&A, 590, A88 [NASA ADS] [CrossRef] [EDP Sciences] [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 [NASA ADS] [CrossRef] [Google Scholar]
 Doyle, T. F., Petit, V., Cohen, D., & Leutenegger, M. 2017, in The Lives and DeathThroes of Massive Stars, eds. J. J. Eldridge, J. C. Bray, L. A. S. McClelland, & L. Xiao, IAU Symp., 329, 395 [NASA ADS] [Google Scholar]
 Drew, J. E., Hoare, M. G., & Denby, M. 1994, MNRAS, 266, 917 [NASA ADS] [CrossRef] [Google Scholar]
 Feldmeier, A., Puls, J., & Pauldrach, A. W. A. 1997, A&A, 322, 878 [NASA ADS] [Google Scholar]
 Friend, D. B., & Abbott, D. C. 1986, ApJ, 311, 701 [NASA ADS] [CrossRef] [Google Scholar]
 Fullerton, A. W., Massa, D. L., & Prinja, R. K. 2006, ApJ, 637, 1025 [NASA ADS] [CrossRef] [Google Scholar]
 Gayley, K. G. 1995, ApJ, 454, 410 [NASA ADS] [CrossRef] [Google Scholar]
 Gräfener, G.,& Hamann, W.R. 2005, A&A, 432, 633 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Hillier, D. J., & Miller, D. L. 1998, ApJ, 496, 407 [NASA ADS] [CrossRef] [Google Scholar]
 Hubeny, I., & Mihalas, D. 2014, Theory of Stellar Atmospheres (Princeton: Princeton University Press) [Google Scholar]
 Huenemoerder, D. P., Oskinova, L. M., Ignace, R., et al. 2012, ApJ, 756, L34 [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., & Kubát, J. 2007, A&A, 464, L17 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Krtička, J., & Kubát, J. 2010, A&A, 519, A50 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Krtička, J., & Kubát, J. 2017, A&A, 606, A31 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Kudritzki, R. P., Pauldrach, A., & Puls, J. 1987, A&A, 173, 293 [NASA ADS] [Google Scholar]
 Langer, N. 2012, ARA&A, 50, 107 [NASA ADS] [CrossRef] [Google Scholar]
 Lucy, L. B. 1971, ApJ, 163, 95 [NASA ADS] [CrossRef] [Google Scholar]
 Lucy, L. B. 2007, A&A, 468, 649 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Lucy, L. B. 2010, A&A, 512, A33 [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., et al. 2005, A&A, 441, 735 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Massa, D., Fullerton, A. W., & Prinja, R. K. 2017, MNRAS, 470, 3765 [NASA ADS] [CrossRef] [Google Scholar]
 Meynet, G., Maeder, A., Schaller, G., Schaerer, D., & Charbonnel, C. 1994, A&AS, 103, 97 [Google Scholar]
 Mihalas, D. 1978, Stellar Atmospheres, 2nd edn. (San Francisco, W. H. Freeman and Co.), 650 [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 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Najarro, F., Hanson, M. M., & Puls, J. 2011, A&A, 535, A32 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Owocki, S. P., & Puls, J. 1999, ApJ, 510, 355 [NASA ADS] [CrossRef] [Google Scholar]
 Owocki, S. P., & Rybicki, G. B. 1984, ApJ, 284, 337 [NASA ADS] [CrossRef] [Google Scholar]
 Owocki, S. P., Castor, J. I., & Rybicki, G. B. 1988, ApJ, 335, 914 [NASA ADS] [CrossRef] [Google Scholar]
 Pauldrach, A. W. A., & Puls, J. 1990, A&A, 237, 409 [NASA ADS] [Google Scholar]
 Pauldrach, A., Puls, J., & Kudritzki, R. P. 1986, A&A, 164, 86 [NASA ADS] [Google Scholar]
 Pauldrach, A. W. A., Hoffmann, T. L., & Lennon, M. 2001, A&A, 375, 161 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Poe, C. H., Owocki, S. P., & Castor, J. I. 1990, ApJ, 358, 199 [NASA ADS] [CrossRef] [Google Scholar]
 Puls, J. 2017, in The Lives and DeathThroes of Massive Stars, eds. J. J. Eldridge, J. C. Bray, L. A. S. McClelland, & L. Xiao, IAU Symp., 329, 435 [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., Markova, N., Scuderi, S., et al. 2006, A&A, 454, 625 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Puls, J., Vink, J. S., & Najarro, F. 2008, A&ARv, 16, 209 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 RamírezAgudelo, O. H., Sana, H., de Koter, A., et al. 2017, A&A, 600, A81 [NASA ADS] [CrossRef] [EDP Sciences] [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]
 SantolayaRey, A. E., Puls, J., & Herrero, A. 1997, A&A, 323, 488 [NASA ADS] [Google Scholar]
 Shenar, T., Oskinova, L., Hamann, W. R., et al. 2015, ApJ, 809, 135 [NASA ADS] [CrossRef] [Google Scholar]
 Smith, N. 2014, ARA&A, 52, 487 [NASA ADS] [CrossRef] [MathSciNet] [Google Scholar]
 Sundqvist, J. O., & Owocki, S. P. 2013, MNRAS, 428, 1837 [NASA ADS] [CrossRef] [Google Scholar]
 Sundqvist, J. O., & Owocki, S. P. 2015, MNRAS, 453, 3428 [NASA ADS] [CrossRef] [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., & Owocki, S. P. 2011, A&A, 528, 64 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Sundqvist, J. O., Puls, J., & Owocki, S. P. 2014, A&A, 568, A59 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Sundqvist, J. O., Owocki, S. P., & Puls, J. 2018, A&A, 611, A17 [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]
 Vink, J. S., de Koter, A., & Lamers, H. J. G. L. M. 2000, A&A, 362, 295 [NASA ADS] [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]
NLTE here means number densities computed assuming statistical equilibrium (see book by Hubeny & Mihalas 2014 for an extensive overview).
Moreover, if we were to extend our models to even deeper layers with higher temperatures (the lower boundary in the prototypical early Ostar model presented in the following section lies at T ≈ 10^{5} K), it would also have to be abandoned in order to capture features of, e.g., the socalled ironbump at T ≈ 1.5–2 × 10^{5} K, where the shift in ionisation of ironlike elements produces a peak also in the quasistatic opacities.
All Tables
All Figures
Fig. 1 Final force balance for the early Ostar model in Table 1. Lower panels of both figures: the solid lines and black squares show Γ = g_{rad} ∕g on the discrete mesh points and the red dashed line values for the rest of the terms in the e.o.m., Eq. (2), i.e., (1 − 2a^{2}∕r + da^{2}∕dr + vdv∕dr(1 − a^{2}∕v^{2}))∕g. Since the e.o.m. is perfectly fulfilled, the solid black and dashed red lines lie directly on top of each other in the figure. The dotted lines mark the sonic point v = a. Upper panels: absolute value of the error f_{err} (Eq. (13)) in the e.o.m. The abscissae in the upper figure show velocity v and in the lower scaled radius r∕r_{c} − 1, with r_{c} the lower boundary radius of the simulation. 

Open with DEXTER  
In the text 
Fig. 2 Upper panel: velocity v (solid line) and sound speed a (dashed line) as a function of scaled radius r∕r_{c} − 1, with r_{c} the radius atthe lower boundary, for the early Ostar model in Table 1. The dotted lines mark the position at which v = a. Lower panel: temperature as a function of spherically modified fluxweighted optical depth . The dotted lines mark the position at which (and thus T = T_{eff}, see text). 

Open with DEXTER  
In the text 
Fig. 3 Iterative evolution of quantities vs. hydrodynamic iteration number i. Uppermost panel: maximum error in the e.o.m.; second panel: relative change in velocity; third: relative change in temperature; fourth: massloss rate; and fifth: terminal wind speed. See text. 

Open with DEXTER  
In the text 
Fig. 4 Upperpanel: maximum error in the e.o.m. (triangles) and the corresponding massloss rate for each hydrodynamic iteration i. Lower panel: relative change in radiative acceleration between two successive NLTE iterations j − 1 and j. The vertical dotted lines mark NLTE iterations j where hydrodynamic iteration updates i are made. Both panels show results from the early Ostar model in Table 1. 

Open with DEXTER  
In the text 
Fig. 5 Comparison of converged temperature (upper panel) and velocity (lower panel) for early Ostar models computed with a simplified temperature structure (solid lines) and a fluxconservation+electron thermal balance method (dashed lines). On the abscissae are radially modified fluxweighted optical depth. 

Open with DEXTER  
In the text 
Fig. 6 Comparison of converged structures for early Ostar models with no wind clumping and no Xrays (solid lines), with wind clumping but without Xrays (dashed lines), and with both wind clumping and Xrays (dasheddotted lines); see text for details. Uppermost panel: velocity vs. scaled radius; lower two panels: Γ vs. velocity (upper) and scaled radius (lower). A clumping factor f_{cl} = 10, starting at v = 0.1v_{∞}, is assumed. See text. 

Open with DEXTER  
In the text 
Fig. 7 Comparison of early Ostar models computed for different values of “microturbulent” velocity v_{turb}. Upper panel:compares the computed terminal wind speeds v_{∞} and lower panel: massloss rates Ṁ. Both panels have been normalized to the results of the standard model with v_{turb} = 10 km s^{−1}. 

Open with DEXTER  
In the text 
Fig. 8 Final results for the “late” Ostar model, with parameters as in Table 1. Upper two panels: Γfactors and lower two panels: velocity and temperature structures, analogous to Figs. 1 and 2, respectively,for the “early” Ostar model. 

Open with DEXTER  
In the text 
Fig. 9 Maximum error in the e.o.m. (triangles) and the corresponding massloss rate for each hydrodynamic iteration i for the “late” Ostar in Table 1. 

Open with DEXTER  
In the text 
Fig. 10 Comparison of the velocity fields in the selfconsistent hydrodynamic simulations (solid lines) to fits assuming “single” (dashed lines) and “double” (dasheddotted lines) βlaws. See text. The dotted background lines compare the fits to simple β = 0.5 and β = 1.0 laws. Upper and lower panels: late and early Ostar models from Table 1, respectively. 

Open with DEXTER  
In the text 
Fig. 11 Effective Q_{eff} values (Eq. (22)) for the converged early Ostar model (Table 1) and α = 2∕3, as a function of velocity (upper panel) and dimensionless radius parameter X (lower panel). 

Open with DEXTER  
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.