Issue 
A&A
Volume 647, March 2021



Article Number  A151  
Number of page(s)  14  
Section  Stellar atmospheres  
DOI  https://doi.org/10.1051/00046361/202039595  
Published online  25 March 2021 
Dynamically inflated wind models of classical WolfRayet stars
^{1}
Institute of Astronomy,
KU Leuven, Celestijnenlaan 200D,
3001
Leuven,
Belgium
email: luka.poniatowski@kuleuven.be
^{2}
Bartol Research Institute, Department of Physics and Astronomy, University of Delaware,
Newark,
DE
19716,
USA
^{3}
Astronomical Institute Anton Pannekoek, Amsterdam University,
Science Park 904,
1098 XH
Amsterdam,
The Netherlands
Received:
5
October
2020
Accepted:
10
December
2020
Context. Vigorous mass loss in the classical WolfRayet (WR) phase is important for the late evolution and final fate of massive stars.
Aims. We develop spherically symmetric timedependent and steadystate hydrodynamical models of the radiationdriven wind outflows and associated mass loss from classical WR stars.
Methods. The simulations are based on combining the opacities typically used in static stellar structure and evolution models with a simple parametrised form for the enhanced line opacity expected within a supersonic outflow.
Results. Our simulations reveal high massloss rates initiated in deep and hot, optically thick layers around T ≈ 200 kK. The resulting velocity structure is nonmonotonic and can be separated into three phases: (i) an initial acceleration to supersonic speeds (caused by the static opacity), (ii) stagnation and even deceleration, and (iii) an outer region of rapid reacceleration (by line opacity). The characteristic structures seen in converged steadystate simulations agree well with the outflow properties of our timedependent models.
Conclusions. By directly comparing our dynamic simulations to corresponding hydrostatic models, we explicitly demonstrate that the need to invoke extra energy transport in convectively inefficient regions of stellar structure and evolution models, in order to prevent drastic inflation of static WR envelopes, is merely an artefact of enforcing a hydrostatic outer boundary. Moreover, the dynamically inflated inner regions of our simulations provide a natural explanation for the oftenfound mismatch between predicted hydrostatic WR radii and those inferred from spectroscopy; by extrapolating a monotonic βtype velocity law from the observable supersonic regions to the invisible hydrostatic core, spectroscopic models likely overestimate the core radius by a factor of a few. Finally, we contrast our simulations with alternative recent WR wind models based on comoving frame (CMF) radiative transfer to compute the radiation force. Since CMF transfer currently cannot handle nonmonotonic velocity fields, the characteristic deceleration regions found here are avoided in such simulations by invoking an ad hoc very high degree of clumping.
Key words: stars: massloss / stars: WolfRayet / stars: winds, outflows / stars: atmospheres
© ESO 2021
1 Introduction
The evolution of stars with initial masses higher than eight times that of the Sun plays an essential role in the chemistry and dynamics of galaxies similar to our Milky Way (Crowther 2007; Doran et al. 2013; Ramachandran et al. 2018; Prantzos et al. 2018). These massive stars are a vital source of heavy elements and UV radiation, enriching their surroundings through strong radiationdriven winds (Lucy & Solomon 1970; Castor et al. 1975; Puls et al. 2008; Vink et al. 2001; Björklund et al. 2021). In the final evolutionary stages, some of those massive stars become hydrogendepleted, usually core Heburning progenitors of neutron stars or black holes (Yoon et al. 2012; Groh et al. 2013). Among these are so called classical WolfRayet (WR) stars (Wolf & Rayet 1867), which are characterised by strong spectral emission lines and a high luminosity to mass ratio, L∕M_{⋆} ~ 10^{4}L_{⊙}∕M_{⊙}, (Crowther 2007). Classical WR stars are distinct from very massive, mainsequence WR stars, which are core Hburning (de Koter et al. 1997; Crowther et al. 2010). Classical WR stars are also different in their composition and evolutionary state from hydrogendeficient WR central stars of planetary nebula (e.g. Todt et al. 2010).
As first suggested by Beals (1929), the prominent emission lines visible in WR spectra indicate strong stellar winds with high terminal speeds (v_{∞}~ 2000−3000 km s^{−1}) and massloss rates (Ṁ ~ 10^{−5}−10^{−3} M_{⊙} yr^{−1}) (Hamann et al. 2019; Sander et al. 2019). However, while the overall wind properties of massive mainsequence OBstars (see Puls et al. 2008 for a review) are rather well reproduced by the linedriven wind theory formulated first by Castor et al. (1975; CAK), this standard theory typically fails to explain the order of magnitude higher massloss rates of classical WR stars (Cassinelli 1991; Lamers & Leitherer 1993).
Nonetheless, WR winds are still thought to be radiation driven, as their high L∕M_{⋆} ratio brings them close to the limit at which the acceleration due to radiation g_{r} balances that of gravity g (i.e. close to Γ = 1, for Eddington factor Γ = g_{r}∕g). In fact, consulting the opacities used in stellar structure calculations, which are generally obtained from tabulations computed assuming a static medium (e.g. Iglesias & Rogers 1996), the Eddington limit for classical WR stars is already reached in subsurface layers with temperatures of about 150−250 kK. Since convection is highly inefficient in these layers (Gräfener et al. 2012), static models then typically display density inversions and highly inflated stellar envelopes (Ishii et al. 1999; Petrovic et al. 2006; Gräfener et al. 2012; Sanyal et al. 2015). For corresponding stellar evolution calculations, such envelope inflation is computationally difficult to treat and various numerical tricks are thus typically required to make computations tractable (Paxton et al. 2013; Ekström et al. 2012). This then leads to very hot (T_{eff} ≳100 kK) and compact (R ~ R_{⊙}) WR stellar surfaces. On the other hand, spectroscopic studies aiming to constrain the WR surface from observationally inferred T_{eff} typically deduce hydrostatic radii that are factors of ~2−3 higher than predicted by such evolution models (e.g. Crowther 2007; Sander et al. 2019; Hamann et al. 2019); this mismatch is sometimes referred to as the coreradius problem of classical WR stars.
Rather than retaining a static envelope, breaching the Eddington limit in subphotospheric layers can initiate an optically thick supersonic outflow. Nugis & Lamers (2002) and Grassitelli et al. (2018) suggested that WR massloss rates can be computed simply by considering only the conditions at the point where the sonic speed is reached (at Γ ≈ 1, since gas pressure terms typically are very small in comparison). Using the OPAL tables of Rosseland mean opacity (Iglesias & Rogers 1996), Ro & Matzner (2016) showed that supersonic velocities are indeed found in deep subphotospheric layers, but that a successful wind solution that could bring the initiated mass flux to infinity could not be found.
However, these models neglect the strong enhancement of the line opacity expected in a supersonic outflow. Previous modelling attempts including the Doppler effect in the line opacity calculations have either relied on a preassumedfixed velocity field to solve for the mass loss (e.g. see Lucy & Abbott 1993; Springmann 1994; Springmann & Puls 1998; de Koter et al. 1997), or attempted an iterative (assuming time independence) solution towards a selfconsistent velocity field and mass loss (Gräfener & Hamann 2005; Sander et al. 2020; Sander & Vink 2020).The latter models have used comovingframe (CMF) radiative transfer for the calculation of g_{r}. Such CMF transfer is a computationally intensive numerical technique that requires the velocity field to remain smooth and monotonic (see, e.g. discussion in Sander et al. 2020). Moreover, all these previous studies of classical WR outflows have been performed in the steadystate limit.
In this paper we present a first attempt to build a model that addresses both (dynamic) envelope inflation and line driving, using a hybrid opacity approach based on combining the opacities used for static stellar structure calculations with a simple variant of the standard parameterisation for line opacities in supersonic flows. This formalism then allows for computation of both steadystate and timedependent WR wind structures. The organisation of the paper is as follows: in Sect. 2, we describe our basic physical setup. We present dynamical models in the steadystate limit and compare these to corresponding static calculations in Sect. 3. Then in Sect. 4 we compare these steadystate models to full timedependent radiationhydrodynamical simulations of dynamically inflated WR outflows. In Sect. 5 we discuss our results and some open questions regarding WR stellar outflows, and finally we summarise our results and provide an outlook for future work in Sect. 6.
2 Physical model
We describe the WR wind outflow by the appropriate hydrodynamical equations of mass and momentum conservation assuming spherical symmetry: (1) (2)
Here ρ, v and P_{g} = k_{B}Tρ∕(μm_{H}) are mass density, velocity and gas pressure, where T is gas temperature, k_{B} is the Boltzmann constant, μ is the mean molecular weight, and m_{H} is the mass of the hydrogen atom, g = GM_{⋆}∕r^{2} is the gravitation acceleration for a constant stellar mass M_{⋆} and (3)
is the acceleration due to stellar radiation, for radiation flux F and fluxweighted opacity (massabsorption coefficient) κ_{F} in cm^{2} g^{−1}. Using the Eddington factor we can also write the equation of motion (e.o.m) as: (4)
The Eddington ratio can be expressed as: (5)
for a stellar luminosity L = 4πr^{2}F. In this paper, we assume that the radiative luminosity remains constant throughout the outflow. This means that we neglect a term corresponding to the work of the radiation field against gravity. We can estimate the corresponding expected luminosity variation by computing the photontiring parameter m = Ṁ∕Ṁ_{max} (Owocki et al. 2017), which is the ratio of the stellar mass loss Ṁ to the maximum amount of mass loss that the stellar luminosity can drive, Ṁ_{max} = LR_{c}∕(M_{⋆}G), where R_{c} is the core radius defined at the fixed lower boundary. Taking the stellar parameters used throughout this paper, M_{⋆} = 10 M_{⊙}, R_{c} = 1R_{⊙}, log_{10}(L∕L_{⊙}) = 5.416, and a typical order of magnitude mass loss, Ṁ ~ 5.0 × 10^{−5} M_{⊙} yr^{−1}, one gets (6)
which demonstrates that the luminosity variation is only a marginal effect for the cases considered in this paper.
In practice, especially for timedependent dynamical computations with possibly nonmonotonic velocity fields, it is not computationally feasible to derive the temperature structure, for example, from radiative equilibrium by means of full solutions to the frequencydependent radiative transfer equations. To simplify,we therefore follow the common approach of replacing the full energy equation by the Lucy (1971) analytic radiative equilibrium model for a grey, spherically symmetric, diluted atmosphere, however, replacing the grey opacity by the actual flux weighted κ_{F} (e.g. see Lucy & Abbott 1993). This allows us to write the temperature structure as: (7)
where is a core effective temperature, defined here from the lower boundary set at R_{c}, σ_{SB} is the StefanBoltzmann constant, is the dilution factor, and τ_{sp} is the spherically modified optical depth (8)
Here again the radius of our fixed lower boundary R_{c} is used as a scaling radius for the spherically modified optical depth scale. This is similar to the WR models by Sander et al. (2020) who also define their stellar radius from the hydrostatic lower boundary. It is different, however, than, for example, Sundqvist et al. (2019) who set the reference radius for the temperature structure in their Ostar models at a photospheric radius R_{ph}. An argument against using R_{ph} as a scaling radius for optically thick WR winds is, however, that R_{ph} ≫ R_{c}, such that the geometric dilution of the radiation field would be neglected for large parts of the outflow (see also discussion in Nugis & Lamers 2002). In principle, the explicit choice of the reference radius will affect the temperature structure throughout our models. However, test calculations have shown that this is a quite marginal effect for the resulting mass loss and velocities of the simulations considered in this paper. In any case, to facilitate comparison with other models we also introduce a stellar photospheric effective temperature T_{ph, eff} at this photospheric radius R_{ph}. Since the definition of a stellar photosphere in a spherically diluted stellar envelope is nontrivial^{1} and depends on the specific computation of the temperature structure, we here simply approximate (9)
such that the photospheric and core effective temperatures and radii are related through (10)
For given stellar parameters L, M_{⋆}, and R_{c}, and a known variation of κ_{F}, Eqs. (1)–(8) form the basic system of coupled differential equations under investigation in this paper. In general, κ_{F} is a complicated function that depends on chemical composition, density, temperature, velocity and redial position of the absorbing and emitting gas. In this study, we focus on two cases. First, we consider static opacities such as those used for stellar structure and evolution computations (see Ro & Matzner 2016). Due to the high continuum optical depth of WR outflows, conditions in the wind launching regions resemble those of radiation diffusion, allowing us to replace the fluxweighted opacity with a Rosseland mean, and so to utilise standard OPAL tables for the opacity mapping as a function of T and ρ (Iglesias & Rogers 1996). However, since these tabulations ignore the influence of Doppler shifts on the line opacity (which are critical for WR wind driving, e.g. Sander et al. 2020), we also consider a hybrid model accounting at least approximately for this effect. More specifically, we use a CAKlike parametrisation to add the cumulative force from an ensemble of spectral lines to the OPAL tables (see, Sect. 3.3).
3 Steadystate approximation
To build physical insight we first examine a simplified steadystate case. Equations (1) and (4) are then written as: (11) (12)
for an isothermal sound speed a^{2} = k_{B}T∕(μm_{H}). Equation (11) defines the spherically symmetric massloss rate Ṁ. Here we assume μ = 4∕3 corresponding to a fully ionised helium plasma.
It is commonly believed that the winds of WR stars are initiated in the deep subphotospheric layers where the gas temperature is T ~ 150−200 kK. In this temperature region, due to iron recombination, a significant increase of the Rosseland mean opacity is observed, sometimes called the iron opacity bump. A simple comparison of scales reveals that throughout the supersonic wind outflow the first two terms of ther.h.s. in Eq. (12) are much smaller than the radiation and gravitation terms, so that we can approximate the e.o.m with (13)
Here, the soundspeed term on the l.h.s. is retained in order to enable a mapping of the supersonic solution onto the subsonic part.
3.1 Boundary conditions
Equations (7)–(13) impose a twopoint boundary value problem, with boundaries at R_{c} and at^{2} r → ∞. The inner boundary at r = R_{c} is defined to be a sonic point, with , where a_{s} is defined to be the sound speed at this sonic point. As can be seen from Eq. (13), this implies Γ(R_{c}) = 1. Below the sonic point as the velocity rapidly diminishes from its sonic point value a_{s} inwards, the density structure is well approximated by the hydrostatic solution (e.g. see Ro & Matzner 2016; Grassitelli et al. 2018). In contrast, in the supersonic part, the density structure diverges from the structure given by the hydrostatic solution. As such, we initiate the calculation at the sonic point and only consider the supersonic part of the outflow, assuming that the subsonic region still can be described by hydrostatic equilibrium. The tabulated OPAL opacities can then be used to locate points for which Γ = 1, constraining the lower boundary density as a function of temperature. In this way, the range of possible massloss rates is set as a function of the lower boundary temperature.
The outer boundary is set by requiring that the wind temperature decreases sufficiently at large radii. Replacing density with the massloss rate using Eq. (11), the spherically modified optical depth (Eq. (8)) is used to compute the temperature structure according to Eq. (7). We compute the optical depth at the outer boundary τ_{out} by radially integrating Eq. (8) in r ∈ [R_{max}, ∞), assuming a terminal velocity v(r ≥ R_{max}) = v(R_{max}) = v_{∞} and κ_{F}(r ≥ R_{max}) = κ_{e}, where κ_{e} is the Thompson scattering mass absorption coefficient. The outer boundary temperature is then set by the maximum between applying τ_{out} in Eq. (8) and a floor value, which is typically set to a few tenths of the stellar effective temperature (see Puls et al. 2005; Sundqvist et al. 2019). By matching the inner and outer boundary constraints, we uniquely constrain the wind structure. In turn this then sets the unique massloss rate of the model from the range of possible solutions at the lower boundary.
3.2 Energy requirement
Additionally, for wind material to escape the stellar potential, the global energy requirement must be satisfied, which means that the integrated mechanical energy of the wind has to be positive. Taking the supersonic limit a^{2} ∕v^{2} ≪ 1 in Eq. (13) yields: (14)
Here v_{0} ~ a_{s} is the initial velocity and is the escape velocity from R_{c}. The first term on the r.h.s. corresponds to the work done by the stellar radiation and the second term is the work done by gravity. We can then introduce as the net energy change in units of : (15)
For the wind to escape to r →∞, this net energy change has to be positive, that is . Alternatively, we can also introduce the ratio of kinetic to core potential energy as a square of the velocity to core escape velocity ratio: (16)
which goes to as r→∞ so that by construction we require w(r →∞) > 0 in order to escape.
Summary of stellar parameters used in this paper.
3.3 Failed Winds in the static opacity limit
We first investigate the possibility of driving the stellar outflow in the static opacity limit, i.e. the limit in which we only consider OPAL opacities, thus neglecting the Doppler effect; κ_{F} ≈ κ_{OPAL} (see also Ro & Matzner 2016). To set up the problem, we choose the stellar mass and the core radius (as given in Table 1) such that they roughly represent the average value of the masstoradius ratio of the observationally inferred WR stellar population from Hamann et al. (2019) and Sander et al. (2019). The luminosity is then set such that the resulting Thompson scattering Eddington ratio Γ_{e} = 0.40, which also corresponds approximately to the observationally inferred average (Hamann et al. 2019; Sander et al. 2019). All models presented inSect. 3 and 4 of this paper use these stellar parameters, and also the OPAL tables for solar metallicity Z = 0.02 and composition by Grevesse & Noels (1993). This setup then gives a range of possible sonicpoint values of such that the condition is fulfilled (see Fig. 1a). These points are used as a single point boundary condition, after which forward integration from the sonic point is carried out. The results of these integrations for various fixed massloss rates Ṁ = 4πR_{c}ρ_{s}a_{s} are shown in Fig. 1a. The figure displays a colour map of Γ with the temperature on the abscissa and the ratio between the radiation pressure (17)
and gas pressure P_{g} on the ordinate (the ratio P_{r}∕P_{g} is a proxy for the density). The integration curves show the density and temperature structure for the different cases, starting from various positions on the Γ = 1 curve at the inner boundary R_{c} and extending outwards. Figure 1b then shows the corresponding net gain of kinetic energy (Eq. (16)), demonstrating that after an initial acceleration and increasing velocities, all of the curves start to decelerate until eventually reaching zero velocity and so terminating the outward integration (some solutions had to be formally terminated before actually reaching zero velocity, due to numerical difficulties in the steep deceleration region). This indicates that none of the potential solutions starting from the Γ = 1 curve are able to escape the stellar gravitation potential. These results are consistent with Ro & Matzner (2016), who also found that OPAL opacities were not able to sustain a radiationdriven mass loss initiated in the deep layers around the iron bump, and indicates that an additional source of opacity is required, raising the question of a missing force.
Fig. 1 (a) Colour map of Γ_{OPAL} for the stellar parameters from Table 1. Temperature is given on the abscissa and radiationtogas pressure ratio on the ordinate. is identified with a white contour. The colour bar on the right corresponds to the colour coding of the plot and gives the numerical value of Γ_{OPAL} for a given temperature and radiationtogas pressure ratio. Different line styles show the radiationtogas pressure ratio and temperature structure for the different fixed massloss rates. (b) Ratio of kinetic to core potential energy of the corresponding solutions from panel a. The abscissa here shows the radius coordinate x. 
3.4 Missing force
The natural candidate for the missing opacity source discussed above comes from the cumulative effect of spectral lines due to the Doppler effect. To estimate this we use the CAK parametrisation for a distribution of spectral lines, using here a simple approximation of taking the absolute value of the velocity gradient: (18)
where is an effective strength of the line ensemble (Gayley 1995), and α sets the CAK line distribution power index. In reality, once regions where d v ∕dr < 0 arise, it is possible for a photon proceeding from the stellar core to be Doppler shifted into resonance with the same spectral line at multiple locations in the outflow, thereby requiring a fully nonlocal computation of the line acceleration. As discussed further below, this could (in principle) be handled in a more selfconsistent manner by, for example, MonteCarlo (MC) lineforce calculations.In this initial study, however, we opt for a simpler approach, in order to make timedependent simulations feasible as well. Two possible approximations that can be taken then are to either take the absolute value of the velocity gradient or use max(0, dv∕dr). The first approach estimates (approximately) an upper limit to the CAKlike force in regions with negative velocity gradients, whereas the second approach would correspond to a lower limit. As most previous studies of linedriven dynamical flows have used (for example, in investigations of linedriven winds from rapidly rotating stars, Petrenz & Puls 2000, and around discs, Proga et al. 1998; Kee et al. 2016), we follow this standard approach in this paper. However, test computations using instead the lower limit max (0, dv∕dr) show that although some details of the radiation force profile then are changed, none of our conclusions are affected. With this approximation, using Eq. (18) for Γ_{CAK} then accounts for the force from spectral lines in supersonic regions, while Γ_{OPAL} handles the force from the continuum and spectral lines in the static limit. As such, the total radiative acceleration can be simply estimated by the sum of the CAKlike contribution and the static OPAL opacities, Γ_{tot} = Γ_{OPAL} + Γ_{CAK}.
Let us note here directly that the CAK line force is originally derived for application in winds that are optically thin for continuum radiation. In the lower layers of the dense WR outflows considered here, this will typically not be the case. A more rigorous treatment would then also need to take into account the change in continuum intensity by solving the transfer equation for line+continuum opacities. While approximations for the diffusion limit exist (Gayley et al. 1995) no general formalism has been developed. As such, in this first study of WR wind dynamics we opt for the simple radialstreaming CAK form above.
Furthermore, due to the explicit dependency of the CAK line force on d v∕dr the sonic point is now formally not a critical point (see Castor et al. 1975). However, our iteration scheme circumvents this issue by always applying the velocity gradient from the previous iteration, meaning that Γ_{tot}(r) technically is no longer an explicit function of dv∕dr. Nevertheless, since we require dv∕dr to be converged between iterations the corresponding feedback upon Γ_{tot}(r) is still (implicitly) accounted for. The agreement between the converged structures of these steadystate models and the timedependent simulations presented in the next section brings further support to this method.
A key question then becomes whether Γ_{CAK} also has a significant impact on the conditions at the optically thick lower boundary at the sonic point. Assuming typical Ostar values α =0.66, (e.g. Puls et al. 2000), approximating the velocity gradient^{3} as d v∕dr ~ Δv_{max}∕ΔR(v = v_{max}), and using the mass density and velocity computed in the static opacity limit for Ṁ = 5.0 × 10^{−5} M_{⊙} yr^{−1}, we find the Eddington ratio for the CAK force at the lower boundary to be Γ_{CAK} = g_{CAK}∕g ~ 0.1. Such a small contribution from the CAK force is not surprising due to the high density at the lower boundary^{4} (Γ_{CAK} ~ 1∕ρ^{α}, with α >0). However, from test calculations it turns out that although the CAK contribution indeed increases as we move outwards in the wind, a simple model with constant α = 0.66, is still not able to provide the force necessary to achieve a positive total wind energy and so drive the wind to infinity (see previous section).
There are, however, several indications that the line force in the outer winds of classical WR stars indeed might be significantly enhanced as compared to Otype stars. For example, the independent MC lineforce calculations by Lucy & Abbott (1993) and Springmann & Puls (1998) both find global momentum rates for WR stars that exceed those of Ostars by an order of magnitude or more. In these MC computations, the outerwind enhancement of the line force stems from the decreasing level of ionisation with increasing distance from the star in the optically thick WR outflows. Due to this shift from hightolow ionisation stages as photons move outwards in the wind, new spectral lines become available for photons to interact with. As these new lines have different interaction frequencies than those previously available in the inner wind parts, this means that the linefrequency gaps typically observed for the nearly constant (frozenin) ionisation conditions of optically thin Ostar winds are effectively closed. This then leads to photon trapping (Lucy & Abbott 1993) and efficient multiline scattering, which significantly increases the line force.
To mimic the effect of an outer wind lineforce enhancement, we introduce here a simple spatial variation α =α(r), which increases the CAK force away from the star by a modest decrease of α. To this end, we assume α(r) as a piecewise linear function with a constant α_{max} from the lower boundary to r_{1} and a constant α_{min} from r_{2} to the outer boundary, interpolating between the two values in the range [r_{1}, r_{2}]. This then allows the wind outflow initiated at the lower boundary to achieve positive total energy, and so be sustained until the upper boundary. While the specific treatment of α is done here in an ad hoc manner, the overall assumed form is quite consistent with the multiscattering lineforce enhancement found in the MC lineforce computations for a CAKlike spectral line distribution by Springmann (1994). The approximate values of α(x) that correspond to the line force computed by Springmann (1994) are given in Fig. 2 as black asterisks, with the α(x) as adapted in the standard model of this paper overplotted as a solid line.
The figure shows that the order of magnitude outerwind lineforce enhancement applied here indeed is on the same order as that found in these MC calculations.
We note, however, that the Springmann (1994) calculations were performed only for parameterised (monotonic) βtype velocity laws. As such, there is not a complete onetoone correspondence between the spatial and velocity scales in these MC calculations and those of the present paper. Nonetheless, the overall similarity with the very simple radial dependence of the α parameter adopted here is encouraging. In a followup work, we plan to couple the hydrodynamic simulations presented here directly to such MC computations, building on the method developed by Puls et al. (2000) for proper estimation of occupation numbers, and accounting fully for multiscattering effects as well as nonmonotonic velocity fields.
As outlined in detail in the following sections (for our standard values listed in Table 2) the formalism outlined above gives rise to dynamically inflated wind solutions, where the structure in the optically thick deep layers near the iron bump is primarily governed by Γ_{OPAL} but where Γ_{CAK} takes over the driving in the diluted outer regions. This then causes a typical wind structure that is characterised by a high massloss rate ignited from the lower boundary, a slow (end even partially negative) acceleration of the deeper layers, and a rapid reacceleration of the outer wind. As such, the typical velocity laws we find not only drastically deviate from the standard βtype laws assumed by most atmospheric models, but they also contain a significant region of negative acceleration (see also Fig. 14 in Pauldrach et al. 1993).
Fig. 2 Plot that demonstrates values of α (see text) required to reproduce the line force computed from MC radiative transfer calculations by Springmann (1994; asterisks). The solid line then overplots the piecewise linear function adopted in this paper to describe the corresponding variation of α. The vertical dashed lines correspond to r_{1} (x_{1}) and r_{2} (x_{2}). 
Summary of standard lineforce parameters used in this paper.
3.5 Dynamically inflated steadystate model
Based on the hybrid formalism introduced above we can now numerically solve Eqs. (7)–(13), iterating towards convergence using a simple scheme based on RungeKutta integrations. As discussed in Sect. 3.1, to constrain the massloss rate it is critical to here take into account the twopoint boundary value nature of the problem. At a basic level, every iteration is performed using two steps: (i) an insideout integration R_{c} → R_{max}, where the velocity, density, and temperature structure are computed, followed by (ii) an outsidein integration R_{max} → R_{c} updating the temperature structure using the correct boundary condition from the optical depth found in the first step. The massloss rate is then updated for the next iteration by finding the correct location on the Γ =1 curve at thelower boundary, using the updated temperature structure from the outsidein integration.
Using this scheme, Fig. 3a shows the final, converged temperature profile for our standard stellar parameters (Table 1) and choice of lineforce parameters (Table 2). The final (also converged) massloss rate for the model is Ṁ = 1.47 × 10^{−5} M_{⊙} yr^{−1}.
In addition to a high predicted massloss rate, the solution displays a nonmonotonic velocity field. The resulting velocity profile shown in Fig. 3b can be split into three basic regions: (1) an initial acceleration from the core followed by (2) stagnation and even a deceleration part (and so a negative velocity gradient, marked with crosses in the figure), and (3) an outer region of fast reacceleration to outflow velocities > 1000 km s^{−1}. The distinct behaviour of these three regions can be understood via the different contributing parts of the total radiation force.
Figure 3c illustrates this, plotting the individual contributions of OPAL and CAKlike opacities to the total radiative acceleration. From the figure, on one hand, it is clear that in the inner regions the radiation force is dominated by OPAL opacities and the CAKlike contribution is small; indeed, the sonic point conditions at the lower boundary can here be quite well approximated using only OPAL opacities, i.e. neglecting any influence from the Dopplershift enhancement of line opacities. On the other hand, in the outer regions Γ_{CAK} becomes thedominant part, causing a reacceleration of the flow and setting a high terminal wind velocity. This reacceleration ensures thatthe energy requirement of Eq. (15) is satisfied. As such, it is this part that ultimately allows the wind to be sustained all the way to the outer boundary, in contrast to the failed outflow and fallback seen in the previous section.
To explore the influence of different opacity contributions on the mass loss of the model, we explore the logarithm of the spherically modified optical depth (Eq. (8)) as a function of radius in Fig. 3d. The dashed vertical line marks the location of the stellar photosphere defined by Eq. (9). We note, however, that at the photosphere τ_{sp} (R_{ph})≠2∕3 as R_{c} is used as ascaling radius in Eq. (8). For most stars, this photospheric radius R_{ph} serves as the hydrostatic stellar radius; by contrast, for the WR stars examined in this paper, the photosphere is located well in the outflowing regions. Inspection of Fig. 3d shows that the majority of τ_{sp} is accumulated within few stellar radii in the regions dominated by OPAL opacities, thus the converged lower boundary temperature, and thus the massloss rate, is primarily controlled by Γ_{OPAL}.
Inspection of the density distribution in Fig. 3e further clarifies the terminology dynamic inflation. In this figure, we have again identified the location R_{ph} with a vertical dashed line. The figure shows a relatively slow decline of the density in the subphotospheric parts, which is somewhat reminiscent of the profiles seen in hydrostatic stellar models near the Eddington limit that undergo envelope inflation (see further below). But unlike such static inflation mod0els, the simulations here allow for a velocity field to develop; the lowdensity envelope now represents a dynamically inflated star.
3.6 Comparison to static inflation
For the same lower boundary conditions as adopted above, we can also compute static models. This can then serve as a direct comparison between the dynamic inflation found in the models above, and the static inflation that often occurs in stellar structure models neglecting the v dv∕dr term in the e.o.m. (and so assuming a hydrostatic stellar envelope throughout, e.g. Petrovic et al. 2006; Gräfener et al. 2012).
To construct such static models, we simply set v ~dv∕dr ~ 0 in our steadystate e.o.m. above. Introducing a column mass (19)
we write the resulting equation of hydrostatic equilibrium as: (20)
We further simplify the temperature calculation by assuming radiative diffusion in these optically thick hydrostatic layers: (21)
for a radiation pressure according to Eq. (17). Combining the equations for gas and radiation pressure, one obtains (22)
where now Γ = Γ_{OPAL} throughout the complete model.
This setup imposes a similar twopoint boundary value problem as previously analysed. To provide a fair comparison with those dynamical simulations, we assume the same lower boundary radius r = R_{⊙} as before and reuse the lower boundary temperature T ≈ 230 kK found from the converged hydrodynamic steadystate model presented above. However, at this inner boundary we now impose Γ < 1 as Γ ≥ 1 would lead to nonzero gas pressure at r →∞ (see also Sect. 2.2 in Gräfener et al. 2012). The outer boundary conditions for gas and radiation pressure (i.e., density and temperature) are also set up in the same way as before. Assuming Γ = const = Γ_{e} for the region r ∈ [R_{max}, ∞), we may integrate the hydrostatic equations to estimate the temperature and the density at the outer boundary r = R_{max}. This fixes the effective temperature and the spherically modified optical depth scale of these hydrostatic (HS) models, requiring that the outer boundary always satisfies . Thus in this model, the stellar photospheric radius is defined as (23)
for the corresponding hydrostatic effective temperature (24)
quite analogous to how these quantities typically are defined in static stellar structure and evolution models.
These boundary conditions are used to numerically solve the new set of equations from the outer boundary inwards. A consistent solution is obtained by shooting for the fixed inner boundary temperature and radius, varying the stellar photospheric radius and assuming innerboundary values as discussed above. In this way we obtain static envelopes corresponding directly to our previous dynamic models. The right panel of Fig. 4 shows the predicted densities for static models with different assumed^{5} Γ_{e} ∈ [0.37, 0.40, 0.41] (for comparison see Fig. 1 in Gräfener et al. 2012). These densities are then compared to those of the steadystate computation in Sect. 3, for which Γ_{e} = 0.40.
We note directly that in the static case an extended region of low density is formed. The photospheric radius of this inflated envelope is extremely sensitive to Γ_{e}, with growing notably as Γ_{e} is increased. For the static model with highest assumed Γ_{e} = 0.41, we already find a large accompanied by a relatively cool effective temperature T_{HS, eff} = 21 kK; further increasing Γ_{e} causes the inferred photospheric radius to inflate to increasingly unrealistic, extremely large values. For the Γ_{e} = 0.40 that we use in our dynamic models we find and T_{HS, eff} = 56 kK. The reason for the inflation here is the same as in the models by Petrovic et al. (2006) and Gräfener et al. (2012). Basically, it is related to the fact that for high values of Γ_{e} the integration moves along a Γ ≈ 1 curve, forcing the density to remain low and almost constant over a large radial extent. Since enough column mass (or equivalently optical depth) still has to be accumulated in order to meet the lower boundary condition, the only way the star can react is by expanding in radius (see also discussion in Owocki 2014). Comparing this now to our dynamic solutions shows that, for the same Γ_{e}, inclusion of the dynamical terms reduces the inflation somewhat. Nonetheless, R_{ph} still lies very far away from the core in such outflow models of dynamically inflated envelopes; indeed, for our standard Γ_{e} = 0.40 we obtain R_{ph} = 4.48 R_{⊙} and T_{ph, eff} = 62 kK. Such a close agreement of photospheric radii in dynamic and static models is however not generally found. Namely, a very small increase to Γ_{e} = 0.41 in the static model leads to an order of magnitude increase of the star’s photospheric radius (see Fig. 4), in contrast to the dynamic model where this barely affects the location of R_{ph}.
The analysis above assumes that all luminosity is carried by radiation. In stellar models, the strong increase in opacity around the iron bump will typically lead to an onset of convection, such that L = L_{rad} + L_{conv}. In principle, this may alleviate the inflation by allowing a portion of the luminosity to be carried by convection, lowering the effective Γ_{e} ~ L_{rad}. However, for the conditions prevailing in WR stars convection in these nearsurface layers should be very inefficient. Following Gräfener et al. (2012), an upper limit to the fraction of the total stellar flux F that can be carried by convection is: (25)
This means that for the T_{c, eff} > 10^{5} kK considered in this paper, less than a percent of the total luminosity will be carried by convection.
This convective inefficiency is confirmed by the computation of WR stars using the stellar structure and evolution code package MESA (Paxton et al. 2019). Namely, using the standard mixinglength theory (MLT) prescription for convection shows that such stellar structure and evolution computations^{6} indeed undergo similar inflation in the WR stage as seen here in Fig. 4. However, in these stellar evolution simulations envelope inflation and the associated density inversions often imply prohibitively short time steps within the model calculations. To prevent this from happening, such stellar models often invoke an additional energy transport that is far greater than that implied by the upperlimit estimate of L_{conv} above (see, for example, the socalled MLT++ prescription in MESA, Sect. 7.2 in Paxton et al. 2013). This then forces L_{rad} to always be low enough such that the calculation can proceed.
To mimic these tricks, in the next step we artificially reduce Γ_{e} ~ L_{rad} by introducing L = L_{conv} + L_{rad}. Then using L_{rad} in Eq. (20), while computing the radiation pressure from the total luminosity in Eq. (21), we solve the differential equations exactly as above. By then assuming different ratios L∕L_{rad} for the same fixed total luminosity, we can directly compare dynamic and static inflation models with and without such a reduction in radiative luminosity. Analogous to the left panel, the right panel of Fig. 4 shows the radial density structure of a case with reduced radiative luminosity L∕L_{rad} = 4 (solid line), comparing this to the original case L∕L_{rad} = 1 (solid, marked line). The ratio of L∕L_{rad} = 4 here is chosen such that the stellar envelope no longer experiences a density inversion. The figure demonstrates clearly that the two solutions have vastly different radial scales; while the original model had and T_{HS, eff} = 56 kK, the model with reduced radiative luminosity has a very hot T_{HS, eff} = 121 kK and a photospheric radius that lies less than a percent above the core R_{c} = R_{⊙}.
The reason for these drastically different stellar photospheric radii is the same as described earlier; including an additional (artificial) energy transport reduces the effective Γ_{e} ~ L_{rad}, so that a dense exponential atmosphere with small scale height is formed. This then avoids inflation and density inversions, and the upshot is a hot and compact WR stellar surface. This is a good example of how envelope inflation and density inversions are often avoided in stellar structure computations by reducing L_{rad} well beyond the limits implied by standard convective energy transport.
A key result of our analysis here is thus that such tricks of the nearsurface regions in WR stars are necessary only because of enforcing a hydrostatic solution. As shown above, when including the dynamical terms the stellar envelope will instead quite naturally develop a dynamically inflated outflow region.
Fig. 3 Solutions of a steadystate model with stellar parameters as in Table 2. (a) Temperature profile. (b) Velocity profile in solid line, the region of stagnation or negative acceleration in X markers. (c) Eddington ratios for full radiation force in solid line, only the OPAL contribution in dotted line, only CAK contribution in dashed line, and Γ = 1 in thin horizontal dashed line. (d) Profile of spherically modified optical depth. (e) Density Profile. Panels d,e: the vertical dotted line at R_{ph} = 4.48R_{c} = 4.48R_{⊙} shows the location of the photosphere. The top axis gives the radial distance in units of R_{⊙}, but since in this paper R_{c} = R_{⊙} this is also equivalent to normalising the radius to R_{c}. 
Fig. 4 Density profiles of different static and steadystate solutions. On both plots the solid lines show the different cases of hydrostatic solutions and dashed lines show the subphotospheric part of the steadystate outflow (see Fig. 3e). The different markers then compare hydrostatic solutions for different assumed Thompson scattering Eddington factor Γ_{e} (left panel) and for different assumed L∕L_{rad} ratios for thefixed total luminosity log_{10}(L∕L_{⊙}) = 5.416 (right panel). The numbers at the right end of each density profile correspond to the photospheric effective temperature of the corresponding model in kiloKelvins, as defined in the text. To better visualise the innermost parts of the simulations, the right panel displays the logarithm of the scaled radius x = 1 − R_{c}∕r on the abscissa. 
4 Timedependent numerical hydrodynamic simulations
The previous section presented calculations in the steadystate limit, using various degrees of approximations for solving the dynamical (and static) equations. Building on the same basic formalism, this section now presents full timedependent numerical hydrodynamics simulations for spherically symmetric WR wind outflows. A key objective is to examine whether such timedependent simulations relax to similar solutions as the simplified steady models of Sect. 3.
The hydrodynamic simulations are performed using the finite volume code MPIAMRVAC (Xia et al. 2018) to solve Eqs. (1), (4), and (7). We use the HLL solver (Harten 1983) with a MINMOD flux limiter (van Leer 1979). The radiative acceleration is computed as in Sect. 3.3 by summing up the contributions from OPAL and CAKlike opacities, Γ_{tot} = Γ_{OPAL} + Γ_{CAK}. Here the OPAL contribution is trivially obtained from the local density and temperature, and the CAK contribution from a simple central finite difference scheme applied to Eq. (18). To avoid information propagating over multiple cells in a single integration we further limit the integration time step d t by:
where dt_{CFL} is the time step set by the standard CourantFriedrichsLewy condition (Courant et al. 1928), and Δr the width of one numerical cell.
As Eq. (7) is not the standard way of computing the energy balance in MPIAMRVAC we have implemented an additional routine to handle this. In this subroutine the optical depth Eq. (8) is first computed by trapezoid integration starting from the outer boundary i = i_{max}, setting τ(i_{max}) = τ_{out} as before, inwards to i = i_{min}. The radial temperaturestructure T_{i} is then updated at each hydrodynamical time step according to Eq. (7). For simplicity, we here define R_{c} to always be at the lower boundary of the simulation.
Fig. 5 Surface plot of velocity as a function of time and reduced radius x for the model as described in text. The figure demonstrates how the velocity profile relaxes from the initially assumed v_{β} to (quasi) steadystate profile over ≈40 ks. 
4.1 Initial and boundary conditions
Analogous to our steadystate simulations, we define the outer boundary to be located at some radius R_{max} and use the same optical depth and temperature conditions as described in the previous section, Eq. (7) and the prespecifiedminimum temperature (see above). As in, for instance, Driessen et al. (2019), the numerical outer boundary conditions on the velocity and density are set by simple extrapolation from the outer most points of the simulation grid. To properly resolve the sonic point we set the inner boundary to a lower radius than the expected sonic point R_{a}, i.e. now R_{min} ≡ R_{c} = 0.99 R_{⊙} < R_{a}. At this R_{min} we then fix the lower boundary density ρ_{⋆} and allow the velocity to adjust to the overlying wind conditions (see, e.g. Owocki et al. 1988). ρ_{⋆} is initially estimated from our corresponding steadystate models by simply assuming hydrostatic equilibrium below the sonic point. The stellar mass and luminosity are still as in Table 1.
The initial wind condition is an analytic outflow structure, in the form of a simple βlaw: (26)
where the terminal velocity , and β = 0.5 are derived from a simple assumed Γ = const > 1. The parameter b sets the minimum initial velocity v_{min} at r = R_{c}, here chosen such that v_{min} ≈ 50 km s^{−1}. The equation of continuity then directly provides the initial density structure. Finally, the initial temperature structure is computed from Eq. (7) with the opacity set according to the assumed input Γ from above. Using these initial and boundary conditions, numerical time integration is carried out until the wind relaxes to a (quasi) steady state. We, however, note that the results of our simulations are not sensitive to these assumed initial conditions.
4.2 Timedependent dynamic inflation model
Starting from the initially assumed structure the wind outflow gradually relaxes to a hydrodynamically consistent structure. The velocity relaxation is displayed in Fig. 5 as a surface plot, where along one horizontal axis the physical time t in ks of the simulation is shown and along the other the radius in x = 1 − R_{c}∕r. The vertical axis then measures velocity (in km s^{−1}). At t = 0 the initial velocity structure, v_{β}, is visible, however as time increases this βtype structure gradually evolves towards the threephased velocity profile discussed in Sect. 3.5. Around at t ≈ 40 ks the simulation has relaxed to a state very similar to that found in our corresponding steadystate model (Sect. 3.5).
Running the simulation for in total 600 ks, Figs. 6a and b show the Eddington ratio Γ and the velocity profile of the final time step. From Fig. 6a we observe that, just as in the steadystate model, the CAKlike force near the lower boundary on average only has a small contribution to the total force; contrarily, near the outer boundary, itis dominant. The radiation force from the OPAL contribution has the opposite behaviour, dominating at the sonic point and inthe inner regions, but becoming secondary in the outer parts. This then again leads to a velocity profile (solid line in Fig. 6b) with a nonmonotonic behaviour, which clearly can not be fit with the simple βtype velocity law assumed in the initial conditions (dotted line).
This velocity profile can be directly compared to that found in our corresponding steadystate model, illustrated by the dashed line in Fig. 6b, which shows close agreement between the two. The mass loss of the timedependent model is Ṁ = 1.64 × 10^{−5} M_{⊙} yr^{−1}, which also agrees well with the massloss rate of Ṁ = 1.47 × 10^{−5} M_{⊙} yr^{−1} from our steadystate model (see Sect. 3). In this way, our timedependent models demonstrate that the dynamicallyinflated solution discussed in the previous section indeed works as a stable attractor for the stellar wind envelope, confirming our basic conclusions regarding a wind outflow controlled primarily by OPAL opacities in the inner regions, and by the CAKlike force in the outer. The small differences in velocity magnitude and mass loss between steadystate and timedependent models are mainly related to differences in setting the lower boundary. In steadystate computations we define the sonic point to be exactly at the lower boundary, whereas in timedependent computations it is allowed to self adjust so that the subsonic part is also resolved.
Fig. 6 Profiles of final time step from the timedependent numerical model (see Sect. 4.2). (a) Eddington ratio, with the definition of line styles the same as on Fig. 3c. (b) Comparison of the initial (dotted line), final time step (solid line) and steadystate (dashed line) velocity fields. 
4.3 Ṁ dependency on CAK α
A general finding of our baseline simulations in Sects. 3 and 4 is that (as long as the wind can escape), the massloss rate is primarily controlled by the conditions in the lower wind; in contrast, the terminal wind speed is controlled mainly by the outer wind conditions. Moreover, while it has turned out that a steadystate model, such as the one presented inSect. 3, is quite difficult to converge, the timedependent simulations of this section are computationally somewhat easier to handle.
Therefore, we have used our timedependent setup to further explore the impact of the CAKlike force parameters on the predictions of Ṁ and v_{∞}. We investigate the impact of these parameters by varying only α_{min} and α_{max}, however we point out that comparable effects could also be achieved by instead introducing corresponding radial variations of the second lineforce parameter . To this end, we test a range of values α_{max} ∈ [0.60, 0.66] and α_{min} ∈ [0.48, 0.53], corresponding roughly to the scatter of α observed in the MC simulations displayed in Fig. 2. The results of these tests are summarised in Tables 3 and 4, along with the massloss rates and terminal velocities found for the different model runs. In addition, we also tested variations in the spatial position of the α_{max} to α_{min} transition for the standard CAK force parameters in the range x ∈ [0.3, 0.7]. However, for these tests, no significant influence on the massloss rates of the models have been observed. Therefore, the following subsections focus exclusively on varying the values of α_{max} and α_{min}.
Summary of massloss rates and terminal wind speeds for models with different values of α_{min}.
Summary of massloss rates and terminal wind speeds for different values of α_{max}.
4.3.1 Dependency on α_{min}
As Γ_{CAK} is the prevalent driver in the outer wind, one might initially expect a similar dependence of the massloss rate on α_{min} as given by the standard CAK model. The CAK expression for the massloss rate Ṁ_{CAK}, typically derived from criticalpoint analyses for Γ_{tot} = Γ_{e} + Γ_{CAK} in the radial streaming limit, is: (27)
From this expression, it directly follows that Ṁ_{CAK} increases steeply if α is decreased while keeping all other parameters constant. However, it can be seen from Table 3 that in our models Ṁ varies only marginally across the range of considered α_{min} values. This distinguishes these optically thick WR simulations from standard CAKtype models applied to the optically thin winds of Otype stars. We further confirm the finding from the steady models that if α_{min} is increased too much, the line force becomes insufficient to sustain the outflow. This happens because the line acceleration in the outer wind must be sufficient to prevent fallback of the mass loss initiated at the Feopacity bump. That is, the line acceleration must provide enough additional energy such that the previously introduced condition is fulfilled. As alluded to above, the CAK massloss rate formula cannot be directly used to predict the mass loss of this hybrid model. Nevertheless, the limiting value of α_{min} that can sustain the wind can be estimated by inverting the above equation, solving for the α at which the CAK massloss rate equals that of our actual model. For the given model with Ṁ = 1.64 × 10^{−5} M_{⊙} yr^{−1}, one finds α_{min} ≈ 0.52. The results of our numerical simulations, summarised in Table 3, indeed suggest that this provides a reasonable estimate for the limiting α_{min} value.
This is also supported by closer inspection of the simulations with α_{min} = 0.52, which results in a flattening of the wind velocity profile as . With a further increase to α_{min} = 0.53, we find thatthe CAKlike force is unable to keep up with the mass flux initiated at the lower boundary. In the 1D simulations presented here, this then leads to a similar failed wind situation as discussed in the previous section; we discuss potential implications of this behaviour for multidimensional versions of our simulations at the end of the next section.
Finally, as summarised in Table 3 the terminal velocity indeed is very sensitive to small changes of α_{min}. Overall, this thus suggests that within this formalism, Ṁ and v_{∞} are determined quite independently. As discussed further in Sect. 5.2, such a basic decoupling of Ṁ and v_{∞} indeed seems to be supported by observational data of WR winds (e.g. see Hamann et al. 2019; Sander et al. 2019). On the other hand, this result also seems to contradict that of the most recent study from Sander & Vink (2020), who find a correlation between Ṁ and v_{∞}. We further discuss possible reasons for this in Sect. 5.4.
4.3.2 Dependency on α_{max}
In contrast to α_{min}, α_{max} has little to no effect on determining the terminal velocity of the outflow. Nevertheless, it can affect Ṁ. This is illustrated in Table 4, where a modest increase in Ṁ is found for decreasing values of α_{max}. We note, however, that again the magnitude of the mass loss variation is small in comparison to that seen in CAKtype wind models of optically thin Ostars.
This rather weak dependence on the CAK force stems from the fact that the dominant contribution to the total radiation force comes from the OPAL opacities in the regions near the stellar core (see also Sect. 3). At the sonic point for α_{max} = 0.6, we find Γ_{CAK} ≈ 0.14 and for the α_{max} = 0.66 we find Γ = 0.02, which is much lower even than our initial order of magnitude estimate Γ ≈ 0.1 in Sect. 3.4. This further supports that the CAK force has a negligible contribution at the sonic point. Moreover, since those inner parts accumulate the majority of the integrated wind optical depth, they also primarily control the lowerboundary temperature to which the overlying wind adjusts, and thereby also the massloss rate.
5 Discussion
5.1 Implication for core radius problem
The dynamic models presented in Sects. 3 and 4 provide a natural explanation to the socalled core radius problem of classical WR stars (see Sect. 1). To understand this, we introduce the standard Rosseland mean optical depth scale (28)
where κ_{Ros} is computed using the Rosseland mean opacities from the OPAL tabulation, such that κ_{Ros} = κ_{OPAL} throughout this work.
Then, following the standard spectroscopic model approach of fitting an assumed v_{β} velocity profile to the entire supersonic region, Fig. 7 directly illustrates how the hydrostatic core (where v_{β} reaches subsonic velocities) in that case would be located at a radius . Indeed, the core radius recovered in this way is a factor ≈2.6 times larger than the true value, demonstrating the basic challenge of determining WR hydrostatic radii from spectroscopic data (e.g. Crowther 2007). The discrepancy essentially stems from the fact that τ_{Ros} = 2∕3 is reached already at very high radii (see the grey area in Fig. 7), so that only the outermost velocity law is visible to us. Observational estimates of the hydrostatic stellar core radius must thus be obtained by some extrapolation of the velocity law assumed in the spectroscopic model, without any knowledge of the true velocity field in the invisible inner regions.
The dynamic models presented in this paper explain this longstanding discrepancy simply by their extended (dynamically inflated) region of lowvelocity outflowing material above the hydrostatic stellar core.
To further illustrate this effect, let us inspect the observational positioning of classical WR stars on the Hertzsprung–Russell diagram (HRD), where the stars are placed according to their observationally inferred luminosities and hydrostatic core effective temperatures. Figure 8 here shows the HRD where the observed galactic sample of WR stars by Hamann et al. (2019) and Sander et al. (2019) is placed. Positions of these stars on the HRD can then be compared to the location of the socalled helium zero age main sequence (HeZAMS; dashed line in figure), computed here with MESA using the standard mixinglength theory for convection with solar metallicity and relative metal fractions as given by Asplund et al. (2009)^{7}. On average, observationally inferred positions of WR stars are then expected to lay to the left of this HeZAMS line. This is because typical stellar evolution models, where an additional energy transport is typically also invoked in order to prevent the WR cores from inflating (see Sect. 3.5), occupy the region to the hotter side of HeZAMS on the HRD. Figure 8, however, shows that the average position of these WR stars (indicated by + in the figure) is located on the cooler side of the HeZAMS instead. To improve the agreement between spectroscopy and evolution predictions a simple correction factor ≈ 3 can be applied to the spectroscopic radii, which then shifts the core effective temperature location (now identified on Fig. 8 by x). And indeed, such a correction factor matches quite well with the correction factor implied by our dynamic inflation models above.
Fig. 7 Steadystate velocity profile (solid) of the stellar model used in this paper. The grey shading marks the region of high photospheric continuum optical depth τ_{Ros} ≥2∕3, with R(τ_{Ros} = 2∕3) ≈ 3.5R_{⋆}. In crosseddashed line a v_{β} velocity fit to the optically thin region is given, which was extrapolated to the optically thick region to recover the stellar core radius. 
Fig. 8 Hertzsprung–Russell diagram including galactic WR stars of nitrogen (WN), carbon (WC), and oxygen sequences (WO). The solid shapes mark inferred core effective temperatures and luminosities from Hamann et al. (2019; WN) and Sander et al. (2019; WO,WC); the mean value is marked with a +. The x marks the location of the modified average luminosity and temperature, where the spectroscopically inferred core radius has been reduced by a factor of 3 (see text). The grey dashed line shows the approximate predicted location of the Helium Zero Age Main Sequence (HeZAMS) computed using MESA (see the text). 
Fig. 9 Distribution of massloss rates and terminal velocities derived from the spectroscopic fitting by Sander et al. (2019) and Hamann et al. (2019) for an observed stellar population of galactic WR stars of nitrogen (WN), carbon (WC), and oxygen sequences (WO). Massloss rates and terminal velocities derived in our steadystate (log_{10} (Ṁ∕ M_{⊙} yr^{−1}) = −4.83, v_{∞} = 2200 km s^{−1}) and timedependent (log_{10}(Ṁ∕ M_{⊙} yr^{−1}) = −4.78, v_{∞} = 2063 km s^{−1}) simulations are consistent with the range of mass loss shown on this figure. 
Summary of derived parameters of our steadystate and timedependent models.
5.2 Empirical constraints and model dependencies of mass loss and terminal velocity
The basic model presented in the previous sections shows an overall quite good agreement with empirical constraints on terminal velocities and massloss rates of WR stars. Figure 9 plots such empirically inferred values of Ṁ and v_{∞}, again from Hamann et al. (2019) and Sander et al. (2019). Compared to the massloss rates and terminal velocities predicted by our hydrodynamic and steadystate models (summarised in Table 5), the latter is located at the high end of the inferred values. However, since v_{∞} is primarily determined by the line force in the outer wind regions, which here is parametrised by the simple radial streaming CAKapproximation outlined in previous sections, a more detailed approach will be necessary to make more quantitative comparisons to empirical v_{∞}. Inspection of Fig. 9 also reveals that the empirical massloss rates and terminal wind speeds show no prominent anticorrelation. This further supports another basic characteristic of our models, namely that Ṁ and v_{∞} are quite decoupled (see Sect. 4.3.1).
5.3 Dependence on metallicity
The model of dynamic inflation presented in this paper relies on the fact that the Γ =1 condition is reached already at the Feopacity bump. As such it is sensitive to the variation of opacity at the Feopacity bump, which is directly controlled by metallicity (e.g. Cantiello et al. 2009). Qualitatively, the model of dynamic inflation should respond to a decrease or increase in metallicity by reducing or enlarging the scale of inflation and the associated mass loss. Moreover, at low enough metallicity the Γ = 1 conditionwill not be fulfilled at the Fe bump, in which case dynamic inflation will no longer take place from the hot Fe bump, and a standard optically thin wind develops instead. This behaviour is also in line with the recently published models by Sander & Vink (2020).
5.4 Comparison to comoving frame models
As mentioned above, recently Sander et al. (2020) and Sander & Vink (2020) published timeindependent wind simulations of classical WR stars. These models solve the same steadystate e.o.m. as here^{8}, however using a different method to calculate the radiative acceleration. There Γ(r) is derived from full solutions of the frequencydependent radiative transfer equation in the comoving frame (CMF), without the use of any parameterisation to estimate the lineforce contribution (see also Sundqvist et al. 2019; Björklund et al. 2021, for similar approaches applied to Ostars). Although these models thus compute Γ(r) for a given velocity and density structure in a more detailed way than here, this CMF method (as currently implemented) has one major disadvantage. Because of the way the CMF transfer equation is numerically solved, it can only be applied to monotonically increasing (or decreasing) velocity fields v(r). As such, this method cannot be used to derive Γ(r) for the type of nonmonotonic velocity structures found here (see Figs. 3b, 6b).
The way Sander et al. (2020) overcome this issue is by assuming an ad hoc very high degree of clumping in the WR outflow (see also Gräfener & Hamann 2005). They parameterise this using a simple model where all wind mass is assumed to be concentrated within overdense clumps, such that the densities of these clumps are ρ_{cl} = ⟨ρ⟩∕f_{vol} = D⟨ρ⟩, where ⟨ρ⟩ is the mean density (assumed to be preserved with respect to a smooth model), f_{vol} the volume fraction of the total wind contained by clumps, and D the clumpoverdensity factor.
By adopting a very high D = 50 (f_{vol} = 0.02), the ionisation balance of the wind is shifted so that more effective driving lines become available, and the line force is thereby increased enough to avoid any deceleration regions. When assuming instead a more modest D = 10 in their models, they do experience the same type ofΓ < 1 regions as here, which leads then to a nonmonotonic v(r) that cannot be handled by the CMF radiative transfer (see Fig. 1 in Sander et al. 2020, and their corresponding discussion). Since observations of electron scattering wings in WR outflows suggest quite modest values D ≈ 4−10 (see overviews by Crowther 2007; Puls et al. 2008), and theoretical windclumping models following the linedriven instability (Owocki et al. 1988) have so far only been developed for OBstars (with again typical values D ≈ 10; Sundqvist et al. 2018; Driessen et al. 2019), it is at present not clear how the high clump densities needed to keep v(r) monotonic in CMFbased WR models might be physically justified.
The models presented here are thus very complementary to such CMFbased simulations. Although our more simplified treatment of Γ involves certain approximations (as discussed in previous sections) it allows us to model regions of stagnation and deceleration (see Sects. 3 and 4), whereas such nonmonotonic flows are incompatible with the methods currently used in CMFbased calculations. Moreover, the method presented here allows for timedependent simulations that otherwise would not be computationally feasible. In this respect, we also recall our finding that if the line force in the outer wind is decreased too much, the wind not only experiences regions of deceleration but may eventually reach zero velocity. In the 1D simulations presented here, this leads to failed wind solutions and model termination. In a multidimensional simulation, one may speculate that it might instead result in a very complex wind structure of coexisting regions of upflows and inflows. In turn, this may then cause shocks and perhaps increased levels of clump formation. These could be simulated in multidimensional hydrodynamics simulations using the method we present here.
The models presented here also demonstrate that even with multiple extended regions where Γ <1, the accumulated net energy of the wind can still be positive (i.e., we can still have ). This suggests that the WR outflow can be driven without invoking such a high degree of clumping, but instead allowing for nonmonotonic velocities. For example, let us consider the lower panel of Fig. 1 in Sander et al. (2020), where the authors provide Γ(r) not only for their final D = 50 model (their solid line) but also for their collapsed D = 10 solution (their dashed line). From these profiles, we estimated the net energy corresponding to both models and found that, indeed, both of the force configurations provide . As such, both solutions may actually escape the stellar potential, but since the D = 10 model would have contained regions of negative acceleration it could not be treated by the assumed CMF radiative transfer framework.
Note that Sander et al. (2020) make the argument that although the assumed very high degree of clumping influences v_{∞} significantly, it does not have a major impact on Ṁ, since clumping is assumed to be present only at supersonic velocities. This is quite consistent with our finding here. On the other hand, Sander & Vink (2020) also find a correlation between Ṁ and v_{∞}, which we do not necessarily observe in our models. This result, however, might simply be related to them forcing a monotonic velocity field in their simulations, since the same correlation was also found in the βlaw WR models by Gräfener et al. (2017).
6 Summary and future work
This paper has presented a model for the wind outflows of classical WolfRayet stars, based on a hybrid opacity method where OPAL opacities are used in combination with a CAKlike line force parametrisation. The latter is used here to account for the effects of Doppler shifts upon spectral line opacity in an approximate way. Our hybrid opacity model essentially results in a two stage WR stellar wind: a supersonic outflow is first initiated in optically thick layers (at the socalled Fe bump) by the OPAL opacities. However, these opacities are not able to prevent fallback upon the core, and so it is the CAKlike force that takes over the driving in the outer regions and ensures that the initiated outflow can also escape the stellar potential.
Our approach leads to deepseated wind initiation with high massloss rates and terminal velocities (see Table 5), in the range of values typically inferred for classical WR stars by spectroscopic analyses. The simulations display dense and slow – dynamically inflated – subphotospheric layers, offering a natural explanation to the socalled coreradius problem of classical WR stars. Direct comparison with hydrostatic computations further shows that neglecting the dynamical terms in the e.o.m. instead leads to drastic inflation of the then assumed static envelope. In this way, we demonstrateexplicitly that the need to invoke huge amounts of ad hoc energy transport in these convectively inefficient regions of stellar structure and evolution calculations is simply an artefact of neglecting the dynamic outer boundary.
The paper here uses a lineforce parameterisation that introduces a very simple radial variation to the CAK powerlaw index α in order to sustain the outer wind outflow. The quite simple nature of this basic setup means that it should be rather straightforward to extend these first simulations towards exploration of a more complete stellar parameter space for classical WR stars. Since the assumed CAK lineforce parameterisation proved to have little effect on the resulting mass loss (as long as we ensure a positive wind energy), such model grids might then offer good alternatives for direct inclusion of massloss rates into stellar evolution models.
The assumed treatment of α does influence the outer wind velocity profile and terminal speed. As such, constraining the morphology of the reaccelerating part of the outflow will require more detailed studies. This is a very challenging task, however, since a nonparametrised calculation of the line force in these regions must be able to deal with nonmonotonic velocity fields. As such, current CMF radiative transfer methods, like those applied in some alternative WR (Gräfener & Hamann 2005; Sander et al. 2020; Sander & Vink 2020) and Otype star wind models (Sundqvist et al. 2019; Björklund et al. 2021), will not be sufficient.
The assumed variation of α can also be varied in multidimensional extensions of the timedependent simulations presented here. In such a study we may allow parts of the wind material to fall back into the deep regions of the outflow. While in 1D simulations such fallback eventually always leads to model termination, in a multidimensional setting it may instead create a highly structured, turbulent flow. This might then still deposit enough momentum into the plasma that a wind can be sustained, but now with clumps of material perhaps mechanically transporting some of the momentum through collisions.
As the wind velocity field also affects spectral line formation, another important follow up regards analysis of spectral features associated with the nonmonotonic velocities found in this paper. In particular, a key question becomes if there are any fundamental differences in the predicted strengths and shapes of strategic spectral lines calculated from our models as compared to those predicted by standard spectroscopic models assuming a monotonic βtype velocity law (e.g. cmfgen, Hillier & Miller 1998). Using the newly developed 3D radiative transfer code by Hennicker et al. (2020), studies of this are already underway and will be presented in an upcoming paper. Finally, a very interesting problem to address concerns stellar structure and evolution, and the manifestation of dynamic inflation in different evolutionary stages. Such a study might be especially important since quite many stars (some already at earlier evolution stages) reach the Eddington limit at opacity bumps in regions that might have insufficient convective energy transport. As demonstrated in Sect. 3, dynamically inflated envelopes should then develop, which in turn could have a significant feedback on the stellar structure as well as on the further evolution of the star.
Acknowledgements
The authors would like to thank J. Puls, T. Shenar and all members of the KUL EQUATION group for fruitful discussion, comments, and suggestions. L.P., J.S., N.K., L.D., A.d.K., L.M., and H.S. acknowledge support from the KU Leuven C1 grant MAESTRO C16/17/007. L.P. further acknowledges the additional support from Shota Rustaveli Georgian National Foundation Grant Project No. FR17391. J.S. further acknowledges additional support by the Belgian Research Foundation Flanders (FWO) Odysseus program under grant number G0H9218N.
References
 Asplund, M., Grevesse, N., Sauval, A. J., & Scott, P. 2009, ARA&A, 47, 481 [NASA ADS] [CrossRef] [Google Scholar]
 Beals, C. S. 1929, MNRAS, 90, 202 [Google Scholar]
 Björklund, R., Sundqvist, J. O., Puls, J., & Najarro, F. 2021, A&A, in press, https://www.aanda.org/10.1051/00046361/202038384 [Google Scholar]
 Cantiello, M., Langer, N., Brott, I., et al. 2009, A&A, 499, 279 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Cassinelli, J. P. 1991, IAU Symp., 143, 289 [Google Scholar]
 Castor, J. I., Abbott, D. C., & Klein, R. I. 1975, ApJ, 195, 157 [NASA ADS] [CrossRef] [Google Scholar]
 Courant, R., Friedrichs, K., & Lewy, H. 1928, Math. Ann., 100, 32 [NASA ADS] [CrossRef] [MathSciNet] [Google Scholar]
 Crowther, P. A. 2007, ARA&A, 45, 177 [NASA ADS] [CrossRef] [Google Scholar]
 Crowther, P. A., Schnurr, O., Hirschi, R., et al. 2010, MNRAS, 408, 731 [NASA ADS] [CrossRef] [Google Scholar]
 de Koter, A., Heap, S. R., & Hubeny, I. 1997, ApJ, 477, 792 [Google Scholar]
 Doran, E. I., Crowther, P. A., de Koter, A., et al. 2013, A&A, 558, A134 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Driessen, F. A., Sundqvist, J. O., & Kee, N. D. 2019, A&A, 631, A172 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Ekström, S., Georgy, C., Eggenberger, P., et al. 2012, A&A, 537, A146 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Gayley, K. G. 1995, ApJ, 454, 410 [NASA ADS] [CrossRef] [Google Scholar]
 Gayley, K. G., Owocki, S. P., & Cranmer, S. R. 1995, ApJ, 442, 296 [NASA ADS] [CrossRef] [Google Scholar]
 Gräfener, G.,& Hamann, W. R. 2005, A&A, 432, 633 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Gräfener, G., Owocki, S. P., & Vink, J. S. 2012, A&A, 538, A40 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Gräfener, G., Owocki, S. P., Grassitelli, L., & Langer, N. 2017, A&A, 608, A34 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Grassitelli, L., Langer, N., Grin, N. J., et al. 2018, A&A, 614, A86 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Grevesse, N., & Noels, A. 1993, Phys. Scr. T, 47, 133 [Google Scholar]
 Groh, J. H., Georgy, C., & Ekström, S. 2013, A&A, 558, L1 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Hamann, W. R., Gräfener, G., Liermann, A., et al. 2019, A&A, 625, A57 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Harten, A. 1983, J. Comput. Phys., 49, 357 [NASA ADS] [CrossRef] [MathSciNet] [Google Scholar]
 Hennicker, L., Puls, J., Kee, N. D., & Sundqvist, J. O. 2020, A&A, 633, A16 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Hillier, D. J., & Miller, D. L. 1998, ApJ, 496, 407 [NASA ADS] [CrossRef] [Google Scholar]
 Iglesias, C. A., & Rogers, F. J. 1996, ApJ, 464, 943 [NASA ADS] [CrossRef] [Google Scholar]
 Ishii, M., Ueno, M., & Kato, M. 1999, PASJ, 51, 417 [NASA ADS] [Google Scholar]
 Kee, N. D., Owocki, S., & Sundqvist, J. O. 2016, MNRAS, 458, 2323 [NASA ADS] [CrossRef] [Google Scholar]
 Lamers, H. J. G. L. M., & Leitherer, C. 1993, ApJ, 412, 771 [NASA ADS] [CrossRef] [Google Scholar]
 Lucy, L. B. 1971, ApJ, 163, 95 [NASA ADS] [CrossRef] [Google Scholar]
 Lucy, L. B., & Abbott, D. C. 1993, ApJ, 405, 738 [NASA ADS] [CrossRef] [Google Scholar]
 Lucy, L. B., & Solomon, P. M. 1970, ApJ, 159, 879 [NASA ADS] [CrossRef] [Google Scholar]
 Nugis, T., & Lamers, H. J. G. L. M. 2002, A&A, 389, 162 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Owocki, S. 2014, ArXiv eprints [arXiv:1409.2084] [Google Scholar]
 Owocki, S. P., Castor, J. I., & Rybicki, G. B. 1988, ApJ, 335, 914 [Google Scholar]
 Owocki, S. P., Townsend, R. H. D., & Quataert, E. 2017, MNRAS, 472, 3749 [Google Scholar]
 Pauldrach, A. W. A., Feldmeier, A., Puls, J., & Kudritzki, R. P. 1993, Space Sci. Rev., 66, 105 [NASA ADS] [CrossRef] [Google Scholar]
 Paxton, B., Cantiello, M., Arras, P., et al. 2013, ApJS, 208, 4 [NASA ADS] [CrossRef] [Google Scholar]
 Paxton, B., Smolec, R., Schwab, J., et al. 2019, ApJS, 243, 10 [Google Scholar]
 Petrenz, P., & Puls, J. 2000, A&A, 358, 956 [NASA ADS] [Google Scholar]
 Petrovic, J., Pols, O., & Langer, N. 2006, A&A, 450, 219 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Prantzos, N., Abia, C., Limongi, M., Chieffi, A., & Cristallo, S. 2018, MNRAS, 476, 3432 [Google Scholar]
 Proga, D., Stone, J. M., & Drew, J. E. 1998, MNRAS, 295, 595 [NASA ADS] [CrossRef] [Google Scholar]
 Puls, J., Springmann, U., & Lennon, M. 2000, A&AS, 141, 23 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Puls, J., Urbaneja, M. A., Venero, R., et al. 2005, A&A, 435, 669 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Puls, J., Vink,J. S., & Najarro, F. 2008, A&ARv, 16, 209 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Ramachandran, V., Hamann, W. R., Hainich, R., et al. 2018, A&A, 615, A40 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Ro, S., & Matzner, C. D. 2016, ApJ, 821, 109 [NASA ADS] [CrossRef] [Google Scholar]
 Sander, A. A. C., & Vink, J. S. 2020, MNRAS, 499, 873 [Google Scholar]
 Sander, A. A. C., Hamann, W. R., Todt, H., et al. 2019, A&A, 621, A92 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Sander, A. A. C., Vink, J. S., & Hamann, W. R. 2020, MNRAS, 491, 4406 [Google Scholar]
 Sanyal, D., Grassitelli, L., Langer, N., & Bestenlehner, J. M. 2015, A&A, 580, A20 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Springmann, U. 1994, A&A, 289, 505 [NASA ADS] [Google Scholar]
 Springmann, U., & Puls, J. 1998, ASP Conf. Ser., 131, 286 [Google Scholar]
 Sundqvist, J. O., Owocki, S. P., & Puls, J. 2018, A&A, 611, A17 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Sundqvist, J. O., Björklund, R., Puls, J., & Najarro, F. 2019, A&A, 632, A126 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Todt, H., Peña, M., Hamann, W. R., & Gräfener, G. 2010, A&A, 515, A83 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 van Leer, B. 1979, J. Comput. Phys., 32, 101 [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]
 Wolf, C. J. E., & Rayet, G. 1867, Acad. Sci. Paris CR, 65, 292 [Google Scholar]
 Xia, C., Teunissen, J., El Mellah, I., Chané, E., & Keppens, R. 2018, ApJS, 234, 30 [Google Scholar]
 Yoon, S. C., Gräfener, G., Vink, J. S., Kozyreva, A., & Izzard, R. G. 2012, A&A, 544, L11 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
