Free Access
Issue
A&A
Volume 632, December 2019
Article Number A126
Number of page(s) 13
Section Stellar atmospheres
DOI https://doi.org/10.1051/0004-6361/201936580
Published online 16 December 2019

© 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 line-driven winds of OB stars are very strong and fast, exhibiting mass-loss 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 grad 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 line-driven winds was given in the seminal paper by Castor et al. (1975; referred to hereafter as CAK). Using the so-called Sobolev approximation1 in combination with an assumed power-law distribution of spectral line-strengths, CAK estimated the total line force, derived a parameterization for grad, and developed a wind-theory 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 CAK-type and Sobolev-based Monte-Carlo line-transfer 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 co-moving frame (CMF), or by means of a non-Sobolev Monte-Carlo 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, steady-state wind, some time-dependent models with grad computed inthe frame of the observer also indicate similar mass loss reductions as compared to CAK-based 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 light-leakage effects of porosity in physical and velocity space, a few first multiwavelength studies of Galactic O-stars (Sundqvist et al. 2011; Šurlan et al. 2013; Shenar et al. 2015) suggest mass-loss rates that are lower than or approximately equal to those predicted by the Vink et al. (2000) standard mass-loss recipe. These results also agree well with X-ray (Cohen et al. 2014)and infra-red (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írez-Agudelo 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 mass-loss rates from hot, massive stars.

Reduced mass-loss 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 massive-star evolution.

In this first paper, we outline the basic methods of our new wind simulations, which were computed using steady-state hydrodynamics and non-local thermodynamic equilibrium (NLTE)2 CMF radiative transfer calculations within the computer-code network FASTWIND (Santolaya-Rey 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 O-stars in the Galaxy. The second paper will present results from a larger grid of such O-star 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 B-stars, examining in particular mass-loss and wind-momentum properties over the so-called 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 follow-up studies will then focus on, for example, time-dependent effects and wind-clumping (see also Sect. 4.3) and potentially even an extension of our current models to the highly evolved, classical Wolf-Rayet 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 grid-solution to the steady-state (time-independent) radiation-hydrodynamical conservation equations of mass and momentum in spherical symmetry. These are supplemented by also considering the energy balance, either by a method based on flux-conservation and the thermal balance of electrons (Puls et al. 2005) or by a simplified temperature structure derived from flux-weighted 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, kb is Bolzmann’s constant, and mh 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 flux-weighted mean opacity κF (cm2 g−1): (4)

Here, opacities κν and radiative fluxes Fν are evaluated in the frame co-moving with the fluid (the co-moving frame, CMF); to order vc, this CMF-evaluated grad 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 flux-weighted optical depth (7)

The stellar “effective temperature” is then here defined as (8)

for the Stefan-Boltzmann 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 flux-weighted) optical depth τRoss = 20.

Finally, mass-conservation provides the density structure for a steady-state mass-loss rate , (9)

For a given (also prescribed) chemical composition and metallicity Z*, grad (r) and a(r) are derived from the density, velocity, and temperature structure by using our radiative transfer and model atmosphere computer-code package FASTWIND (Santolaya-Rey et al. 1997; Puls et al. 2005; Carneiro et al. 2016; Puls 2017; Sundqvist & Puls 2018); FASTWIND solves for the population number-density 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 YHe = nHenH = 0.1, T(r) is derived either as before (from a flux-correction 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 grad(r) term responsible for the wind driving – is now self-consistently computed from full solutions of the radiative transfer equations in the CMF for all included chemical elements and contributing lines (Puls 2017); the resulting grad is then obtained by direct integrations of the CMF frequency-dependent 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 Santolaya-Rey et al. 1997). Regarding the latter (standard) option, the key point is that there, a quasi-hydrostatic atmosphere below a transition velocity v ≈ 0.1a(Teff) 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 grad. 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 low-velocity, 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 grad(r) and a(r), the ordinary differential equation (ODE) Eq. (2) is solved by a simple Runge-Kutta 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 rR* = 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 Kramer-like parameterization for the flux-weighted opacity in the deepest quasi-static 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 a2 is in cgs-units, κe is the opacity due to Thomson scattering, and ka and kb the corresponding fit-coefficients. This is essentially the same method as in the standard version of the FASTWIND code (Santolaya-Rey 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 negligible3.

For all layers above v ≈ 0.1 km s−1, the grad(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 mass-loss rate other than for example the singularity and regularity conditions applied in CAK theory. This mass-loss 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 time-dependent models of line-driven O-star winds, we tried fixing the density at some lower boundary radius r0 deep in the atmosphere and let the corresponding velocity v0 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 rcr (v = a) after a new grad has been computed for a given velocity and sound-speed structure. Defining (11)

this factor should equal Γ = gradg 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 δΓ = frc − Γ is added to the current estimate of Γ. Since it is well known (Castor et al. 1975) that for a line-driven wind grad ∝ 1∕α, with α being some power, we may thus update our mass-loss 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 iteration-cycle to simply assume α = 1. We stress, however, that this (of course) does not mean that really grad ∝ 1∕, but only that this provides a good way of updating the iterative mass-loss rate towards convergence.

thumbnail Fig. 1

Final force balance for the early O-star model in Table 1. Lower panels of both figures: the solid lines and black squares show Γ = gradg 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 − 2a2r + da2∕dr + vdv∕dr(1 − a2v2))∕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 ferr (Eq. (13)) in the e.o.m. The abscissae in the upper figure show velocity v and in the lower scaled radius rrc − 1, with rc 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 self-consistently by means of a flux-conservation+electron thermal balance method. In this case, the calculation is done in parallel with computing a new estimate of grad, 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 time-consuming, 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 floor-temperature of typically T ≈0.4 Teff in order to prevent excessive cooling in the outermost wind (see also Puls et al. 2005). Given current estimates for grad(r) and μ(r), Eq. (12) can now be formulated as an ODE for da2∕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 grad. 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 flux-conservation+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 quasi-continuum 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 well-established 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 grad 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 frequency-dependent opacities and fluxes have been obtained, the final grad 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 grad is converged. This grad 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 grad, 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 vturb, set here to a standard O-star 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 vturb 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 grad(r) has been calculated from a given density and velocity structure. Re-writing the e.o.m. Eq. (2) as (13)

we note that ferr in Eq. (13) is zero in a dynamically consistent model. We thus require that (14)

be below a certain threshold, typically set to 0.51 × 10−2 (and neglecting the deepest quasi-static layers where the Kramer-like approximation is used for grad, 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 higher-mass and lower-mass O-star 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 O-star model

For the more luminous, higher-mass model, M ≈ 60 M, Fig. 1 displays the final, converged force balance of the simulation, demonstrating a perfect agreement between left- and right-hand sides in the e.o.m Eq. (2). We note further that even the lowermost part where the Kramer-like approximation (see above) is applied gives a force balance accurate to within a few percent (see upper panel) for the final fit-parameters ka = 1013.19 and kb = 4.611.

The predicted mass-loss rate of the model is 1.51 × 10−6 M yr−1, the terminal wind speed is 2480 km s−1, and the scaled wind-momentum rate η = ṀvcL* = 0.25. Above the quasi-static 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 grad 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 flux-weighted 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 non-Sobolev simulations based on simplified methods to compute grad (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 mass-loss rate and terminal wind speed. The figure demonstrates how, as the error in the e.o.m. decreases, the mass-loss 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 mass-loss rate remains fairly constant while decreases towards convergence. We note again that for each such hydrodynamic iteration-update i, the NLTE and CMF radiative transfer computations iterate toward a new converged radiative acceleration for the given , , and . Figure 4 demonstrates this, showing Δgrad,j over the NLTE iteration cycle j; in this figure, the valleys in Δgrad,j correspond to the places at which grad 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 log10ΔTi = −4.3, log10Δvi = −2.5, log10Δv,i = −3.7, log10Δ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 O-star in Table 1, but now using a full flux-conservation+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” start-model.Nonetheless, this early O-star model converges to the same mass-loss 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 grad in the optically thick, quasi-static 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” 2a2r −da2∕dr are less than 1% of the corresponding values of grad (i.e., for the supersonic part the e.o.m. is completely dominated by inertia and the counteracting radiative and gravitational accelerations).

Table 1

Input stellar and predicted wind parameters of models.

thumbnail Fig. 2

Upper panel: velocity v (solid line) and sound speed a (dashed line) as a function of scaled radius rrc − 1, with rc the radius atthe lower boundary, for the early O-star model in Table 1. The dotted lines mark the position at which v = a. Lower panel: temperature as a function of spherically modified flux-weighted optical depth . The dotted lines mark the position at which (and thus T = Teff, see text).

Open with DEXTER

4.3 Influence of wind clumping and X-rays

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). Time-dependent models (Owocki et al. 1988; Feldmeier et al. 1997; Sundqvist et al. 2018) following the nonlinear evolution of this line-deshadowing instability (LDI) show a characteristic two-component-like structure consisting of spatially small and dense clumps separated by large regions of very rarified material, accompanied by strong thermal shocks and a highly non-monotonic 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 mass-loss 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 high-energy radiation resulting from the wind shocks may still affect the overall wind ionization balance (e.g., Bouret et al. 2005) and grad, and so might induce feedback-effects also upon global quantities.

Here we examine such feedback effects by computing models with identical input parameters as for the early O-star in Table 1, but now including parametrized forms of such wind clumping and X-ray emissions. X-rays 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πvr2). 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 clumping-law should suffice for our test-purposes here. We further note that within the approximations used in this paper, fcl is the inverse to the fractional wind volume fvol occupied by clumps, for example, here (but see discussion in Sundqvist & Puls 2018). X-ray 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 X-ray luminosity LxL*≈ 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 X-rays and (ii) one including both clumping and X-rays. 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 start-velocity for clumping v = 0.1v corresponds to a rather low radius rR*≈ 1.05, implying clumping also in near-star wind regions. However, the mass-loss rates are only marginally different, since grad around the critical point is barely affected by any clumping or X-ray 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 X-ray 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 mass-loss rates.

thumbnail 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: mass-loss rate; and fifth: terminal wind speed. See text.

Open with DEXTER
thumbnail Fig. 4

Upperpanel: maximum error in the e.o.m. (triangles) and the corresponding mass-loss 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 O-star model in Table 1.

Open with DEXTER
thumbnail Fig. 5

Comparison of converged temperature (upper panel) and velocity (lower panel) for early O-star models computed with a simplified temperature structure (solid lines) and a flux-conservation+electron thermal balance method (dashed lines). On the abscissae are radially modified flux-weighted optical depth.

Open with DEXTER
thumbnail Fig. 6

Comparison of converged structures for early O-star models with no wind clumping and no X-rays (solid lines), with wind clumping but without X-rays (dashed lines), and with both wind clumping and X-rays (dashed-dotted lines); see text for details. Uppermost panel: velocity vs. scaled radius; lower two panels: Γ vs. velocity (upper) and scaled radius (lower). A clumping factor fcl = 10, starting at v = 0.1v, is assumed. See text.

Open with DEXTER

4.4 Influence of turbulent speed

The turbulent speed vturb applied in our models (see Sect. 3.3) effectively broadens the line profiles and so also affects the computation of grad. 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 O-star value vturb = 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 de-shadowed so that the corresponding (fore-aft symmetric) line resonance zones do not extend into the stellar core. In turn, the reduced grad may then also reduce the predicted mass-loss 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 O-star in Table 1, but with varying turbulent velocities vturb. The figure shows resulting mass-loss rates and terminal wind speeds, normalized to the standard model with vturb = 10 km s−1. The clearly visible trend in the figure confirms the above discussion; higher turbulent velocities mean lower predicted mass-loss 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 O-star model are a reduction in mass loss by approximately a factor of two when increasing vturb from 10 to 18 km s−1, accompanied by an increase in terminal wind speed by approximately 50%.

thumbnail Fig. 7

Comparison of early O-star models computed for different values of “micro-turbulent” velocity vturb. Upper panel:compares the computed terminal wind speeds v and lower panel: mass-loss rates . Both panels have been normalized to the results of the standard model with vturb = 10 km s−1.

Open with DEXTER

4.5 Late O-star model

For the lower-mass model (M ≈ 27M), the two upper panels of Fig. 8 show the converged force balance analogous to Fig. 1. At first glance, the lower-mass model displays similar characteristics to the higher-mass model, however closer inspection shows that the dip in grad in subsonic layers is now somewhat more prominent, and that once the Doppler-shifted line-opacity kicks in near the sonic point, grad very quickly shoots up tovalues higher than 20 times the local gravity (compare Fig. 1 for the higher-mass model where this rise rather results in Γ ≈ 3). The steep acceleration gives rise to a fast, low-density wind with a low final mass-loss rate = 2.55 × 10−8 M yr−1 and high terminal wind speed v = 5320 km s−1, giving a scaled wind-momentum rate η = ṀvcL* = 0.05.

We can readily understand the character of the fast wind by considering the e.o.m. Eq. (2) in the supersonic, zero sound-speedlimit: (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 self-consistent simulations to models with such analytic β fields.

Further for this low-density 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 O-star 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 O-star 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 log10ΔTi = −4.3, log10Δvi = −2.5, log10Δv,i = −2.5, log10Δi = −2.4, and .

thumbnail Fig. 8

Final results for the “late” O-star 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” O-star model.

Open with DEXTER
thumbnail Fig. 9

Maximum error in the e.o.m. (triangles) and the corresponding mass-loss rate for each hydrodynamic iteration i for the “late” O-star 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 O-star wind models and results that appear in the literature.

5.1 Comparison to other steady-state 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 steady-state hydrodynamics, grid-based grad depending explicitly only on radius). As such, we computed an additional model adopting the same stellar parameters, Teff = 42 000 K, R*R = 15.9, log10LL = 5.85, vturb = 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 η = ṀvcL* = 0.3. For the mass-loss 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 Sobolev-based 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 mass-loss scaling relation Eq. (11) in Krtička & Kubát (2017) yields log10 = −5.9 and log10 = −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 mass-loss 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 Monte-Carlo 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 mass-loss rate. Using ZZ = 1 and the recommended in the Vink et al. (2000) mass-loss recipe, we find log10 = −5.3 and log10 = −6.6 for our early and late O-stars, respectively. These mass-loss 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 base-abundance 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 O-star in Table 1, but now assuming solar abundances according to Anders & Grevesse (1989). This simulation converged to a mass-loss 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 O-star 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 near-star 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 Monte-Carlo 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 near-star CMF radiation force also have upon the supersonic velocity field and the global energy balance. Regarding the non-Sobolev models by Lucy (2007, 2010), these parametrize the Monte-Carlo line force as purely a function of velocity, and further solve the e.o.m. only in the plane-parallel limit up to a velocity ≈4a (i.e., focusing on the wind initiation only). In any case, Lucy (2010) also reported significantly lower mass-loss rates than those obtained by the corresponding Vink et al. recipe, in general agreement with the above.

Finally, regarding the dependence of the mass-loss 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 vturb.

thumbnail Fig. 10

Comparison of the velocity fields in the self-consistent hydrodynamic simulations (solid lines) to fits assuming “single” (dashed lines) and “double” (dashed-dotted 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 O-star 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 quasi-static photosphere in order to approximate the outflowing wind velocity field. Figure 10 compares models computed with such β-laws to the self-consistent velocity structures above. To allow for simple visual inspections, it is useful here to re-cast the radial coordinate in a dimensionless parameter X = 1 − rcr, such that β in the high-velocity parts simply becomes the slope on a log-log plot showing velocity vs. radius, that is, log10vv = βlog10X.

The figure shows fits to the self-consistent velocity structures using a “single” β-law of the form: (18)

with transition radius rtr and derived from a transition velocity vtr (see also Sect. 3, first paragraph). For comparison, we also display fits using “double” β-laws of the form: (19)

for u2 = 1 − rtr b2r and . Such double β-laws are sometimes also used in spectroscopic wind studies, for example in attempts to better capture the different acceleration regions in WR-stars.

Using Eq. (18) (Eq. (19)), Levenberg-Marquardt best-fits for the parts with v > vtr were created by varying β (and β2) for a range of assumed vtr. For the early O-star model, the best-fit single-law model shows β = 0.68 and vtra(Teff) = 0.54. This velocity field is compared to the self-consistent 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(Teff) is significantly higher than the v ≈ 0.1a(Teff) 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 self-consistent models. Indeed, for the late O-star simulation the best-fit single-law has β = 0.57 and an even higher vtra(Teff) = 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 (dashed-dotted 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 well-constrained, our fits also signal that while for the early O-star some constraints can be placed on the β2 parameter, for the late O-star β2 is unconstrained (no visual improvements are seen to the fits either; see the dashed and dashed-dotted curves in Fig. 10). This suggests that this double β-law is not an optimal form for characterization of the present O-star 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 near-star regions. In this respect, some first tests suggest that a promising way may be to also account for a quasi-static exponential rise v ~ eΔrh in the low-velocity wind regions, with some effective scale height h, in the above parameterizations. This idea will be further explored in the follow-up paper, where a larger grid of hydrodynamic models will be presented.

Overall,these first comparisons suggest that for early O-stars a model with a photosphere to wind transition point v ≈ 0.5a(Teff) 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 O-stars with low-density 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 CAK-theory (Castor et al. 1975), the radiation line force in a spherically symmetric wind is: (20)

where α is the CAK power-law index, which physically describes the ratio of the optically thick line force contribution to the total one (Puls et al. 2000), and fd is the finite-disk 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 ge in the case where all lines are optically thin (here, α = 0). For a given α, we may use these CAK expressions to compute the “effective” Qeff that the CMF line force gcmf (grad reduced by the continuum acceleration) corresponds to: (22)

Figure 11 shows this Qeff as a function ofvelocity and dimensionless radius X (see above),for the converged early O-star 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 Qeff 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 Qeff 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 Qeff 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 radiation-hydrodynamical simulations that for early O-stars in the Galaxy typically use a radially fixed when modeling the time-dependent dynamical evolution of the wind (e.g., Sundqvist et al. 2014, 2018).

We caution, however, that a similar analysis for the lower-mass, late-type O-star would result in significantly lower Qeff values (reaching maximum values of only a few hundred if again assuming α = 2∕3), due primarily to the very low mass-loss 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 gCAK in CAK theory (not using the Sobolev approximation or any assumptions about underlying line-distribution functions). In particular, gcmf here has noexplicit dependence on velocity gradient or velocity, in contrast to gCAK (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.

thumbnail Fig. 11

Effective Qeff values (Eq. (22)) for the converged early O-star 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 semi-cyclic behavior, wherein low values of can be observed during semi-regular 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 mass-loss rates with % lay in the range ≈ (1.5−2.0) × 10−6 M yr−1.

Indeed, since the steady-state e.o.m. Eq. (2) is nonlinear, it is not a priori certain that a unique and stable solution can be found. For time-dependent 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 point4. 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 X-type 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 time-independent e.o.m. will not converge.

In this respect, it is important to point out that since the models here compute grad 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 X-type branch, irrespective of what the “true” underlying solution actually is. Indeed, such forcing toward the X-type 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 follow-up work targeting detailed comparisons of our new steady-statemodels with time-dependent hydrodynamic calculations.

5.5 The “weak-wind” problem

The low mass-loss rate derived for the late O-star in Sect. 4.5 is generally consistent with the empirical finding that observationally inferred mass-loss rates of late O-dwarfs typically are much lower than suggested by previous theoretical models (the so-called “weak-wind” problem, see, e.g., Puls et al. 2008; Marcolino et al. 2009). Indeed, for our late-O model star we obtain a rate that is approximately nine times lower than that predicted by the theoretical mass-loss 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 low-density winds from Galactic O-dwarfs (or at least make them invisible to UV observations). In this respect, one speculation is that LDI-generated shocks (see Sect. 4.3) might not be able to efficiently cool in such low-density environments, meaning that also the bulk wind is left much hotter than expected. Such a shock-heated bulk wind seems to be suggested by the X-ray observations of low-density winds by Huenemoerder et al. (2012) and Doyle et al. (2017), who also suggested that the “weak-wind” problem might simply be related to insufficient cooling (see also Lucy 2012). Indeed, from soft X-ray ROSAT observations already Drew et al. (1994) suggested that large portions of hot, uncooled gas might exist in the outer regions of low-density 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 O-star wind with a higher density. Specifically, although the included X-rays 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 low-density winds and including a proper treatment of the shock-heated gas and associated cooling processes.

6 Summary and future work

We have developed new radiation-driven wind models for hot, massive stars, suitable for basic predictions of global parameters like mass-loss and wind-momentum rates. The simulations are based on an iterative grid solution to the spherically symmetric, steady-state 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 O-stars 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 mass-loss 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 steady-state 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 mass-loss rates we predict are significantly lower than those (Vink et al. 2000, 2001) normally included in models of massive-star 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 O-star mass-loss rates may be overestimated in current massive-star 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 radiation-driven mass loss for direct use in stellar evolution calculations (Björklund et al., in prep.). Other follow-ups include investigations of the mass-loss vs. metallicity dependence in our models, comparison to time-dependent simulations, and extension of the current O-star models to B-stars (in particular supergiants across the bistability jump) and potentially also Wolf-Rayet stars.

Acknowledgements

We thank the referee for useful comments on the manuscript. We also thank Stan Owocki for so many fruitful discussions on line-driving. J.O.S. and R.B. also thank the KU Leuven EQUATION-group 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 Marie-Sklodowska-Curie grant agreement No 656725. F.N. acknowledges financial support through Spanish grants ESP2015-65597-C4-1-R and ESP2017-86582 C4-1-R (MINECO/FEDER).

References

  1. Anders, E., & Grevesse, N. 1989, Geochim. Cosmochim. Acta, 53, 197 [NASA ADS] [CrossRef] [Google Scholar]
  2. Asplund, M., Grevesse, N., Sauval, A. J., & Scott, P. 2009, ARA&A, 47, 481 [NASA ADS] [CrossRef] [Google Scholar]
  3. Bouret, J.-C., Lanz, T., & Hillier, D. J. 2005, A&A, 438, 301 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  4. Carneiro, L. P., Puls, J., Sundqvist, J. O., & Hoffmann, T. L. 2016, A&A, 590, A88 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  5. Castor, J. I., Abbott, D. C., & Klein, R. I. 1975, ApJ, 195, 157 [NASA ADS] [CrossRef] [Google Scholar]
  6. Cohen, D. H., Wollman, E. E., Leutenegger, M. A., et al. 2014, MNRAS, 439, 908 [NASA ADS] [CrossRef] [Google Scholar]
  7. Doyle, T. F., Petit, V., Cohen, D., & Leutenegger, M. 2017, in The Lives and Death-Throes of Massive Stars, eds. J. J. Eldridge, J. C. Bray, L. A. S. McClelland, & L. Xiao, IAU Symp., 329, 395 [NASA ADS] [Google Scholar]
  8. Drew, J. E., Hoare, M. G., & Denby, M. 1994, MNRAS, 266, 917 [NASA ADS] [CrossRef] [Google Scholar]
  9. Feldmeier, A., Puls, J., & Pauldrach, A. W. A. 1997, A&A, 322, 878 [NASA ADS] [Google Scholar]
  10. Friend, D. B., & Abbott, D. C. 1986, ApJ, 311, 701 [NASA ADS] [CrossRef] [Google Scholar]
  11. Fullerton, A. W., Massa, D. L., & Prinja, R. K. 2006, ApJ, 637, 1025 [NASA ADS] [CrossRef] [Google Scholar]
  12. Gayley, K. G. 1995, ApJ, 454, 410 [NASA ADS] [CrossRef] [Google Scholar]
  13. Gräfener, G.,& Hamann, W.-R. 2005, A&A, 432, 633 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  14. Hillier, D. J., & Miller, D. L. 1998, ApJ, 496, 407 [NASA ADS] [CrossRef] [Google Scholar]
  15. Hubeny, I., & Mihalas, D. 2014, Theory of Stellar Atmospheres (Princeton: Princeton University Press) [Google Scholar]
  16. Huenemoerder, D. P., Oskinova, L. M., Ignace, R., et al. 2012, ApJ, 756, L34 [NASA ADS] [CrossRef] [Google Scholar]
  17. Keszthelyi, Z., Puls, J., & Wade, G. A. 2017, A&A, 598, A4 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  18. Krtička, J., & Kubát, J. 2007, A&A, 464, L17 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  19. Krtička, J., & Kubát, J. 2010, A&A, 519, A50 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  20. Krtička, J., & Kubát, J. 2017, A&A, 606, A31 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  21. Kudritzki, R. P., Pauldrach, A., & Puls, J. 1987, A&A, 173, 293 [NASA ADS] [Google Scholar]
  22. Langer, N. 2012, ARA&A, 50, 107 [NASA ADS] [CrossRef] [Google Scholar]
  23. Lucy, L. B. 1971, ApJ, 163, 95 [NASA ADS] [CrossRef] [Google Scholar]
  24. Lucy, L. B. 2007, A&A, 468, 649 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  25. Lucy, L. B. 2010, A&A, 512, A33 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  26. Lucy, L. B. 2012, A&A, 544, A120 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  27. Lucy, L. B., & Solomon, P. M. 1970, ApJ, 159, 879 [NASA ADS] [CrossRef] [Google Scholar]
  28. Marcolino, W. L. F., Bouret, J., Martins, F., et al. 2009, A&A, 498, 837 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  29. Martins, F., Schaerer, D., Hillier, D. J., et al. 2005, A&A, 441, 735 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  30. Massa, D., Fullerton, A. W., & Prinja, R. K. 2017, MNRAS, 470, 3765 [NASA ADS] [CrossRef] [Google Scholar]
  31. Meynet, G., Maeder, A., Schaller, G., Schaerer, D., & Charbonnel, C. 1994, A&AS, 103, 97 [Google Scholar]
  32. Mihalas, D. 1978, Stellar Atmospheres, 2nd edn. (San Francisco, W. H. Freeman and Co.), 650 [Google Scholar]
  33. Mokiem, M. R., de Koter, A., Vink, J. S., et al. 2007, A&A, 473, 603 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  34. Muijres, L. E., de Koter, A., Vink, J. S., et al. 2011, A&A, 526, A32 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  35. Najarro, F., Hanson, M. M., & Puls, J. 2011, A&A, 535, A32 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  36. Owocki, S. P., & Puls, J. 1999, ApJ, 510, 355 [NASA ADS] [CrossRef] [Google Scholar]
  37. Owocki, S. P., & Rybicki, G. B. 1984, ApJ, 284, 337 [NASA ADS] [CrossRef] [Google Scholar]
  38. Owocki, S. P., Castor, J. I., & Rybicki, G. B. 1988, ApJ, 335, 914 [NASA ADS] [CrossRef] [Google Scholar]
  39. Pauldrach, A. W. A., & Puls, J. 1990, A&A, 237, 409 [NASA ADS] [Google Scholar]
  40. Pauldrach, A., Puls, J., & Kudritzki, R. P. 1986, A&A, 164, 86 [NASA ADS] [Google Scholar]
  41. Pauldrach, A. W. A., Hoffmann, T. L., & Lennon, M. 2001, A&A, 375, 161 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  42. Poe, C. H., Owocki, S. P., & Castor, J. I. 1990, ApJ, 358, 199 [NASA ADS] [CrossRef] [Google Scholar]
  43. Puls, J. 2017, in The Lives and Death-Throes of Massive Stars, eds. J. J. Eldridge, J. C. Bray, L. A. S. McClelland, & L. Xiao, IAU Symp., 329, 435 [NASA ADS] [Google Scholar]
  44. Puls, J., Springmann, U., & Lennon, M. 2000, A&AS, 141, 23 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  45. Puls, J., Urbaneja, M. A., Venero, R., et al. 2005, A&A, 435, 669 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  46. Puls, J., Markova, N., Scuderi, S., et al. 2006, A&A, 454, 625 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  47. Puls, J., Vink, J. S., & Najarro, F. 2008, A&ARv, 16, 209 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  48. Ramírez-Agudelo, O. H., Sana, H., de Koter, A., et al. 2017, A&A, 600, A81 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  49. 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]
  50. Santolaya-Rey, A. E., Puls, J., & Herrero, A. 1997, A&A, 323, 488 [NASA ADS] [Google Scholar]
  51. Shenar, T., Oskinova, L., Hamann, W. R., et al. 2015, ApJ, 809, 135 [NASA ADS] [CrossRef] [Google Scholar]
  52. Smith, N. 2014, ARA&A, 52, 487 [NASA ADS] [CrossRef] [MathSciNet] [Google Scholar]
  53. Sundqvist, J. O., & Owocki, S. P. 2013, MNRAS, 428, 1837 [NASA ADS] [CrossRef] [Google Scholar]
  54. Sundqvist, J. O., & Owocki, S. P. 2015, MNRAS, 453, 3428 [NASA ADS] [CrossRef] [Google Scholar]
  55. Sundqvist, J. O., & Puls, J. 2018, A&A, 619, A59 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  56. Sundqvist, J. O., Puls, J., Feldmeier, A., & Owocki, S. P. 2011, A&A, 528, 64 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  57. Sundqvist, J. O., Puls, J., & Owocki, S. P. 2014, A&A, 568, A59 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  58. Sundqvist, J. O., Owocki, S. P., & Puls, J. 2018, A&A, 611, A17 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  59. Šurlan, B., Hamann, W.-R., Aret, A., et al. 2013, A&A, 559, A130 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  60. Vink, J. S., de Koter, A., & Lamers, H. J. G. L. M. 2000, A&A, 362, 295 [NASA ADS] [Google Scholar]
  61. Vink, J. S., de Koter, A., & Lamers, H. J. G. L. M. 2001, A&A, 369, 574 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]

1

The Sobolev approximation assumes number densities and source functions are constant over a few Sobolev lengths LSob = vth∕(dvs∕ds) in direction s for ion thermal speed vth. This then allows for a local treatment of the radiative line transfer.

2

NLTE here means number densities computed assuming statistical equilibrium (see book by Hubeny & Mihalas 2014 for an extensive overview).

3

Moreover, if we were to extend our models to even deeper layers with higher temperatures (the lower boundary in the prototypical early O-star model presented in the following section lies at T ≈ 105 K), it would also have to be abandoned in order to capture features of, e.g., the so-called iron-bump at T ≈ 1.52 × 105 K, where the shift in ionisation of iron-like elements produces a peak also in the quasi-static opacities.

4

Since the velocity gradient does not appear explicitly in such observer’s frame calculations of grad, the sonic point remains the only critical point in the steady limit.

All Tables

Table 1

Input stellar and predicted wind parameters of models.

All Figures

thumbnail Fig. 1

Final force balance for the early O-star model in Table 1. Lower panels of both figures: the solid lines and black squares show Γ = gradg 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 − 2a2r + da2∕dr + vdv∕dr(1 − a2v2))∕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 ferr (Eq. (13)) in the e.o.m. The abscissae in the upper figure show velocity v and in the lower scaled radius rrc − 1, with rc the lower boundary radius of the simulation.

Open with DEXTER
In the text
thumbnail Fig. 2

Upper panel: velocity v (solid line) and sound speed a (dashed line) as a function of scaled radius rrc − 1, with rc the radius atthe lower boundary, for the early O-star model in Table 1. The dotted lines mark the position at which v = a. Lower panel: temperature as a function of spherically modified flux-weighted optical depth . The dotted lines mark the position at which (and thus T = Teff, see text).

Open with DEXTER
In the text
thumbnail 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: mass-loss rate; and fifth: terminal wind speed. See text.

Open with DEXTER
In the text
thumbnail Fig. 4

Upperpanel: maximum error in the e.o.m. (triangles) and the corresponding mass-loss 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 O-star model in Table 1.

Open with DEXTER
In the text
thumbnail Fig. 5

Comparison of converged temperature (upper panel) and velocity (lower panel) for early O-star models computed with a simplified temperature structure (solid lines) and a flux-conservation+electron thermal balance method (dashed lines). On the abscissae are radially modified flux-weighted optical depth.

Open with DEXTER
In the text
thumbnail Fig. 6

Comparison of converged structures for early O-star models with no wind clumping and no X-rays (solid lines), with wind clumping but without X-rays (dashed lines), and with both wind clumping and X-rays (dashed-dotted lines); see text for details. Uppermost panel: velocity vs. scaled radius; lower two panels: Γ vs. velocity (upper) and scaled radius (lower). A clumping factor fcl = 10, starting at v = 0.1v, is assumed. See text.

Open with DEXTER
In the text
thumbnail Fig. 7

Comparison of early O-star models computed for different values of “micro-turbulent” velocity vturb. Upper panel:compares the computed terminal wind speeds v and lower panel: mass-loss rates . Both panels have been normalized to the results of the standard model with vturb = 10 km s−1.

Open with DEXTER
In the text
thumbnail Fig. 8

Final results for the “late” O-star 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” O-star model.

Open with DEXTER
In the text
thumbnail Fig. 9

Maximum error in the e.o.m. (triangles) and the corresponding mass-loss rate for each hydrodynamic iteration i for the “late” O-star in Table 1.

Open with DEXTER
In the text
thumbnail Fig. 10

Comparison of the velocity fields in the self-consistent hydrodynamic simulations (solid lines) to fits assuming “single” (dashed lines) and “double” (dashed-dotted 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 O-star models from Table 1, respectively.

Open with DEXTER
In the text
thumbnail Fig. 11

Effective Qeff values (Eq. (22)) for the converged early O-star 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 (full-text 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 48-96 hours after online publication and is updated daily on week days.

Initial download of the metrics may take a while.