For example, in Sect. 3.6 the stellar photoshere is located at the τ_{sp} = 4∕3 surface in the radiation diffusion approximation, rather than the typically assumed photoshere location at the τ_{sp} = 2∕3 surface.
The numerical simulations presented in Sects. 3–4 a posteriori confirm the simple order of magnitude estimate made here, showing explicitly that the force contribution from the CAKlike opacities near the lower boundary indeed is very small for our considered cases (see, e.g. Fig. 3c). This thus resolves some potential issues of our basic approach (simple summation of OPAL and CAKlike opacities) regarding double counting the effect of spectral lines in the static limit.
Detailed description of these models as well as Input files to reproduce our MESA simulations are provided at https://doi.org/10.5281/zenodo.4054811
The computed location of the HeZAMS as well as Input files to reproduce our MESA simulations are provided at https://doi.org/10.5281/zenodo.4054811
All Tables
Summary of massloss rates and terminal wind speeds for models with different values of α_{min}.
Summary of massloss rates and terminal wind speeds for different values of α_{max}.
All Figures
Fig. 1 (a) Colour map of Γ_{OPAL} for the stellar parameters from Table 1. Temperature is given on the abscissa and radiationtogas pressure ratio on the ordinate. is identified with a white contour. The colour bar on the right corresponds to the colour coding of the plot and gives the numerical value of Γ_{OPAL} for a given temperature and radiationtogas pressure ratio. Different line styles show the radiationtogas pressure ratio and temperature structure for the different fixed massloss rates. (b) Ratio of kinetic to core potential energy of the corresponding solutions from panel a. The abscissa here shows the radius coordinate x. 

In the text 
Fig. 2 Plot that demonstrates values of α (see text) required to reproduce the line force computed from MC radiative transfer calculations by Springmann (1994; asterisks). The solid line then overplots the piecewise linear function adopted in this paper to describe the corresponding variation of α. The vertical dashed lines correspond to r_{1} (x_{1}) and r_{2} (x_{2}). 

In the text 
Fig. 3 Solutions of a steadystate model with stellar parameters as in Table 2. (a) Temperature profile. (b) Velocity profile in solid line, the region of stagnation or negative acceleration in X markers. (c) Eddington ratios for full radiation force in solid line, only the OPAL contribution in dotted line, only CAK contribution in dashed line, and Γ = 1 in thin horizontal dashed line. (d) Profile of spherically modified optical depth. (e) Density Profile. Panels d,e: the vertical dotted line at R_{ph} = 4.48R_{c} = 4.48R_{⊙} shows the location of the photosphere. The top axis gives the radial distance in units of R_{⊙}, but since in this paper R_{c} = R_{⊙} this is also equivalent to normalising the radius to R_{c}. 

In the text 
Fig. 4 Density profiles of different static and steadystate solutions. On both plots the solid lines show the different cases of hydrostatic solutions and dashed lines show the subphotospheric part of the steadystate outflow (see Fig. 3e). The different markers then compare hydrostatic solutions for different assumed Thompson scattering Eddington factor Γ_{e} (left panel) and for different assumed L∕L_{rad} ratios for thefixed total luminosity log_{10}(L∕L_{⊙}) = 5.416 (right panel). The numbers at the right end of each density profile correspond to the photospheric effective temperature of the corresponding model in kiloKelvins, as defined in the text. To better visualise the innermost parts of the simulations, the right panel displays the logarithm of the scaled radius x = 1 − R_{c}∕r on the abscissa. 

In the text 
Fig. 5 Surface plot of velocity as a function of time and reduced radius x for the model as described in text. The figure demonstrates how the velocity profile relaxes from the initially assumed v_{β} to (quasi) steadystate profile over ≈40 ks. 

In the text 
Fig. 6 Profiles of final time step from the timedependent numerical model (see Sect. 4.2). (a) Eddington ratio, with the definition of line styles the same as on Fig. 3c. (b) Comparison of the initial (dotted line), final time step (solid line) and steadystate (dashed line) velocity fields. 

In the text 
Fig. 7 Steadystate velocity profile (solid) of the stellar model used in this paper. The grey shading marks the region of high photospheric continuum optical depth τ_{Ros} ≥2∕3, with R(τ_{Ros} = 2∕3) ≈ 3.5R_{⋆}. In crosseddashed line a v_{β} velocity fit to the optically thin region is given, which was extrapolated to the optically thick region to recover the stellar core radius. 

In the text 
Fig. 8 Hertzsprung–Russell diagram including galactic WR stars of nitrogen (WN), carbon (WC), and oxygen sequences (WO). The solid shapes mark inferred core effective temperatures and luminosities from Hamann et al. (2019; WN) and Sander et al. (2019; WO,WC); the mean value is marked with a +. The x marks the location of the modified average luminosity and temperature, where the spectroscopically inferred core radius has been reduced by a factor of 3 (see text). The grey dashed line shows the approximate predicted location of the Helium Zero Age Main Sequence (HeZAMS) computed using MESA (see the text). 

In the text 
Fig. 9 Distribution of massloss rates and terminal velocities derived from the spectroscopic fitting by Sander et al. (2019) and Hamann et al. (2019) for an observed stellar population of galactic WR stars of nitrogen (WN), carbon (WC), and oxygen sequences (WO). Massloss rates and terminal velocities derived in our steadystate (log_{10} (Ṁ∕ M_{⊙} yr^{−1}) = −4.83, v_{∞} = 2200 km s^{−1}) and timedependent (log_{10}(Ṁ∕ M_{⊙} yr^{−1}) = −4.78, v_{∞} = 2063 km s^{−1}) simulations are consistent with the range of mass loss shown on this figure. 

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.