Issue 
A&A
Volume 614, June 2018



Article Number  A86  
Number of page(s)  17  
Section  Stellar atmospheres  
DOI  https://doi.org/10.1051/00046361/201731542  
Published online  19 June 2018 
Subsonic structure and optically thick winds from Wolf–Rayet stars
^{1}
ArgelanderInstitut für Astronomie, Universität Bonn,
Auf dem Hügel 71,
53121
Bonn, Germany
email: luca@astro.unibonn.de
^{2}
Dublin Institute for Advanced Studies,
31 Fitzwilliam Place,
Dublin
2, Ireland
^{3}
Centre for AstroParticle Physics and Astrophysics, DIAS Dunsink Observatory,
Dunsink Lane,
Dublin
15, Ireland
^{4}
Department of Physics and Astronomy, University of Sheffield, Hicks Building, Hounsfield Rd, Sheffield, S3 7RH, UK
^{5}
MaxPlanckInstitut für Astronomie,
Königstuhl 17,
69117
Heidelberg, Germany
Received:
10
July
2017
Accepted:
11
February
2018
Mass loss by stellar wind is a key agent in the evolution and spectroscopic appearance of massive main sequence and postmain sequence stars. In Wolf–Rayet stars the winds can be so dense and so optically thick that the photosphere appears in the highly supersonic part of the outflow, veiling the underlying subsonic part of the star, and leaving the initial acceleration of the wind inaccessible to observations. Here we investigate the conditions and the structure of the subsonic part of the outflow of Galactic Wolf–Rayet stars, in particular of the WNE subclass; our focus is on the conditions at the sonic point of their winds. We compute 1D hydrodynamic stellar structure models for massive helium stars adopting outer boundaries at the sonic point. We find that the outflows of our models are accelerated to supersonic velocities by the radiative force from opacity bumps either at temperatures of the order of 200 kK by the iron opacity bump or of the order of 50 kK by the heliumII opacity bump. For a given massloss rate, the diffusion approximation for radiative energy transport allows us to define the temperature gradient based purely on the local thermodynamic conditions. For a given massloss rate, this implies that the conditions in the subsonic part of the outflow are independent from the detailed physical conditions in the supersonic part. Stellar atmosphere calculations can therefore adopt our hydrodynamic models as ab initio input for the subsonic structure. The close proximity to the Eddington limit at the sonic point allows us to construct a sonic HR diagram, relating the sonic point temperature to the luminositytomass ratio and the stellar massloss rate, thereby constraining the sonic point conditions, the subsonic structure, and the stellar wind massloss rates of WNE stars from observations. The minimum stellar wind massloss rate necessary to have the flow accelerated to supersonic velocities by the iron opacity bump is derived. A comparison of the observed parameters of Galactic WNE stars to this minimum massloss rate indicates that these stars have their winds launched to supersonic velocities by the radiation pressure arising from the iron opacity bump. Conversely, stellar models which do not show transonic flows from the iron opacity bump form lowdensity extended envelopes. We derive an analytic criterion for the appearance of envelope inflation and of a density inversion in the outer subphotospheric layers.
Key words: hydrodynamics / stars: atmospheres / stars: massive / stars: Wolf–Rayet / stars: winds, outflows / stars: massloss
© ESO 2018
1 Introduction
A star loses energy from its surface not only in the form of radiation, but also in the form of a stream of particles, the stellar wind (Parker 1958). Stellar winds are outflows of material from a star which play a major role especially in massive star evolution. Stellar winds can in fact lead massive stars to lose a significant fraction of their initial mass, thus enriching the interstellar medium and influencing the structure of the outer layers and the spectroscopic appearance of these stars (Chiosi & Maeder 1986; Lamers et al. 1991; Langer 2012).
The proximity of massive stars to the Eddington limit in the upper part of the Hertszprung–Russell (HR) diagram can lead to dense outflows. High massloss rates can veil the hydrostatic layers and result in the appearance of emission lines in the spectra of luminous stars (e.g. Hamann et al. 2006; Gräfener & Hamann 2008; Bestenlehner et al. 2014). In this case part of the radiation is scattered back from the radially expanding atmosphere, giving feedback in the form of a backwarming of the outer hydrostatic layers. This feedback is often described as proportional to the frequency independent mean optical depth τ, which is the average number of photon meanfree paths along the line of sight (Mihalas 1978). For optically thin winds the photosphere arises in the subsonic part of the outflow. More massive stars have winds that can become optically thick, where the continuum – or part of it – is produced in the supersonic part of the outflow.
The highest steady massloss rates by stellar wind are found in the late stages of the evolution of massive stars, in what is known as the Wolf–Rayet phenomenon. Wolf–Rayet (WR) stars are very luminous stars close to the Eddington limit and characterized by highly supersonic dense winds which can shroud the hydrostatic layers from direct observation (Langer et al. 1988; Langer 1989; Crowther 2008). The mass loss is thought to be driven by radiation pressure, i.e. momentum transfer via absorption and scattering of photons in the partially optically thick wind. High momentum transfer efficiency is necessary to explain the high massloss rates observed in WR stars, implying that multiple scattering events and enhanced opacities due to lineblanketing are necessary to accelerate these stellar winds up to the observed terminal wind velocities of the order of thousands of km s^{−1} (Abbott & Conti 1987; Lamers & Cassinelli 1999; Owocki 2013; Bestenlehner et al. 2014).
Empirical determinations of the hydrostatic radii of Galactic WR stars in recent decades has revealed a disagreement between the spectral analysis based on wind models (Hamann et al. 2006; Sander et al. 2012) and stellar structure calculations of mostly helium burning stars (Langer 1989; Heger & Langer 1996; Petrovic et al. 2006; Gräfener et al. 2012; Georgy et al. 2012; Grassitelli et al. 2016a). The difficulties in reconciling this disagreement, where there is a difference of up to a factor 10 in the derived hydrostatic radii (known as the “WR radius problem”), has been ascribed to the location of the wavelengthdependent photosphere (Hillier 1987, 2015; Langer et al. 1988; Kato & Iben 1992; Heger & Langer 1996) which forms in the supersonic part of the wind. Consequently, the detailed dynamics and density structure in the expanding envelope and supersonic wind are uncertain, which relates also to the complex opacities and their interplay with the accelerating outflow (Schaerer 1996; Heger & Langer 1996) and to the presence of density inhomogeneities in the wind (Moffat et al. 1988; StLouis et al. 2009; Chené et al. 2011; Michaux et al. 2014).
There appears to be a lack of radiative force to drive the high massloss rates observed in WR stars (e.g. Gräfener & Hamann 2005; Sander et al. 2015). This is true even considering the line opacity contribution of the absorption lines (especially UV) exposed to unattenuated radiative flux by the progressive Doppler shift of the stellar wind, which allows for the absorption of the unattenuated continuum, reducing thus the effect of line selfshadowing (Lucy & Solomon 1970; Castor et al. 1975; Abbott & Lucy 1985; Lucy & Abbott 1993; Owocki & Gayley 1997). Consequently, an approximate relation is commonly adopted to describe thedynamics of the winds (i.e. the velocity and density profiles), namely a betavelocity law (Castor et al. 1975; Langer 1989; Hamann et al. 2006; Hillier 2015; Sander et al. 2017). This approximation contributes to the uncertainty in the determination of the stellar radius and the detailed temperature stratification within the nonLTE supersonic outflow.
Hydrostatic models of helium stars above ~10 M_{⊙} with Galactic metallicity show the formation of a lowdensity inflated envelope having a density inversion below their plane parallel grey atmosphere (Petrovic et al. 2006; Gräfener et al. 2012; Grassitelli et al. 2016a; McClelland & Eldridge 2016). Inflation and density inversion arise as these models approach a local Eddington factor of unity in the proximity of the iron opacity bump at a temperature of ≈ 200 kK. This hydrostatic limit leads to the flattening of the temperature and density gradients in the outer layers, up to building a positive gas pressure gradient to counterbalance the strong radiative force (Joss et al. 1973; Langer 1997; Owocki 2015; Sanyal et al. 2015).
Heger & Langer (1996), Schaerer (1996), Petrovic et al. (2006), and more recently Ro & Matzner (2016) and Nakauchi & Saio (2018) have investigated the envelope configuration of massive stars while considering the effects of mass loss by stellar winds on the structure of the outer layers. The need for a strong radiative acceleration to launch these dense winds points towards the role of the iron opacity bump as a source of momentum for the flow already from the subsonic quasihydrostatic part of the star. Petrovic et al. (2006) showed that sufficiently high massloss rates can lead to configurations in which the inflation of the envelope is inhibited, and where the flow is accelerated by the radiative forces in the proximity of the iron opacity bump.
In this paper we investigate the interplay between an optically thick wind and the subsonic structure of massive helium stars, focusing especially on the conditions at the sonic point. In Sect. 2 we explore some theoretical considerations about optically thick winds and the importance of the sonic point. In Sect. 3 we apply our considerations to computing hydrodynamic stellar structure models, in Sect. 4 we show selected results and applications, in Sect. 5 we develop our results to constrain predicted massloss rates by introducing new useful tools, in Sect. 6 these results are discussed, and in Sect. 7 we draw our conclusions.
2 Optically thick winds
The requirement of a smooth, radiationpressuredriven, steadystate, transonic flow sets specific conditions that the physical quantities and their gradients must satisfy at the sonic point for optically thick winds. Based on the reasonable assumption of a radiationdriven wind (analogous to the O stars, Vink 2015), Nugis & Lamers (2002) estimated that the sonic point of hot WR stars is located at high optical depth (i.e. τ ≈3–30) and that it has to occur within the temperature ranges where the Rosseland opacity κ increases as a function of radius, i.e. where the flow accelerates as a result of the increase in the local opacity (d κ∕d r > 0), because energy and momentum input are necessary at the sonic point to accelerate the outflow.
From the momentum equation, the velocity gradient in Eulerian coordinates describing the hydrodynamics of a steadystate stellar wind can be written as (Bjorkman 1995; Lamers & Cassinelli 1999; Nugis & Lamers 2002; Hubeny & Mihalas 2014; Sander et al. 2017) $$\frac{1}{v}\frac{\text{d}v}{\text{d}r}=\left(g{g}_{\text{rad}}2\frac{{c}_{\text{s}}^{2}}{r}+\frac{\text{d}{c}_{\text{s}}^{2}}{\text{d}r}\right)/({v}^{2}{c}_{\text{s}}^{2}),$$(1)
where g and g_{rad} are the gravitational and radiative accelerations, c_{s} the local isothermal sound speed, v the flow velocity, and r the radial coordinate. Equation (1), known as the “Bondi equation” (Bondi 1952), implies that when v = c_{s}, d v∕d r diverges unless the numerator is null as well. The requirement of a null numerator at the sonic point implies that $${g}_{\text{rad}}=g2\frac{{c}_{\text{s}}^{2}}{r}+\frac{\text{d}{c}_{\text{s}}^{2}}{\text{d}r},$$(2)
where g = Gm∕r^{2} with m being the mass coordinate. In this equation the radiative acceleration can be expressed as $${g}_{\text{rad}}=\frac{\kappa L}{4\pi {r}^{2}c},$$(3)
with κ the Rosseland mean opacity from Thompson scattering and from boundbound and boundfree absorption, which is assumed in this context to be consistent with the fluxweighted mean opacity, and with the constants G and c holding their usual meaning. Assuming that the opacities are not affected by the velocity and the velocity gradient, Eq. (1) shows the presence of a critical point for the momentum equation at the sonic point (see Appendix C for further discussion).
Equation (2) can be written in terms of the Eddington factor Γ = g_{rad}∕g, and becomes $$\Gamma =1\frac{2\frac{{c}_{\text{s}}^{2}}{r}\frac{\text{d}{c}_{\text{s}}^{2}}{\text{d}r}}{g},$$(4)
which showsthat the condition $\frac{\text{d}{c}_{\text{s}}^{2}}{\text{d}r}2\frac{{c}_{\text{s}}^{2}}{r}\ll g$ leads to Γ ≈ 1 (or equivalently to a local effective gravity g_{eff} ≈ 0) at the sonic point. This has a crucial implication: for $\frac{\text{d}{c}_{\text{s}}^{2}}{\text{d}r}2\frac{{c}_{\text{s}}^{2}}{r}\ll g$, a radiationdriven, smooth, steadystate transonic outflow from a star needs to have its sonic point located at Γ ≈ 1. A transonic radiationpressuredriven stellar wind has to selfadjust its structure to meet this requirement.
Temperature stratification at the sonic point
Formal studies of the dynamics of stellar winds and their stability in spherical symmetry were first conducted by Parker (1958, 1966). However, from a mathematical and conceptual point of view, solving the physical problem of the dynamics of a gas outflow is analogous to studying spherically symmetric mass accretion (Bondi 1952; Tamazawa et al. 1975; Moncrief 1980; Thorne et al. 1981; Bjorkman 1995; Visser 1998). The topology of the steadystate solutions for the stellar wind equation shows two characteristic solutions: one that starts subsonic and reaches supersonic finite velocity (known as the wind solution) and one that starts supersonic and becomes subsonic for increasing distance (the accretion solution). The problem of stellar wind and accretion are conceptually analogous and, mutatis mutandis, the results obtained in one context can be applied to the other.
Accretion onto a black hole can also be divided into optically thin and optically thick accretion depending on the location of the sonic point radius with respect to the photosphere (Tamazawa et al. 1975; Thorne et al. 1981; Meier 1982; Flammang 1984; Nobili et al. 1991). In the optically thin case the hydrodynamics of the flow and the radiative transfer are effectively decoupled (Mihalas 1978; Thorne 1981; Flammang 1982). The nonlocal coupling of the radiation field to the structure of the atmosphere as a whole implies only a weak dependence of the energy transfer problem on the local thermodynamical conditions (Mihalas 1978). However, as the photon mean free path λ, defined via $$\lambda ={(\rho \kappa )}^{1}.$$(5)
becomes smaller than the length scale on which the macroscopic quantities change, the system enters into the photon diffusion regime (Mihalas 1978; Thorne et al. 1981; Kippenhahn & Weigert 1990). The radiative transport equation in the asymptotic limiting case of large optical depth and small mean free paths yields the diffusive approximation, which can be expressed as (Mihalas 1978; Kippenhahn & Weigert 1990) $$\frac{\text{d}{P}_{\text{rad}}}{\text{d}\tau}=\frac{F}{c},$$(6)
with P_{rad} the radiation pressure, F the flux, and dτ = −κρ dr.
Similarly, the stellar wind problem can be considered in the optically thin and optically thick regime. In Appendix B we investigate in more detail the behaviour of the radiation field in the proximity of the sonic point at the conditions typically met in stellar models of WNE stars. Backed by these considerations, we assume the conditions at the sonic point of WNE stars to closely match local thermodynamical equilibrium (LTE).
The LTE assumption implies that the temperature gradient can be defined by local quantities only, with the energy flux set by the need to transport the stellar luminosity outwards. At the high optical depth assumed for the sonic point, in fact, the temperature gradient becomes independent of the detailed global structure of the atmosphere, losing its explicit dependence on the optical depth characteristic of the optically thin situation (see Appendix B).
For a given massloss rate, the sonic point conditions of an optically thick wind (density and temperature) are uniquely defined, andthus the subsonic structure becomes independent from the detailed conditions above it, which only define the velocity profile and the terminal wind velocity of the outflow. Setting the massloss rate, in fact, sets the temperature at the sonic point. This can be seen combining the condition $$v=\sqrt{\frac{{k}_{\text{B}}T}{\mu {m}_{\text{H}}}}={c}_{\text{s}},$$(7)
where μ is the mean molecular weight, m_{H} the mass of a proton, and k_{B} the Boltzman constant, with the definition of steadystate massloss rate $$\dot{M}=4\pi {r}^{2}\rho v,$$(8)
which leads to $${T}_{\text{S}}=\frac{\mu {m}_{\text{H}}}{{k}_{\text{B}}}{\left(\frac{\dot{M}}{4\pi {r}^{2}\rho}\right)}^{2}.$$(9)
Therefore, having the massloss rate imposing the blanketed sonic point temperature while coupled to the subsonic stellar structure, hydrodynamic models can be computed independently from wind calculations without the need of a global computation (as we show in Sect. 4).
Envelope inflation
Both hydrostatic and hydrodynamic stellar models with high luminositytomass ratios close to the Eddington limit show diluted and extended envelopes, often characterized by very inefficient convective energy transport and the presence of a density inversion (Langer 1997; Ishii et al. 1999; Gräfener et al. 2012; Köhler et al. 2015; Sanyal et al. 2015, 2017; Grassitelli et al. 2015; Owocki 2015). These stellar models are said to be inflated, having a subphotospheric envelope configuration with large pressure and density scale heights and the local Eddington factor close to unity. This happens in the temperature range where an opacity bump arises, most notably the iron opacity bump.
Assuming that the temperature stratification is given only by radiative transport in the diffusive approximation, thus neglecting the contribution from convective energy transport (appropriate for the envelopes of massive helium star models, Grassitelli et al. 2016a), for chemically homogeneous layers with constant luminosity, using Eq. (6) the temperature profile can be written as $$\frac{\text{d}T}{\text{d}r}=\frac{\rho T}{4{P}_{\text{rad}}}\frac{\kappa L}{4\pi {r}^{2}c},$$(10)
from which $$\frac{\text{d}{c}_{\text{s}}^{2}}{\text{d}r}={c}_{\text{s}}^{2}\frac{\text{d}\text{\hspace{0.17em}}\mathrm{ln}(T)}{\text{d}r}=\frac{\beta}{4(1\beta )}\frac{\kappa L}{4\pi {r}^{2}c},$$(11)
with β the ratio of gas to total pressure. Combining then Eq. (10) with Eq. (1) and with the differential form of the steadystate continuity equation $$\text{d}\text{\hspace{0.17em}}\mathrm{ln}(v)+\text{d}\text{\hspace{0.17em}}\mathrm{ln}(\rho )+2\text{d}\text{\hspace{0.17em}}\mathrm{ln}(r)=0,$$(12)
the slope of the density profile follows $$\frac{\text{d}\text{\hspace{0.17em}}\mathrm{ln}(\rho )}{\text{d}r}=\left(g\frac{\kappa L}{4\pi {r}^{2}c}\left(1+\frac{{P}_{\text{gas}}}{4{P}_{\text{rad}}}\right)+\frac{2{v}^{2}}{r}\right)/({v}^{2}{c}_{\text{s}}^{2}).$$(13)
Here we candistinguish two regimes: the subsonic and the supersonic. In the subsonic, subEddington region, an increase in opacity or generally any outwarddirected force (e.g. the centrifugal force) tends to reduce the slope of the density profile, thus increasing the density scale height and giving rise to the envelope inflation encountered in the helium star models of e.g. Petrovic et al. (2006), Gräfener et al. (2012), and Grassitelli et al. (2016a). Consequently, although counterintuitive, from Eq. (13) an increase in the radiative force in the subsonic region implies a less steep increase in flow velocity, while making the density profile more and more flat; the layers aim to preserve hydrostatic equilibrium counterbalancing the local gravitational force, and the radiative acceleration effectively acts as a reduction in the local gravity. As the gas pressure gradient needed to counterbalance this reduced effective gravity is smaller, the increased opacity can lead, together with the acceleration of the flow, to an almost flat density profile when the local Γ approaches unity (Sanyal et al. 2017).
However, reaching and exceeding the local Eddington limit does not necessarily guarantee that the outflow reaches supersonic velocities. In case of a (still) subsonic steadystate flow, if the numerator on the righthand side of Eq. (13) becomes negative, a density inversion has to arise associated with a decrease in velocity, whereas a further increase in radiative force leads to a steep decrease in the velocity and an increase in the gas to total pressure. This hydrodynamic envelope inflation solution is similar to the “breeze” solution for winds which do not asymptotically reach supersonic velocities (Lamers & Cassinelli 1999; Hubeny & Mihalas 2014).
A clear definition, and therefore a criterion, of envelope inflation was not available prior to publication of this work (Sanyal et al. 2017). In Appendix A we discuss more extensively the appearance of inflated envelopes. We do so by performing hydrostatic stellar structure calculations of massive helium star models using different opacities and thus comparing how the outer layers react to radiative forces of different intensity. Noninflated hydrostatic massive helium star models show density and gas pressure scale height which rapidly decrease as they approach the surface. This is due to the gas pressure gradient being the leading force counterbalancing gravity and preserving hydrostatic equilibrium. However, in the presence of high Eddington factors in stars with high L∕M ratios, an inflection point in density and gas pressure can appear, and thus density and gas pressure scale heights start to increase.
We show that an inflection point can appear as a consequence of a rapid increase in opacity (or equivalently in the Eddington factor) when the radiation pressure gradient becomes the primary agent counterbalancing gravity. Analytically, an inflection point in the gas pressure profile is found when (see Eq. (A.8)) $$\frac{\text{d}\text{\hspace{0.17em}}\mathrm{ln}(\kappa )}{\text{d}\text{\hspace{0.17em}}\mathrm{ln}(T)}\lesssim 1\frac{1}{\Gamma},$$(14)
showing both the need for an Eddington factor close to unity and for a steep increase in the opacity with temperature. In view of this discussion and of Appendix A, we can thus understand envelope inflation as the appearance of large gas pressure scale heights in response to a steep increase in the radiative opacity in radiationpressuredominated layers close to the Eddington limit. The low gas pressure gradient prevents the density from decreasing steeply for extended regions of space, and implies that the surface is found at significantly larger radii compared to noninflated stellar models.
3 Model assumptions
Computing the stellar structure is a task that requires solving a hyperbolic set of differential equations simultaneously with welldefined boundary conditions (Kippenhahn & Weigert 1990). We adopt the Bonn evolutionary code (BEC), a Lagrangian hydrodynamic onedimensional stellar evolution code (Langer et al. 1988, 1994; Heger & Langer 1996; Heger 1998; Heger et al. 2000; Petrovic et al. 2006; Yoon et al. 2006; Kozyreva et al. 2014). It solves the set of coupled nonlinear partial differential stellar structure equations in the form (Kippenhahn & Weigert 1990) $$\begin{array}{l}{\left(\frac{\partial m}{\partial r}\right)}_{t}=4\pi {r}^{2}\rho ,\\ \%{\left(\frac{\partial r}{\partial t}\right)}_{m}=v,\\ {\left(\frac{a}{4\pi {r}^{2}}\right)}_{t}=\frac{Gm}{4\pi {r}^{4}}+\frac{\partial P}{\partial m},\\ \%,{\left(\frac{\partial T}{\partial m}\right)}_{t}=\frac{Gm}{4\pi {r}^{4}}\frac{T}{P}\nabla {\left(1+\frac{{r}^{2}}{Gm}\frac{\partial v}{\partial t}\right)}_{m},\\ {\left(\frac{\partial L}{\partial m}\right)}_{t}={\u03f5}_{n}{\u03f5}_{g}{\u03f5}_{\nu},\end{array}$$
where m and t, the mass coordinate and time, are the two independent variables; ρ is the density; v is the velocity; a is the acceleration; r is the radial coordinate; T is the temperature; L is the luminosity; $\nabla :=\frac{\text{d}\text{\hspace{0.17em}}\mathrm{log}(T)}{\text{d}\text{\hspace{0.17em}}\mathrm{log}(P)}$ is the temperature gradient; P is the total pressure given by the sum of gas and radiation pressure; G is the gravitational constant; and ϵ is the energy production/loss per unit mass and unit time related to nuclear processes (subscript n), to gravitational contraction/expansion (subscript g), and neutrino loss (subscript ν). These equations express: the conservation of mass (Eq. (15)), the definition of velocity (Eq. (16)), the conservation of momentum (Eq. (17)), the energy transport (Eq. (18)), and the energy conservation (Eq. (19)). These equations together with the network of nuclear reaction rates, the set of equations of the mixing length theory (BöhmVitense 1958), and the OPAL opacity tables (Iglesias & Rogers 1996), define the structure and evolution of a stellar model. In addition, mass loss by stellar wind can also be applied.
The set of nonlinear coupled partial differential equations above has five dependent variables (namely ρ, L, v, T, and r) and requires a number of boundary conditions equal to the degrees of freedom of the system, i.e. five. The central boundary conditions are trivially set by the requirement of zero r, L, and v in the centre of the stellar models, while the outer boundary conditions are usually determined by the assumption of plane parallel grey atmosphere (Langer 1989; Heger 1998).
The boundary conditions set constraints on the family of possible solutions (if any) of the set of partial differential equations describing a star. We modify BEC adopting as surface boundary conditions $$\dot{M}=4\pi {r}^{2}\rho v,$$(20)
and $$v=\sqrt{\frac{{k}_{\text{B}}T}{\mu {m}_{\text{H}}}}={c}_{\text{s}},$$(21)
i.e. we impose the boundary of the stellar model at the sonic point while preserving continuity. The massloss (or massaccretion) rate is a free parameter in our calculations, i.e. it can be chosen either manually or according to mass loss by stellar wind prescriptions (such as those in e.g. Nugis & Lamers 2000; Vink et al. 2001). Mass loss is treated adopting a pseudoLagrangian scheme in the outer part of the stellar model, i.e. the independent variable becomes q = m∕M_{tot} (the relative mass coordinate) and consequently the Lagrangian operator time derivative for a steady state becomes $\frac{\text{d}}{\text{d}t}=\frac{q\dot{M}}{{M}_{\text{tot}}}\frac{\partial}{\partial q}$ (Neo et al. 1977; Heger 1998).
We adopt BEC to compute the subsonic structure of massive helium stars. Models are computed with a chemical composition as in Heger et al. (2000), with metallicity Z = 0.02 and helium fraction 0.98, and stationarity is achieved with large time steps (i.e. greater than 10^{3} seconds) while inhibiting chemical evolution. The velocity profile is computed selfconsistently with the radiative acceleration derived via the velocityindependent Rosseland opacity and the imposed massloss rate. Convection is included in the calculations, but as shown by Grassitelli et al. (2016a), the convective flux in the outer layers is several orders of magnitude smaller than the radiative flux. We compute models for chemically homogeneous massive helium zero age main sequence stars of mass 10, 15, and 20 M_{⊙} at solar metallicity and under various applied massloss rates by stellar wind. Rotation and magnetic field are not included in the calculations;however, they are not expected to affect the general theoretical considerations in Sect. 2. We also neglect the structural effects from the turbulent pressure on the stellar models. For comparison purposes, we also compute models with classical plane parallel grey atmosphere boundary conditions (cf. Heger 1998; Petrovic et al. 2006).
4 Results
Figure 1 shows the velocity profiles of the outflow in the outer subsonic part of a set of 15 M_{⊙} stellar models computed with sonic point boundary conditions and with different adopted massloss rates. For comparison the velocity profile of a plane parallel grey atmosphere BEC model is also shown.
All the models presented in Fig. 1 show a rapid increase in the radial velocity starting from R ≈ 1.13 R_{⊙}. This increase is associated with the increase in radiative force in correspondence with the iron opacity bump (Febump) in the temperature range 5.2 ≲log(T) ≲ 5.5. The more compact, black models in Fig. 1 show that the velocity profile is steeper for the higher massloss rate applied, with a flow velocity monotonically increasing until it reaches the sound speed in the temperature range of the hot rising part of the Febump. The model with the highest massloss rate (5 × 10^{−5} M_{⊙} yr^{−1}) achieves a transonic flow from the smallest radius, while the other stellar models show larger sonic point radii as the adopted massloss rates decrease (see Table 1 for the parameters of the models). All the monotonic solutions are thus solutions where the outflow is accelerated by the iron opacity (see also Fig. 2). Compared to the electron scattering opacity, thisincrease in opacity provides a stronger outwarddirected radiative force, leading to higher local Eddington factors and flows that can reach supersonic velocities already at 200 kK. The higher the massloss rate, the steeper the velocity profile, the sooner (in terms of radius) the flow meets the sound speed, and the hotter and more compact is the resulting model.
However, solutions with nonmonotonic velocity profiles are also present. While the group of more compact monotonic solutions has massloss rates of Ṁ ≥ 5 × 10^{−6} M_{⊙} yr^{−1} (see Table 1), the pink model with a massloss rate of Ṁ = 4 × 10^{−6} M_{⊙} yr^{−1} shows an outflow that does not become supersonic in the proximity of the iron opacity bump. In this case the flow slows down to velocities of the order of a few km/s after being accelerated by the increased opacity of the Febump. It then reaccelerates and finds the sonic point in the hot rising part of the helium opacity bump (Hebump) when log (T) ≈ 4.6−5. This kind of extended solution implies the formation of a density inversion following the decrease in flow velocity. This solution approaches and exceeds Γ =1 in the envelope (in correspondence with the positive density gradient), but shows an outflow that does not exceed the local sound speed before the opacity peak of the Febump (log(T) ≈ 5.2). For the 15 M_{⊙} model, however, even in these more extended solutions, the sonic point radius does not increase by more than 10% compared to the more compact solutions.
For comparison, a stellar model with plane parallel grey atmosphere boundary conditions and a massloss rate of 10^{−5} M_{⊙} yr^{−1} is plotted in Fig. 1, which shows an outflow accelerating to supersonic velocities near the Febump, but then decelerating to give rise to a density inversion below the peak temperature of the iron opacity bump. This model also shows that the subsonic structure of the model computed with plane parallel grey atmosphere boundary conditions and the one with sonic point boundary conditions do not differ. This is the case even if the first has layers above the sonic point while the second has no information concerning the conditions above the sonic point other than the prescribed massloss rate.
Figure 2 shows the density inversion in the plane parallel and Hebump solutions, with a peak density one order of magnitude higher than the underlying layers. The inflection point in the density profile (i.e. d ^{2} ρ∕dr^{2} = 0) is also visible at R ≈ 1.13 R_{⊙} , indicating the location where inflation begins (see Sect. 2). This radius is defined as the core radius. For the compact, Febump solutions, the larger the radius, the lower the density of the sonic point, and the larger the density scale height. Figure 2 also shows the steep increase in opacity due to the recombination of iron encountered at the base of envelope. This Febump increases the opacity with respect to the electron scattering opacity by almost a factor of 3, leading to the appearance of the inflection point and providing the momentum to accelerate the flows up to the sonic point from temperatures as high as a few hundred thousand kelvins.
As a consequence of the need for a dκ∕dr > 0 (Nugis & Lamers 2002), for the 15 M_{⊙} helium star models we find two types of configurations depending on the applied massloss rate:

hot, Febump stellar models having corresponding sonic temperatures which are higher than log(T_{S}) ≳ 5.2 and radii similar to the coresize;

cooler, slightly more extended Hebump models with sonic points in the range log(T_{S}) ≈ 4.6−4.8.
Similarly, in Fig. 3 we show how a 10 M_{⊙} helium star model readjusts to the different massloss rates. As in Fig. 1, higher massloss rates lead to more compact solutions, with steeper velocity profiles. Differently from the 15 M_{⊙} model, the more compact and less luminous 10 M_{⊙} helium star models find the sonic point in the helium opacity bump already starting from Ṁ = 1 × 10^{−5} M_{⊙} yr^{−1} (see Table 1). This is due to the lower luminosities and higher densities in the envelopes of these models.
The opposite is true for the 20 M_{⊙} model in Fig. 4, for which the solutions reach the sonic point in the hot part of the Febump at ≈ 1.4R_{⊙} for all the applied massloss rates. In this case the high luminosities (and low densities in the envelope) of the models allow a transonic flow, even for relatively low massloss rates. Figure 4 also shows how the model with a plane parallel grey atmosphere behaves in the supersonic part of the outflow (for Ṁ = 8 × 10^{−6} M_{⊙} yr^{−1}), which rapidly reaches flow velocities of more than 100 km s^{−1}, to then slow down again following the density inversion due to the decrease in opacity between the Fe and Hebump.
Sonic point conditions
In Table 1, the physical quantities computed at the sonic point are reported for our stellar models. The mean free paths estimated at the sonic point are small, of the order of 0.01–1% of the sonic radius and an order of magnitude smaller than the pressure scale height at the sonic point. They tend to increase for the lower mass loss by stellar wind applied and for the Hebump solutions, mostly due to the lower sonic point densities. From Appendix B and the analysis concerning the radiation field at the sonic point, mean free paths smaller than the scale at which the thermodynamic quantities change suggest that local conditions are not far from LTE.
For the diffusive approximation to be valid, the mean free path should also be considered together with the optical depth of the sonic point. We estimate the optical depth of the wind via the integral $${\tau}_{\text{S}}={\displaystyle {\int}_{{R}_{\text{S}}}^{\infty}\text{d}}r\text{\hspace{0.17em}}\text{\hspace{0.17em}}\kappa (\rho ,T)\text{\hspace{0.17em}}\rho ,$$(22)
from the sonic point radius R_{S} to infinity. In the supersonic region ρ is computed via the continuity equation (Eq. (20)) while adopting a betavelocity law (Langer 1989, Eq. (11)) connected continuously on top of the stellar model, having the exponent unity and a typical velocity at infinity of 1600 km/s. The opacity κ is instead computed via the OPAL opacity tables throughout the wind. To evaluate the opacity, the OPAL tables need as input not only the local density, but also the temperature, approximated in this case via the differential form of the T − τ relation from Nugis & Lamers (2002) with the sonic point temperature as reference temperature. These values, which should be considered only a rough estimate of the actual optical depth, are shown in Table 1 and indicate that the optical depth of the sonic point is of the order of 10 in case of massloss rates greater than log (Ṁ) ≳−5. Instead, for log(Ṁ) ≲−5 the optical depth becomes close to unity and the mean free path at the sonic point becomes quite large, indicating that the assumption of diffusive energy transport and LTE might not be valid for these massloss rates.
Another key difference between the Hebump and the Febump models can be seen considering the optical depth parameter (Castor et al. 1975; Mihalas 1978; Lamers & Cassinelli 1999) $$t=\frac{{v}_{\text{th}}}{\lambda}/\frac{\text{d}v}{\text{d}r},$$(23)
with v_{th} being the thermal velocity of the particles (Lamers & Cassinelli 1999). The optical depth parameter tends to be of theorder of 10 in the hot Febump solutions, while it drops to less than 1 in the Hebump case. Considering that the optical depth parameter indicates how well the velocity independent Rosseland mean opacity can reproduce the fluxweighted opacity also in the presence of a velocity gradient, low values of t point to the inadequacy of the use of the OPAL opacity tables in this context. For low t, the radiative acceleration due to the progressive linedeshadowing associated with the velocity gradient could instead already be important for launching these winds.
A key aspect that emerges from Table 1 is that the sonic point coincides to a very good approximation with an Eddington factor of unity. In general the Eddington limit can be rewritten in the form of an Eddington opacity κ_{EDD}, $${\kappa}_{\text{EDD}}:=4\pi cGM/L,$$(24)
only proportional to the masstoluminosity ratio; this is the opacity necessary for a star in hydrostatic equilibrium to reach the Eddington limit.
Figure 5 shows the Rosseland opacity from the OPAL tables (Iglesias & Rogers 1996) as a function of temperatureand density, together with the structure of the outer layers of the stellar models shown in Fig. 1. Associating an Eddington opacity with the opacity contours in Fig. 1, we can define the contour of Γ =1 relative to a given L∕M. This is thus the contour where we expect to find the sonic points of our stellar models. In Fig. 5 the Eddington opacities relative to our 10, 15, and 20 M_{⊙} models are plotted. They follow the increase in opacity around log(T) ≈ 5.2, i.e. the Febump, and at log(T) ≈ 4.6, i.e. the HeII opacity bump. The more compact stellar models in Figs. 1, 3, and 4 have their sonic points located almost exactly on their respective opacity contours and in the hot part of the iron opacity bump, i.e. 5.2 < log(T) < 5.5. The two more extended models showing a density inversion, i.e. the Hebump sonic point model and the plane parallel grey atmosphere model, have their outer layers follow closely the contour of Γ ≈ 1 for log (T) ≳ 5.2 (see Sect. 2), to then cross this limit within their density inversions in the proximity of the peak temperature of the Febump. The plane parallel model finds its surface back into the subEddington region, while the Hebump models continues to cooler temperatures to then find its sonic point close to its Eddington opacity contour at log (T) ≈ 4.7. The grey shaded regions in Fig. 5 mark instead combinations of parameters for which we do not expect to find the sonic point of radiationdriven stellar winds, due to the decrease in opacity as a function of temperature.
Fig. 1 Velocity as a function of radius for a set of 15 M_{⊙} stellar models with different massloss rates and boundary conditions. The black continuous lines (from right to left) indicate models with Ṁ equal to 5, 8,10, 20, 30, 50 × 10^{−6} M_{⊙} yr^{−1} and the pink dashed line indicates the model with Ṁ = 4 × 10^{−6} M_{⊙} yr^{−1}, all having sonic point boundary conditions. The orange dotdashed lines indicates the model with Ṁ = 1 × 10^{−5} M_{⊙} yr^{−1} computed with plane parallel grey atmosphere boundary conditions. The green dashed line indicates the local isothermal sound speedfrom the pink dashed model. 
Fig. 2 Density profile (top, in units of g cm^{−3}) and opacity profile (bottom, in units of cm^{2} g^{−1}) as a functionof radius for the same 15 M_{⊙} stellar models with different massloss rates and boundary conditions as in Fig. 1. The colours are as in Fig. 1, and the red circles indicate the location of the sonic point and the green dashed line the Eddington opacity (as defined in Eq. (24) in Sect. 4) of the pink dashed model. 
Fig. 3 Velocity profile as a function of radius for a set of 10 M_{⊙} stellar models with different massloss rates (as in Fig. 1). The black continuous lines indicate models with Ṁ ≥ 2 × 10^{−5} M_{⊙} yr^{−1} (2, 3, 5 × 10^{−5}), while the pink dashed lines indicate models with Ṁ ≤ 1 × 10^{−5} M_{⊙} yr^{−1} (5, 8, 10 × 10^{−6}). The green dashed line and the orange dotdashed line indicate the local isothermal sound speed of the model with Ṁ = 5 × 10^{−6} M_{⊙} yr^{−1} and the plane parallel grey atmosphere model with Ṁ = 8 × 10^{−6} M_{⊙} yr^{−1}, respectively. 
Fig. 4 Velocity profile as a function of radius for a set of 20 M_{⊙} stellar models with different massloss rates (as in Fig. 1). The black continuous lines indicate models with Ṁ ≥ 5 × 10^{−6} M_{⊙} yr^{−1} (5, 10, 20, 30, 50 × 10^{−6}), while the orange dotdashed lines indicate models with Ṁ = 8 × 10^{−6} M_{⊙} yr^{−1} computed with plane parallel grey atmosphere. The green dashed line indicates the local isothermal sound speed. 
Fig. 5 Contour plot showing the opacity (in cm^{2} g^{−2} ) from OPAL tables (colourcoded, see scale at right) as a function of temperature (in K) and of density (in g cm^{−3} ) for a chemical composition of [X, Y, Z] = [0, 0.98, 0.02], with X hydrogen mass fraction, Y helium mass fraction, and Z metallicity. The black lines correspond to the outer layers of the hot, more compact 15 M_{⊙} stellar models with Ṁ ≥ 8 × 10^{−6} M_{⊙} yr^{−1} from Fig. 1 (see also Table 1), while the pink line indicates the cooler model with Ṁ = 4 × 10^{−6} M_{⊙} yr^{−1}, and the orange dotdashed line the plane parallel grey atmosphere model with Ṁ = ×10^{−5} M_{⊙} yr^{−1}. The yellow circles indicate the locations of the sonic points of this set of models, while the white lines indicate the isocontours of the Eddington opacity for the 10 M_{⊙} (dotted), 15 M_{⊙} (continuous), and 20 M_{⊙} (dashed) stellar models. Grey areas indicate regions in the log(ρ) − log(T) diagramwhere the sonic point of a radiationdriven wind cannot be located (d κ∕d r ≤ 0). 
5 Massloss rate constraints
We can make use of Eq. (9), (21), and (24) to obtain the stellar wind massloss rates of WR stars as a function of their sonic point temperatures and luminositytomass ratios. This can be done via the opacity tables, assuming as before that the value of Γ is exactly one at the sonic point (see Fig. 5 and Table 1). From this assumption it follows that the sonic point opacity is equal to the Eddington opacity, which is only proportional to the ratio M∕L (Eq. (24)). Manipulating the OPAL tables we can unequivocally associate a sonic point density ρ_{S} with a sonicpoint temperature and a given Eddington opacity (as in Fig. 5). Combining the sonic point density and the isothermal sound speed (which is a function of the temperature only) the mass flux at the sonic point can be obtained. Once the mass flux is defined, a radius still needs to be specified to obtain a massloss rate (Eq. (20)). In Fig. 6 we assume a typical sonic point radius of 1 R_{⊙} for all luminositytomass ratios, obtaining a sonic HR diagram^{1} from which the massloss rate can be derived as a function of the sonic point temperature and the luminositytomass ratio.
Two assumptions are used to construct this diagram: the validity of the velocity independent OPAL opacities and the adoption of a specific radius to derive massloss rates from massfluxes. However, as previously discussed, the sonic point Rosseland opacities are not expected to differ significantly due to the effects of a velocity gradient, at least in the case of the Febumpdriven flows. Concerning the adopted radius, as shown in Sect. 4, solutions with the sonic point in the Febump are compact, due to the lack of a strongly inflated envelope. The core and sonic point radii in the Febump tend to be very similar and are not expected to differ by more than a factor of 2 from 1 R_{⊙} , at least not in the mass range of most Hfree WR stars, i.e. 10–20 M_{⊙} ^{2}.
In Fig. 6 it can be seen how the massloss rate contours follow the opacity profile of the iron opacity bump, as in Fig. 5. Above log(T_{S}) ≈ 5.2, the higher the sonic point temperature, the higher the expected massloss rate at fixed L∕M (see Fig. 1). Fixing instead the sonic point temperature, the higher the L∕M, the lower the massloss rate, which indicates that as the L∕M ratio increases, an Eddington factor of unity can be more easily reached, and thus flows can be accelerated to supersonic velocities already for lower massloss rates. Equivalently, given Ṁ (thus a contour in the diagram), the flow can be accelerated until it reaches the Eddington limit, and therefore transonic velocities, already from higher temperatures for higher L∕M ratios.
In contrast, below log(T_{S}) ≈ 5.2 the predictions tend to be much more uncertain, mostly because sonic point temperatures lower than the iron bump peak imply that the stellar wind is not accelerated to supersonic velocities by the iron opacity bump, but rather the opacity bump leads to a lowdensity inflated envelope configuration with a density inversion. In this configuration the assumption of a typical radius starts to break down together with the assumption of LTE in the subsonic part of the outflow (see Table 1).
Similarly, we can derive the L∕M ratio as a function of the sonic point temperature and the massloss rate. This is done in Fig. 7, where (thanks again to Eq. (24)) from the sonic point temperature and log (Ṁ) we can obtain κ_{EDD} and thus the L∕M ratio. For log (T_{S}) ≳ 5.2, higher sonic point temperatures lead to higher massloss rates at fixed L∕M, or at fixed Ṁ a higher luminositytomass ratio requires a higher sonic point temperature. The same trend can also be seen in the temperature range of the raising helium opacity bump, i.e. 4.6 ≲ log(T_{S}) ≲ 5.0, but with L∕M ratios which are lower at constant Ṁ compared to the temperatures in the hot part of the iron opacity bump. This can be interpreted as stars needing to have higher L∕M in order to drive a certain Ṁ while having the sonic point located in the iron bump compared to stars having their winds driven by the helium opacity bump^{3} . Moreover, in the temperature range of the helium opacity bump, flows can be accelerated to supersonic velocities already in stars with significantly lower L∕M than the values necessary in the iron opacity bump.
Fig. 6 Contour plot showing the massloss rate (colourcoded in units of M_{⊙} yr^{−1}) as a function of sonic point temperature (in K) and of the luminositytomass ratio (in L_{⊙} / M_{⊙} ) of a star. Contours are derived based on the OPAL tables for a chemical composition [X,Y,Z]=[0,0.98,0.02], with X hydrogen mass fraction, Y helium mass fraction, and Z metallicity. 
Fig. 7 Contour plot showing the luminositytomass ratio (colourcoded in units of L_{⊙} ∕ M_{⊙} ) as a function of sonic point temperature (in K) and of the massloss rate (in M_{⊙} yr ^{−1}) of a star. Contours are derived based on the OPAL tables for a chemical composition [X,Y,Z]=[0,0.98,0.02], with X hydrogen mass fraction, Y helium mass fraction, and Z metallicity. 
Comparison to observations
Figure 6 shows that a lowmass, lowluminosity stellar model does not have an outflow that reaches transonic velocities within the iron opacity bump unless a certain minimum massloss rate is applied. Lower massloss rates lead instead to inflated envelope configurations such as the low massloss rate cases in Figs. 1 and 3.
A minimum massloss rate as a function of the luminositytomass ratio can therefore be derived, and indicates which massloss rate implies a radiationdriven supersonic flow starting from the iron opacity bump (minimum Ṁ _{Fe}). This minimum massloss rate marks the separation between compact and extended solutions in Figs. 1 and 3. It is shown in Fig. 8, where the derived minimum L∕M above which the flow becomes transonic in the hot part of the iron opacity bump is plotted as a function of the massloss rate. The minimum Ṁ _{Fe} decreases at higher luminosities, and it can be closely approximated by a parabola $$\mathrm{log}\left(\frac{L}{M}\right)=1.690.80\text{\hspace{0.17em}}\mathrm{log}\left(\dot{M}\right)0.06\text{\hspace{0.17em}}\mathrm{log}{\left(\dot{M}\right)}^{2},$$(25)
or equivalently in terms of luminosity $$\mathrm{log}\left(L\right)=0.441.35\text{\hspace{0.17em}}\mathrm{log}\left(\dot{M}\right)0.097\text{\hspace{0.17em}}\mathrm{log}{\left(\dot{M}\right)}^{2}.$$(26)
In Fig. 8 this minimum massloss rate as a function of L∕M is comparedto the sample of observed Galactic WNE analysed by Hamann et al. (2006). We collected L∕M and Ṁ exclusively for the Hfree WNE stars, as even a small mass fraction of hydrogen can significantly affect the structure of the envelopes (Schootemeijer & Langer 2017). The vast majority of the observed stars exceed the minimum massloss rate. They are therefore consistent with models having outflows that are radiation pressure driven to supersonic velocities by the iron bump.
Assuming the massluminosity relation from Langer (1989) and, as before, a typical radius of 1 R_{⊙} , we can construct another sonic HR diagram relating luminosity and sonic point temperature. This is done in Fig. 9, where regions are colourcoded according to the expected massloss rate. In this diagram we locate the observed WNE stars from Fig. 8. Their location clusters at sonic point temperatures around log (T_{S}) ≈ 5.3. Therefore, based solely on their observed luminosities and massloss rates, all the WNE stars of this sample can be understood as having their outflows accelerated to supersonic velocities by the radiation pressure consequenceof the high number of transitions associated with iron and irongroup elements at around log (T_{S}) ≈ 5.2. The sonic point temperature as a function of luminosity for the best lintear fit of the observed WNE stars in Fig. 9 follows the relation $$\mathrm{log}(L)=3.18\text{\hspace{0.17em}}\mathrm{log}({T}_{\text{S}})11.36,$$(27)
suggesting that, in general, WNE stars with higher luminosities have higher sonic point temperatures. A similar relation as in Eq. (27) can be derived also in terms of massloss rate and sonic point temperature, $$\mathrm{log}(\dot{M})=5.66\text{\hspace{0.17em}}\mathrm{log}({T}_{\text{S}})\mathrm{34.65.}$$(28)
This time Eq. (28) suggests that higher massloss rates lead to higher sonic point temperatures, or in other words that stars with high massloss rates have their sonic point located deeper inside their atmospheres (see also Table 1).
Fig. 8 Minimum luminositytomass ratio (in units of L_{⊙} / M_{⊙} , left vertical axis) and corresponding luminosity according to the mass–luminosity relation from Langer (1989, right vertical axis), as a function of massloss rate (in M_{⊙} yr ^{−1}) for an outflow to be driven to supersonic velocities by the radiative acceleration from the iron opacity bump. The green symbols indicate Hfree WNE stars from Hamann et al. (2006), with the different shapes indicating different spectral subclasses. 
Fig. 9 Sonic HR diagram relating the luminosity (in L_{⊙} ) to the sonic point temperature (in K). The background is colourcoded according to the massloss rates by stellar wind (legend on the right, in units of M_{⊙} yr ^{−1}) predicted via the use of the OPAL opacity tables and the proximity of the sonic point to the Eddington limit (see Fig. 6). The black symbols are observed Galactic WNE stars from Hamann et al. (2006), located in this diagram viatheir observed luminosity and massloss rates and with symbols indicating their spectral subclass, as in Fig. 8. The white dashed line is the best fitting linear relation for the observed WNE stars. 
6 Discussion
6.1 WR radius problem
For typical Galactic WNE stars the continuum originates in the dense outflow at a significant fraction of the terminal velocities, i.e. well beyond the hydrostatic domain (Crowther 2008). Therefore, the velocity law at the base of the wind is not well constrained by observations (Schmutz et al. 1992; Hillier 2015).
As shown in Sect. 4, the subsonic structure of helium star models suggests that the location of the sonic point of WNE stars is in the range 0.9–1.5 R_{⊙} , with supersonic flows originating from temperatures around 200 kK (see Fig. 9).
For the supersonic layers, Sander et al. (2017) has recently outlined how for hydrodynamically consistent models of hot, radiationdriven stellar winds the velocity profile is shaped by the different ionization levels of the elements present in the wind. This suggests that the adoption of a single betavelocity law for the supersonic wind, independently of the detailed temperature stratification, available opacities, and subsonic structure is probably an oversimplification (see also Gräfener & Hamann 2005). Hydrodynamic models for WR winds might bring the hydrostatic radii in the range predicted by our stellar structure calculation (see e.g. Fig. 2 in Sander et al. 2015). A nonconstant velocity gradient was found in the supersonic part of the flow of WR111, a WC5 star, for which a hydrodynamic model for the optically thick wind was available (Gräfener & Hamann 2005). We were able to continuously connect the subsonic structure of one of our sonic point boundary conditions models and the supersonic wind model at the sonic point (see Appendix D).
Still, the apparent homogeneity in the location of the observed WNE stars in the sonic HR diagram from Fig. 9 does not yetexplain the scattered distribution found in the classical HR diagram (cf. Fig. 8 in Hamann et al. 2006). From Langer (1989), at a given metallicity the properties of Hfree postmain sequence stars are expected to be mostly defined by one single parameter, the luminosity of the star. Instead WNE stars with comparable luminosities are found to have different spectroscopic subclasses. In this respect in Fig. 9 we note how the predicted massloss rate depends sensitively on the sonic point temperature. As such, small differences in the metallicity or in the subsonic force balance (e.g. due to rotation or magnetic fields), or the exact evolutionary phase and chemical composition of the core, might result in significantly different massloss rates, which consequently could lead to differences in the photospheric stellar parameters, which might reconcile observations with theoretical expectations.
A decelerating flow above the sonic point could possibly lead to configurations with multiple critical points and the appearance of a density inversion similar to that observed in inflated envelopes. The presence in massive main sequence stars of an increase in density shortly after the peak temperature of the iron opacity bump has recently been confirmed via 3D simulations by Jiang et al. (2015), implying that the density inversion is not necessarily erased, especially not by subsonic flows. The stagnation could then possibly imply the appearance of instabilities such as photon bubbles (Dessart & Owocki 2003; Blaes & Socrates 2003; Owocki 2015), linedeshadowing instabilities (Owocki & Rybicki 1984; Sundqvist & Owocki 2013), or streamlike lateral flows possibly connected with the corotating interaction regions (Mullan 1984; StLouis et al. 2009). Hydrodynamic and timedependent multidimensional simulations are necessary to assess the dynamical stability of such configurations and might shed light on what ultimately sets the massloss rate of these objects.
6.2 Optical depth of the sonic point
Integrating the diffusive temperature gradient (Eq. (10)) from the sonic point to infinity $${\int}_{\infty}^{{r}_{\text{S}}}\kappa}\rho \text{\hspace{0.17em}}\text{d}r={\displaystyle {\int}_{{T}_{\infty}}^{{T}_{\text{S}}}\frac{4ac}{3F}}\text{\hspace{0.17em}}\text{d}{T}^{4},$$(29)
with a the radiation density constant, and assuming that the sonic point temperature ${T}_{\text{S}}^{4}\gg {T}_{\text{eff}}^{4}>{T}_{\infty}^{4}$, valid in the case of optically thick winds, leads to a tautological relation between optical depth and sonic point temperature $${\tau}_{\text{S}}\approx \frac{ac}{3F}{T}_{\text{S}}^{4}.$$(30)
This, combined with Eq. (9), sets a relation which needs to be fulfilled at the sonic point for a given massloss rate, i.e. $${\left(\frac{3F}{ac}{\tau}_{\text{S}}\right)}^{1/4}\approx \frac{\mu {m}_{\text{H}}}{{k}_{\text{B}}}{\left(\frac{\dot{M}}{4\pi {r}^{2}\rho}\right)}^{2}.$$(31)
The massloss rate is a free parameter in our hydrodynamic stellar models. Equation (31) can then be used to constrain the only massloss rate which is consistent with the expected optical depth if combined with a prescribed velocity law for the wind (as in Gräfener et al. 2017), or assuming an optical depth for the sonic point.
6.3 Turbulence
Atmosphere calculations of hot stellar winds often include an extra turbulent broadening of the lines, associated with the turbulent motion in the atmosphere (e.g. Dessart & Owocki 2005; Sander et al. 2017). This turbulent velocity contributes to the equation of state of stellar matter and, via its gradients, to the structure of the outer layers^{4} . From the stellar structure calculations point of view, the structural effects of the inclusion of the turbulent pressure terms in the convective zones are marginal (Grassitelli et al. 2015). Moreover, subsurface convection in the subsonic layers might be inhibited in helium star models (Ro & Matzner 2016). Therefore, we neglect the turbulent terms from our analysis, expecting at most a difference of a fraction of the gas pressure in the equation of state (Grassitelli et al. 2015, 2016b). If present, this turbulent motion might alter our results, and introduce a systematic difference between our analysis and stellar wind models.
6.4 Previous works
Compact stellar structure models were also favoured by wind models from Ro & Matzner (2016). Considering the sonic point at the Febump as a boundary condition for their wind models, these authors show that while using velocity independent Rosseland opacities, winds fail to accelerate up to the terminal wind velocities. However, the dynamics in the supersonic part of the flow show the presence of more compact, strong wind solutions and extended, weak stagnating solutions. The bifurcation between these solutions that takes place depends on the ratio between the temperature scale height (computed in the diffusive limit within the supersonic part of the wind) and the local radius. In particular, their more compact wind solutions are expected to be more prominentlyaffected by lineforce amplification due to Doppler enhancement from flow velocities of the order 150–200 km s^{−1}, suggesting that the Rosseland mean opacities become rapidly inadequate above the sonic point. The question is whether such an increase in opacity is sufficient to accelerate the outflow monotonically up to the escape speed of the star, or whether stagnation and a complex velocity profile might appear.
After this manuscript was submitted, Nakauchi & Saio (2018) constructed hydrostatic Hestar models connected to a single betavelocity law, optically thick stellar wind models. They confirm the crucial role played by the iron opacity bump in launching the winds of WNE stars, in agreement with our fully hydrodynamic models. Less satisfactory is the comparison between theobserved photospheric temperature and the temperature values estimated by their wind models, further pointing towards the inadequacy of a prescribed velocity law with a single beta exponent in reproducing the initial acceleration of these stellar winds. Similarly, Gräfener et al. (2017) made use of the proximity to Γ ≈1 at the sonic point (see Sect. 2) while adopting a prescribed stellar wind dynamics to investigate WR massloss rates and properties. Their Ṁ − L trend at Galactic metallicity resembles the observed one, favouring therefore more compact solutions with outflows driven by the iron opacitybump. Both Nakauchi & Saio (2018) and Gräfener et al. (2017) make use of a prescribed single betavelocity law for the supersonic layers to constrain the massloss rate of WR stars with optically thick winds, while Ṁ is a free parameter in our hydrodynamic models.
In the supersonic part of WR outflows, the presence of instabilities and inhomogeneities should be taken into account. Invoking density inhomogeneities, or “clumping”, in stellar models has an influence on the mean opacity, due to the enhanced density in clumps (Moffat et al. 1988; Hamann & Koesterke 1998). However, a porous structure could also imply lowered mean opacity, counteracting the effect of smallscale clumping (Shaviv 1998; Oskinova et al. 2007). Assuming that the material near the Feopacity peak is clumped, Gräfener et al. (2012) were able to extend the surface radii of their hydrostatic helium star models, in practical terms, by significantly enhancing the iron bump opacity (see also Appendix A). The surface temperatures of such hydrostatic, strongly inflated models computed with plane parallel grey atmospheres are then compared with the fictitious effective temperatures at τ = 20 of the atmosphere calculations performed by Hamann et al. (2006). This effective temperature at τ = 20 should not be confused with the actual blanketed temperature at the base of WNE wind models, which is about a factor of 2 larger (see e.g. Appendix D). As we show in Sect. 4, strongly inflated hydrodynamic solutions imply supersonic flows already at the base of the inflated envelopes for the typical massloss rates of WNE stars (cf. the plane parallel atmosphere model in Fig. 4). Moreover, it would be inconsistent to think that the iron opacity bump is simultaneously responsible for both the inflation of the envelope and the acceleration of the flow. While the detailed effects of clumping and porosity remain a subject for future research, we concentrated in this work on the homogeneous case as those instabilities are expected to initiate above the sonic point (Sundqvist et al. 2013; Owocki 2015). If this is not the case and clumping is already present at the relatively high densities and optical depth of the sonic point of massive helium star models, it might affect to some extent the local opacity.
Standing waves in helium stars (Glatzel et al. 1993) are also not expected. As already pointed out by Grassitelli et al. (2016a), while investigating pulsations in massive helium star models, mass loss has an inhibitory effect on pulsations already in the presence of inflated envelopes.
7 Conclusions
We investigated the conditions at the sonic point and the subsonic structure of hydrodynamic models for chemically homogeneous massive helium stars, thought to be representative of WR stars in the WNE phase.
For typical massloss rates we find that the outflows of our stellar models, computed with newly implemented boundary conditions at the sonic point, are accelerated to transonic velocities by the momentum provided by the increase in opacity associated with the recombination of the iron and irongroup elements, the iron opacity bump, in the temperature range 160–220 kK or by the helium opacity bump in the range 40–60 kK.
The Eddington factor at the sonic point of our models is very close to unity. This allows us to build a sonic HR diagram relating the luminositytomass ratio of stars with optically thick winds to their sonic point temperature. Knowing the sonic point temperature and relying on available opacity tables, it is possible to predict the massloss rates of stars with optically thick winds from this diagram. We use this to derive a simple relation for the minimum massloss rate for WNE stars necessary to have the outflow accelerated to supersonic velocities by the iron opacity bump only as a function of the luminositytomass ratio. The vast majority of the observed Galactic WNE stars exceed this minimum massloss rate, suggesting therefore that the radiative acceleration in the hot part of the iron opacity bump is responsible for the launch of the outflows of these stars to supersonic velocities. This leads to sonic point radii of the order of one solar radius for WNE stars.
In this paper we shed new light on the WR radius problem and show that the iron opacity bump is key to driving outflows through the sonic point in WNE stars, and hence to determining the massloss rate of the star. Our models may provide a useful reference while adopting the sonic HR diagram to compare stellar structure and atmosphere models of optically thick outflows with observations, and may provide solid inner boundary conditions for stellar atmosphere codes, reducing the complexity of the calculations with simple physical arguments and providing ab initio subsonic structure models.
We have also derived analytic relations to better understand the mechanisms at work when extended lowdensity envelopes appear in the stellar models, characterized by the proximity to the Eddington limit, the socalled envelope inflation. In the envelope of a stellar model in hydrostatic equilibrium the emergence of a large gas pressure scale height is the response to the rapid increase in radiative opacity when approaching an opacity bump in radiationpressuredominated layers close to the Eddington limit.
Acknowledgements
L. G. thanks Andreas Sander, Stan Owocki, Alina Istrate, Luca Fossati, and JeanClaude Passy for the constructivediscussions. J. M. acknowledges funding from a Royal Society–Science Foundation Ireland University Research Fellowship.
Appendix A Envelope inflation criterion
Here we attempt to characterize the structure of inflated envelopes. With BEC we compute three versions of the same 15 M_{⊙} helium star model, but adopting three different OPAL opacity tables corresponding to the opacity of stellar matter at three different metallicities. For illustrative purposes, we impose hydrostatic equilibrium and no mass loss by stellar wind. For this computation we also recover the plane parallel grey atmosphere boundary conditions (Heger 1998; Yoon et al. 2006). The advantage of studying a helium zero age main sequence massive helium star model is that such a model is hot enough that only the iron opacity bump is present below its surface. This allows us to more clearly investigate the response of the outer layers subject onlyto the radiative force arising by this opacity bump.
The density profiles of these newly computed stellar models are shown in Fig. A.1 where the blue model is computed with the opacity table for a metallicity Z = 0, while the orange and the black models are computed for Z = 0.02 and Z = 0.05, respectively.The Z = 0 blue model clearly shows a density profile that monotonically decreases as a function of radius, showing moreover a density scale height H_{ρ} increasingly smaller as one approaches the surface. This model shows how the density profile behaves in absence of a pronounced increase in opacity in the outer subsurface layers, showing no sign of envelope inflation and no corehalo structure below its surface. It isconsidered the reference model.
Fig. A.1 Density profiles in units of g cm^{−3} of a 15 M_{⊙} helium star model with different iron bump opacities. The blue dotted model is computed with the opacity table corresponding to a zero metallicity composition (i.e. without the iron opacity bump), the orange continuous model is computed with the opacity table corresponding to Z = 0.02, while the black dashed model is computed with the opacity table corresponding to Z = 0.05. 
This is not the case for the helium star models with opacities corresponding to the metallicity Z = 0.02 and Z = 0.05. Their density profiles in Fig. A.1 show the characteristic envelope inflation configuration encountered by e.g. Gräfener et al. (2012), Sanyal et al. (2015), and Grassitelli et al. (2016a). In these models the higher iron bump opacity brings the outer layers to the Eddington limit; this leads to the formation of an extended, lowdensity envelope which, for the Z = 0.02, increases the surface radius by approximately 10%, while for the Z = 0.05 case, the increase in radius is more than a factor of 2 compared to the noninflated case. Consequently these two models also find their surfaces at lower effective temperatures, T_{eff}≈130 kK, 120 kK, 80 kK for the metallicities Z = 0.00, 0.02, 0.05, respectively. At first the mildly and the more inflated models in Fig. A.1 both show density profiles with density scale heights decreasing in the core region. However, around 1 R_{⊙} , the two inflated models start to differ from the noninflated model. Their density scale heights start to increase due to the increase in the outwarddirected radiative force in the proximity of the iron opacity bump. This can be understood by writing Eq. (13) in the general hydrostatic case (v = 0) $${c}_{\text{s}}^{2}\frac{\text{d}\text{\hspace{0.17em}}\mathrm{ln}(\rho )}{\text{d}r}=g+{g}_{\text{rad}}{c}_{\text{s}}^{2}\frac{\text{d}\text{\hspace{0.17em}}\mathrm{ln}(T)}{\text{d}r}.$$(A.1)
In the force balance an increase of g_{rad} leads to a larger density scale height, and thus the highopacity models have more shallow density gradients. As the opacity keeps rising in the hot part of the iron opacity bump, the two inflated models flatten their density profiles and reach a minimum, followed by a rapid rise in density. From Eq. (A.1) (or Eq. (13) in the case of v^{2} ≪ gr ), a positive density gradient (i.e. a density inversion) appears when $$\Gamma \gtrsim \frac{1\beta}{1\frac{3}{4}\beta},$$(A.2)
or equivalently $$\frac{{P}_{\text{gas}}}{{P}_{\text{rad}}}\gtrsim 4\frac{1\Gamma}{\Gamma}.N1$$(A.3)
Equaling the sides in the inequality (A.2) gives $$\beta \simeq \frac{1\Gamma}{1\frac{3}{4}\Gamma},$$(A.4)
which defines the condition when either $\frac{\text{d}\rho}{\text{d}r}=0$ (and the flow stays subsonic) or, in the hydrodynamic case, the sonic point β of transonic outflows (see Sect. 2). Moreover, the inequality (A.2) can be compared with the criterion for convection (Langer 1997), $$\Gamma \ge (1\beta )\frac{3224\beta}{3224\beta 3{\beta}^{2}},$$(A.5)
showing that a density inversion is necessarily convective (see also Joss et al. 1973).
Equation (A.4) shows that a null density gradient can appear only for specific sets of β and Γ , with lower β (and thus densities) for higher Γ, as we can see in Fig. A.1. A modified Eddington factor could then be introduced, namely $${\Gamma}_{\text{S}}=\frac{\kappa L}{4\pi cGM}\left(1+\frac{{P}_{\text{gas}}}{4{P}_{\text{rad}}}\right),$$(A.6)
which includes the outwarddirected force associated with the gas temperature gradient, and which defines (when equalto unity) either the beginning of a density inversion or the condition at the sonic point of stars for which convection is inefficient (see Sect. 2 and Eq. (4)). This equation also shows that the local Eddington factor cannot be exactly unity either at the sonic point or when $\frac{\text{d}\rho}{\text{d}r}=0$. In either case, these two quantities tend to differ by only few percent in the radiationpressuredominated outer layers of massive stars, which is in fact equivalent to the approximation $\frac{\text{d}{c}_{\text{s}}^{2}}{\text{d}r}2\frac{{c}_{\text{s}}^{2}}{r}\ll g$ made in Sect. 2 and below.
From Fig. A.1 it is evident that the inflated stellar models differ markedly from the noninflated model when the profiles show an inflection point, i.e. a minimum in d ln (P_{gas})∕dr, or almost equivalently d ln(ρ)∕dr. An increase in the density scale height follows the inflection point, contrary to the monotonically decreasing density scale height of stellar models which do not inflate. This is due to an increasing contribution from g_{rad} in the momentum equation, affecting the balance between the gravitational force and the gas pressure gradient that defines the condition of hydrostatic equilibrium. Analytically, the second derivative of the gas pressure profile writes $$\begin{array}{l}\frac{{\text{d}}^{2}\mathrm{ln}({P}_{\text{gas}})}{\text{d}{r}^{2}}=\frac{{\text{d}}^{2}\mathrm{ln}(\rho )}{\text{d}{r}^{2}}+\frac{{\text{d}}^{2}\mathrm{ln}(T)}{\text{d}{r}^{2}}\\ \text{\hspace{1em}}\text{\hspace{1em}}\text{\hspace{1em}}\text{\hspace{1em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}=\frac{(g{g}_{\text{rad}})}{{c}_{\text{s}}^{2}}\frac{\text{d}\text{\hspace{0.17em}}\mathrm{ln}(T)}{\text{d}r}+\frac{{g}_{\text{rad}}}{{c}_{\text{s}}^{2}}\frac{\text{d}\text{\hspace{0.17em}}\mathrm{ln}(\kappa )}{\text{d}r},\end{array}$$(A.7)
where for simplicity we have neglected the geometrical term as 2∕r ≪dln(T)∕dr (which might not be always possible for main sequence stars). In writing Eq. (A.7) we have assumed a fixed chemical composition and no energy source or sink in the outer layers (strictly true for hydrostatic models in thermal equilibrium, and true to a high degree for steadystate hydrodynamic models).
In Eq. (A.7) the first term on the righthand side is always negative in subEddington regions. Instead, the second term on the righthand side depends on both the radiative acceleration and the opacity profile, and it is positive when the opacity increases with increasing radius. This is the term that becomes dominant at and above the inflection point. For a constant opacity, an inflection point would be present only when Γ = 1. However, in the case of more complex, nonconstant opacity profiles, a rapid increase in opacity at temperatures where the recombination of some elements takes place can lead to an increase in the gas pressure scale height, and thus density scale height, as the second derivative of gas pressure changes sign.
From Eq. (A.7) a criterion for the appearance of the inflection point, i.e. for the beginning of an envelope inflation can be derived: $$\frac{\text{d}\text{\hspace{0.17em}}\text{ln}(\kappa )}{\text{d}\text{\hspace{0.17em}}\text{ln}(T)}\equiv \frac{4\text{\hspace{0.17em}}\text{d}\text{\hspace{0.17em}}\text{ln}(\Gamma )}{\text{d}\text{\hspace{0.17em}}\text{ln}({P}_{\text{rad}})}\lesssim 1\frac{1}{\Gamma},$$(A.8)
We can nowunderstand envelope inflation as the appearance of low gas pressure gradients in response to a steep increase in the radiative opacity in radiationpressuredominated layers close to the Eddington limit. Above the inflection point, the increase in opacity increases the gas pressure scale height ${H}_{{P}_{\text{gas}}}$ following the increase in g_{rad} as Γ → 1 (Sanyal et al. 2017): $$\frac{1}{{H}_{{P}_{\text{gas}}}}:=\frac{\text{dln}({P}_{\text{gas}})}{\text{d}r}=\frac{g{g}_{\text{rad}}}{{c}_{\text{s}}}=\frac{g}{{c}_{\text{s}}}(1\Gamma ).$$(A.9)
From here it follows that as the gas pressure scale height increases, the location of the photosphere is found at larger radii. Analytical expressions for the latter in the case of envelope inflation originating from the iron bump in very massive helium stars and neglecting convection have been given in Gräfener et al. (2012). We note here that stellar models can be inflated, i.e. have extended lowdensity envelopes close to the Eddington limit, but do not show a density inversion (or a gas pressure inversion; see e.g. the 85 M_{⊙} main sequence model in Sanyal et al. 2015). On the other hand, a density inversion appearing when Γ ≈ 1 is always an indication of envelope inflation.
Effects of convection on inflated envelopes
The presence of convection can affect the structure of the inflated envelope (Sanyal et al. 2015). It effectively reduces the radial extent of the inflated layers by affecting the temperature stratification without contributing directly to the force balance. In convective layers, the ratio of the luminosity transported by radiation to the total luminosity can be expressed as $$\frac{{L}_{\text{rad}}}{L}=\frac{\nabla}{{\nabla}_{\text{rad}}},$$(A.10)
where ∇_{rad} is the temperature gradient required to transport all the luminosity by radiation, while ∇ is the actual temperature gradient which includes the contribution from convection. In a convective layer, according to the Schwartzschild criterion, ∇_{rad} > ∇ (Kippenhahn & Weigert 1990). Together with Eq. (A.10), the radiative acceleration from Eq. (3) writes as $${g}_{\text{rad}}^{\text{conv}}=\frac{\kappa L}{4\pi {r}^{2}c}\frac{\nabla}{{\nabla}_{\text{rad}}}<{g}_{\text{rad}},$$(A.11)
where the inequality arises from ∇_{rad} > ∇. Therefore, convection effectively reduces the radiative acceleration by not contributing to the force balance or, equivalently, has the same effect of a reduction of the local opacity as $${\kappa}_{\text{conv}}=\frac{\nabla}{{\nabla}_{\text{rad}}}\kappa ,$$(A.12)
leading to small gas pressure scale height (Eq. (A.9)). The criterion in Eq. (A.8) is still valid, as far as the opacity gradient takes into account the effect of convection.
Appendix B Radiation field at the sonic point
In general, the temperature stratification in a stellar atmosphere is a global problem due to the intrinsic coupling of the radiation field to different regions of the atmosphere. However, we discuss here the limiting case in which the radiative energy transport turns from being a global problem, to a local one, and apply it to our analysis of the conditions at the sonic point of WNE stars.
The intensity of radiation I at the sonic point optical depth τ_{S} > 1 in the Rosseland approximation can be expressed in terms of a Taylor–McLaurin series around the equilibrium Planckian value B(T) of the source function. The source function S can be written as (Mihalas 1978; Hansen & Kawaler 1994) $$S(t)={\displaystyle \sum _{n=0}^{\infty}\frac{{(t{\tau}_{\text{S}})}^{n}}{n!}}\frac{{\text{d}}^{n}B({\tau}_{\text{S}})}{\text{d}{\tau}^{n}},$$(B.1)
with B(τ_{S}) := B(T[τ_{S}]) the Planck function at the temperatureof the sonic point and t a dummy variable. Integration of the equation of radiative transfer $$\frac{\partial I{e}^{\tau /\mu}}{\partial \tau}=\frac{1}{\mu}S(t){e}^{\tau /\mu},$$(B.2)
for the outgoing radiation I^{+}, i.e. from τ = ∞ to τ = τ_{S}, leads to $${I}^{+}({\tau}_{\text{S}},\mu \ge 0)={\displaystyle \sum _{n=0}^{\infty}{\mu}^{n}}\frac{{\text{d}}^{n}B({\tau}_{\text{S}})}{\text{d}{\tau}^{n}},$$(B.3)
with μ := cos(θ) indicating the direction of the beam. The expression for the incoming radiation I^{−} from τ = 0 to τ = τ_{S} writes instead $${I}^{}({\tau}_{\text{S}},\mu \le 0)={\displaystyle \sum _{n=0}^{\infty}{\mu}^{n}}\frac{{\text{d}}^{n}B({\tau}_{\text{S}})}{\text{d}{\tau}^{n}}\left(1{e}^{\frac{{\tau}_{\text{S}}}{\mu}}{\displaystyle \sum _{k=0}^{n}{\left(\frac{{\tau}_{\text{S}}}{\mu}\right)}^{k}}\frac{1}{k!}\right).$$(B.4)
The inwarddirected radiation differs from the outwarddirected radiation by terms of the order of ${e}^{\frac{\tau}{\mu}}$.
We can now more quantitatively refer to WNE stars. At first we can say that, assuming as reference a typical optical depth for the sonic point of the order of τ_{S} ≈ 5 (see e.g. Table 1 and Gräfener & Hamann 2005), the exponential term in Eq. (B.4) contributes for less than 1%. For low n therefore, the term ${e}^{\frac{\tau}{\mu}}{\displaystyle \sum _{k=0}^{n}{\left(\frac{\tau}{\mu}\right)}^{k}}\frac{1}{k!}$ is negligible,as usually assumed in order to derive the diffusive approximation (Mihalas 1978). Considering the momenta of the radiative energy transport equation, namely the energy density, the energy flux, and the radiation pressure, an approximation of the derivatives by appropriate differences shows that the ratio of successive terms in the series is of order O(1/${\tau}_{\text{S}}^{2}$) (see Mihalas 1978). As such, we can study the convergence of the series via $$O\left(\frac{1}{{\tau}_{\text{S}}^{2}}\right)\approx O\left(\frac{{\lambda}^{2}}{{H}_{\text{P}}^{2}}\right),$$(B.5)
where we have taken the pressure scale height H_{P} as the characteristic scale of the system. From Table 1 the pressure scale heights are an order of magnitude larger than the mean free path at the sonic point, assuring the rapid convergence of the series to the Planckian value, with the highorder terms contributing less than 1%. We can thus assert that the conditions at the sonic point of WNE stars, especially those with stellar winds driven by the iron opacity bump, closely approach the LTE conditions.
Consequently, the radiation field at the sonic point of these stars is close to being isotropic, with the radiative transfer losing the explicit dependence on the optical depth in Eq. (B.4). We can further estimate the level of anisotropy of the radiation field at the sonic point in case of outflows driven by the iron opacity bump as (Mihalas 1978) $$\frac{\text{Anisotropicterm}}{\text{Isotropicterm}}\propto {\left(\frac{{T}_{\text{eff}}{}_{\text{S}}}{{T}_{\text{S}}}\right)}^{4}\approx 5\%,$$(B.6)
with the T_{eff}_{S} fictitious effective temperature at the sonic point derived from the Stefan–Boltzmann law, having adopted a radius of 1 R_{⊙} and a typical luminosity of 10^{5} L_{⊙} . We can also estimate the proximity to the LTE conditions and the validity of the diffusive approximation by comparing the radiation timescale associated with the photon thermalization $${\tau}_{t}=\frac{\lambda}{c},$$(B.7)
to the other macroscopic timescales, i.e. the dynamical and thermal timescales of the envelopes (Thorne et al. 1981; Thorne 1981). The timescale for photon thermalization, i.e. the time needed by the radiation field to approach its LTE value, is of the order of seconds, while the dynamical and thermal timescales were estimated to be of the order of hundreds of seconds for massive helium stars (Grassitelli et al. 2016a).
All these considerations lead us to conclude that the states of the gas and radiation are close to their LTE conditions. In LTE, it is assumed that the local T and ρ are sufficient to determine the ionization state of the gas, all microscopic processes being close to detailed balance. This approximation is justified as far as the thermalization depth of radiation, i.e. the depth at which the source function approaches its equilibrium value, is smaller than the sonic point optical depth (Mihalas 1978; Pistinner & Eichler 1995).
Appendix C Nozzle analogy and critical point
The critical solution describing the kinematics of stellar winds can be understood in analogy with an ideal rocket (or de Laval) nozzle, i.e. a nozzle with walls that narrow at first, reach a throat, and then rapidly expand (Abbott 1980; Bjorkman 1995; Lamers & Cassinelli 1999; Shore 2007). Given a compressible steadystate flow of gas, the subsonic flow accelerates in the converging part of the nozzle due to mass conservation. If the initial conditions are such that the flow can reach supersonic velocities before passing the throat, the sonic point is always located exactly at the throat of the nozzle. This system, common to almost all modern rockets (Raga & Canto 1989), shows how a smooth steadystate flow can selfadjust such that a transonic flow can be established via information that travels at the speed of sound. Once the transonic flow is established, the conditions in the diverging supersonic part define the local kinematics of the flow, such as its velocity and temperature, but they can influence neither the subsonic part of the flow nor the massflow rate (Lamers & Cassinelli 1999; Shore 2007; Maciel 2014).
The analogy with a stellar wind lies in the form of the momentum equation of a de Laval nozzle, that is (Lamers & Cassinelli 1999) $${v}^{2}\frac{\text{dln}v}{\text{d}x}={c}_{\text{s}}^{2}\frac{\text{dln}\rho}{\text{d}x}\frac{\text{d}{c}_{\text{s}}^{2}}{\text{d}x},$$(C.1)
for an ideal gas flow, with x spatial coordinate. Combining this with the continuity equation dln(v) + dln(ρ) + dln(S) = 0, where S corresponds to the area of a section of the nozzle, the force balance is given by $$\left({v}^{2}{c}_{\text{s}}^{2}\right)\frac{\text{dln}v}{\text{d}x}={c}_{\text{s}}^{2}\frac{\text{dln}S}{\text{d}x}\frac{d{c}_{\text{s}}^{2}}{\text{d}x}.$$(C.2)
Comparing Eq. (C.2) with Eq. (1), the pseudonozzle term in the context of radiationpressuredriven winds is $${c}_{\text{s}}^{2}\frac{\text{dln}S}{\text{d}x}=g\left(1\Gamma +\frac{2{c}_{\text{s}}^{2}}{gr}\right),$$(C.3)
which showshow, in the limit of $2{c}_{\text{s}}^{2}\ll gr$, the Eddington limit acts similarly to a nozzle term for a radiationdriven stellar outflow.
The common property between the hydrodynamic of a flow through a pipe and a stellar wind is the readjustment to the characteristic solution. It is the solution that starts subsonic, goes through the critical point of the momentum equation (where the massflow is defined, Lamers & Cassinelli 1999), and has a finite velocity at infinity (Abbott 1980; Nugis & Lamers 2002; Shore 2007). In the nozzle the critical point also defines the last point of the system which can still communicate with all the regions of the flow. By reducing the perturbative analysis in e.g. Abbott (1980) to a case where the line opacity and acceleration are independent from the velocity profile, we can easily show that the characteristic speed of the hyperbolic system of stellar wind equations is the local sound speed^{5} . However, as noticed by Thorne et al. (1981) while investigating the event horizon, the diffusive approximation for radiative transfer is notoriously acausal (due to the parabolic nature of the partial differential equation governing the phenomenon). This shortcoming implies that information can be transmitted instantaneously. Hence a relativistic formulation governing diffusion together with the introduction of the second sound speed^{6} would be more appropriate (Flammang 1982; Mihalas & Mihalas 1984; Mandelis et al. 2001; Ali & Zhang 2005). If the speed of heat were, in this context, the same as or smaller than the speed of sound, a sonic horizon would arise in this class of stars, with the sonic point being the last point that could communicate with all the regions of the flow.
Effects of velocity gradients on the critical point
Equation 1 shows that the sonic point is the critical point of a stellar wind, but this is true only if no other term with a dependency on dv∕dr contributes to the righthand side of Eq. (1). In the context of linedriven winds, Castor et al. (1975) and later applications of their theory (the CAK theory) showed the importance of lines unsaturated due to Doppler shift in driving the winds of massive stars. In fact, in the presence of steep velocity gradients the rest wavelength of each accelerating element in the outflow is Doppler shifted such that unattenuated continuum can be absorbed and momentum transfer is very efficient. This in turn allows massive stars to have winds with very high terminal wind velocities. In this context the opacity arising from the spectral lines includes a dependency on the velocity gradient, i.e. κ(ρ, T, dv∕dr, ...), and therefore Abbott (1980) argue that the sonic point is no longer, strictly speaking, the critical point of the outflow.
However, Nugis & Lamers (2002) suggest that the contribution to the opacity by the Doppler shifted lines in hot WR stars is usually less than 1% at typical sonic point conditions. This can be quantified by investigating the inverse of the optical depth parameter (Table 1). This quantity compares the Doppler shift of the lines as being due to the velocity gradient through a mean free path to the thermal velocity of the protons (Nugis & Lamers 2002; Ro & Matzner 2016). They estimate this quantity to be typically of the order of 100 for WR stars. Independently of this, Lucy (2007a,b) also suggested that, in contrast to the standard CAK ansatz, the sonic point retains its role of critical point for the system of differential equations describing the system (see also Lucy 1975, 1998; Poe et al. 1990; Müller & Vink 2008; Sander et al. 2017). The importance of the unsaturated lines in the highly supersonic outflow is instead clear (Castor et al. 1975; Abbott 1980; Abbott & Lucy 1985; Kudritzki et al. 1989; Gayley & Owocki 1995; Gräfener & Hamann 2005; Vink & de Koter 2005; Puls et al. 2008; Gräfener & Hamann 2008).
Appendix D WR111
We adopt the hydrodynamically selfconsistent wind structure derived for WR111, a WC5 WR star, by Gräfener & Hamann (2005) and connect it smoothly to our hydrodynamic model. This is done in order to show how our calculations of the subsonic structure can be smoothly matched to the sonic point conditions of a hydrodynamical wind model of a postmain sequence WR star and also to give a picture of what the internal profile of such a star looks like. This is done in Fig. D.1, where the velocity and temperature profiles of a 12 M_{⊙} stellar model with sonic point boundary conditions are plotted together with the supersonic structure of the wind (Gräfener & Hamann 2005). The stellar model has been built such that its sonic point conditions, namely sonic point temperature, density, radius, luminosity, and massloss rate, match those derived for WR111. Luminosity and mass loss for the stellar model are, following Gräfener & Hamann (2005), log (L∕ L_{⊙} ) = 5.45 and log (Ṁ) ≈−5.1, respectively.Homogeneous zero age main sequence helium star models with such luminosity were found to have a core radius that was too large to reproduce the radius of this wind model. Therefore, given that WR111 is of the WC class and it is consequently thought to be in a more advanced stage of its evolution, we used a more evolved model (with approximately 0.5 carbon mass fraction in its centre) instead of a homogeneous helium star.The wind model by Gräfener & Hamann (2005) found its sonic point at an optical depth of τ = 5.4 with temperature T_{S} ≈ 200 kK and log (ρ_{S}) ≈ 8.5. For this star Gräfener & Hamann (2005) finds a T_{*}, i.e. the effective temperature at τ = 20, equal to 140 kK. Instead, for the same object, spectral analysis conducted by Gräfener et al. (2002) and Sander et al. (2012) assuming a betavelocity law led to T_{*} of the orderof 85 kK. The corresponding radius derived via the Stefan–Boltzmann law would thus be ≈ 2.5 R_{⊙} , different by more than a factor of 2 from the actual quasihydrostatic radius in Fig. D.1. This shows how misleading it might be to use the Stefan–Boltzmann law at τ = 20 and, more importantly, how the adoption of a betavelocity law to describe the dynamics of the optically thick winds of WR stars might lead to a significantly different radius estimate compared to hydrodynamic calculations.
Fig. D.1 Temperature (in K) and velocity profile (in km s^{−1}) as a function of radius for the WR111 stellar structure and wind model. The orange profile indicates the subsonic structure derived with our newly developed models, while the blue profile is the supersonic wind profile from Gräfener & Hamann (2005). 
References
 Abbott, D. C. 1980, ApJ, 242, 1183 [NASA ADS] [CrossRef] [Google Scholar]
 Abbott, D. C., & Conti, P. S. 1987, ARA&A, 25, 113 [NASA ADS] [CrossRef] [Google Scholar]
 Abbott, D. C., & Lucy, L. B. 1985, ApJ, 288, 679 [NASA ADS] [CrossRef] [Google Scholar]
 Ali, Y. M., & Zhang, L. C. 2005, Int. J. Heat Mass Trans., 48, 2397 [CrossRef] [Google Scholar]
 Bestenlehner, J. M., Gräfener, G., Vink, J. S., et al. 2014, A&A, 570, A38 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Bjorkman, J. E. 1995, ApJ, 453, 369 [NASA ADS] [CrossRef] [Google Scholar]
 Blaes, O., & Socrates, A. 2003, ApJ, 596, 509 [NASA ADS] [CrossRef] [Google Scholar]
 BöhmVitense, E. 1958, ZAp, 46, 108 [NASA ADS] [Google Scholar]
 Bondi, H. 1952, MNRAS, 112, 195 [NASA ADS] [CrossRef] [MathSciNet] [Google Scholar]
 Castor, J. I., Abbott, D. C., & Klein, R. I. 1975, ApJ, 195, 157 [NASA ADS] [CrossRef] [Google Scholar]
 Chené, A.N., Moffat, A. F. J., Cameron, C., et al. 2011, ApJ, 735, 34 [NASA ADS] [CrossRef] [Google Scholar]
 Chiosi, C., & Maeder, A. 1986, ARA&A, 24, 329 [NASA ADS] [CrossRef] [Google Scholar]
 Crowther, P. A. 2008, in Massive Stars as Cosmic Engines, eds. F. Bresolin, P. A. Crowther, & J. Puls (IAU Symposium), 250, 47 [NASA ADS] [Google Scholar]
 Dessart, L., & Owocki, S. P. 2003, A&A, 406, L1 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Dessart, L., & Owocki, S. P. 2005, A&A, 432, 281 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Flammang, R. A. 1982, MNRAS, 199, 833 [NASA ADS] [CrossRef] [Google Scholar]
 Flammang, R. A. 1984, MNRAS, 206, 589 [NASA ADS] [CrossRef] [Google Scholar]
 Gayley, K. G., & Owocki, S. P. 1995, ApJ, 446, 801 [NASA ADS] [CrossRef] [Google Scholar]
 Georgy, C., Ekström, S., Meynet, G., et al. 2012, A&A, 542, A29 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Glatzel, W., Kiriakidis, M., & Fricke, K. J. 1993, MNRAS, 262, L7 [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., & Hamann, W.R. 2008, Mass Loss from Stars and the Evolution of Stellar Clusters, eds. A. de Koter, L. J. Smith, & L. B. F. M. Waters, ASP Conf. Ser., 388,21 [NASA ADS] [Google Scholar]
 Gräfener, G., Koesterke, L., & Hamann, W.R. 2002, A&A, 387, 244 [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., Fossati, L., SimónDiáz, S., et al. 2015, ApJ, 808, L31 [NASA ADS] [CrossRef] [Google Scholar]
 Grassitelli, L., Chene’, A.N., Sanyal, D., et al. 2016a, A&A, 590, A12 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Grassitelli, L., Fossati, L., Langer, N., et al. 2016b, A&A, 593, A14 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Hamann, W.R., & Koesterke, L. 1998, A&A, 335, 1003 [NASA ADS] [Google Scholar]
 Hamann, W.R., Gräfener, G., & Liermann, A. 2006, A&A, 457, 1015 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Hansen, C. J.,& Kawaler, S. D. 1994, Stellar Interiors. Physical Principles, Structure, and Evolution (Berlin Heidelberg: SpringerVerlag) [Google Scholar]
 Heger, A. 1998, PhD thesis, MaxPlanck Institute for Astrophysics [Google Scholar]
 Heger, A., & Langer, N. 1996, A&A, 315, 421 [NASA ADS] [Google Scholar]
 Heger, A., Langer, N., & Woosley, S. E. 2000, ApJ, 528, 368 [NASA ADS] [CrossRef] [Google Scholar]
 Hillier, D. J. 1987, ApJS, 63, 947 [NASA ADS] [CrossRef] [Google Scholar]
 Hillier, D. J. 2015, in WolfRayet Stars: Proceedings of an International Workshop held in Potsdam, eds. W.R. Hamann, A. Sander, & H. Todt (Potsdam: Potsdam University Press), 65 [Google Scholar]
 Hubeny, I., & Mihalas, D. 2014, Theory of Stellar Atmospheres [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]
 Jiang, Y.F., Cantiello, M., Bildsten, L., Quataert, E., & Blaes, O. 2015, ApJ, 813, 74 [NASA ADS] [CrossRef] [Google Scholar]
 Joss, P. C., Salpeter, E. E., & Ostriker, J. P. 1973, ApJ, 181, 429 [NASA ADS] [CrossRef] [Google Scholar]
 Kato, M., & Iben, Jr. I. 1992, ApJ, 394, 305 [NASA ADS] [CrossRef] [Google Scholar]
 Kippenhahn, R., & Weigert, A. 1990, Stellar Structure and Evolution, XVI (Berlin Heidelberg: SpringerVerlag) [CrossRef] [Google Scholar]
 Köhler, K., Langer, N., de Koter, A., et al. 2015, A&A, 573, A71 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Kozyreva, A., Yoon, S.C., & Langer, N. 2014, A&A, 566, A146 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Kudritzki, R. P., Pauldrach, A., Puls, J., & Abbott, D. C. 1989, A&A, 219, 205 [NASA ADS] [Google Scholar]
 Lamers, H. J. G. L. M., & Cassinelli, J. P. 1999, Introduction to Stellar Winds (Cambridge: Cambridge University Press) [CrossRef] [Google Scholar]
 Lamers, H. J. G. L. M., Maeder, A., Schmutz, W., & Cassinelli, J. P. 1991, ApJ, 368, 538 [NASA ADS] [CrossRef] [Google Scholar]
 Langer, N. 1989, A&A, 210, 93 [NASA ADS] [Google Scholar]
 Langer, N. 1997, in Luminous Blue Variables: Massive Stars in Transition, eds. A. Nota & H. Lamers, ASP Conf. Ser., 120, 83 [NASA ADS] [Google Scholar]
 Langer, N. 2012, ARA&A, 50, 107 [NASA ADS] [CrossRef] [Google Scholar]
 Langer, N., Kiriakidis, M., El Eid, M. F., Fricke, K. J., & Weiss, A. 1988, A&A, 192, 177 [NASA ADS] [Google Scholar]
 Langer, N., Hamann, W.R., Lennon, M., et al. 1994, A&A, 290, 819 [NASA ADS] [Google Scholar]
 Langer, N., & Kudritzki, R. P. 2014, A&A, 564, A52 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Lucy, L. B. 1975, Mem. Soc. Roy. Sci. Liège, 8, 359 [NASA ADS] [Google Scholar]
 Lucy, L. B. 1998, in Cyclical Variability in Stellar Winds, ed. L. Kaper & A. W. Fullerton, 16 [CrossRef] [Google Scholar]
 Lucy, L. B. 2007a, A&A, 468, 649 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Lucy, L. B. 2007b, A&A, 474, 701 [NASA ADS] [CrossRef] [EDP Sciences] [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]
 Maciel, W. J. 2014, Hydrodynamics and Stellar Winds: An Introduction (Berlin: SpringerVerlag) [CrossRef] [Google Scholar]
 Mandelis, A., Nicolaides, L., & Chen, Y. 2001, Phys. Rev. Lett., 87, 020801 [NASA ADS] [CrossRef] [Google Scholar]
 McClelland, L. A. S., & Eldridge, J. J. 2016, MNRAS, 459, 1505 [NASA ADS] [CrossRef] [Google Scholar]
 Meier, D. L. 1982, ApJ, 256, 681 [NASA ADS] [CrossRef] [Google Scholar]
 Michaux, Y. J. L., Moffat, A. F. J., Chené, A.N., & StLouis, N. 2014, MNRAS, 440, 2 [NASA ADS] [CrossRef] [Google Scholar]
 Mihalas, D. 1978, Stellar atmospheres, 2nd edn. (San Francisco: W. H. Freeman and Co.) [Google Scholar]
 Mihalas, D., & Mihalas, B. W. 1984, Foundations of radiation hydrodynamics (New York: Oxford University Press) [Google Scholar]
 Moffat, A. F. J., Drissen, L., Lamontagne, R., & Robert, C. 1988, ApJ, 334, 1038 [NASA ADS] [CrossRef] [Google Scholar]
 Moncrief, V. 1980, ApJ, 235, 1038 [NASA ADS] [CrossRef] [Google Scholar]
 Mullan, D. J. 1984, ApJ, 283, 303 [NASA ADS] [CrossRef] [Google Scholar]
 Müller, P. E., & Vink, J. S. 2008, A&A, 492, 493 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Nakauchi, D., & Saio, H. 2018, ApJ, 852, 126 [NASA ADS] [CrossRef] [Google Scholar]
 Neo, S., Miyaji, S., Nomoto, K., & Sugimoto, D. 1977, PASJ, 29, 249 [NASA ADS] [Google Scholar]
 Nobili, L., Turolla, R., & Zampieri, L. 1991, ApJ, 383, 250 [NASA ADS] [CrossRef] [Google Scholar]
 Nugis, T., & Lamers, H. J. G. L. M. 2000, A&A, 360, 227 [NASA ADS] [Google Scholar]
 Nugis, T., & Lamers, H. J. G. L. M. 2002, A&A, 389, 162 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Oskinova, L. M., Hamann, W.R., & Feldmeier, A. 2007, A&A, 476, 1331 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Owocki, S. 2013, Stellar Winds, ed. T. D. Oswalt & M. A. Barstow, 735 [Google Scholar]
 Owocki, S. P. 2015, in Astrophysics and Space Science Library, ed. J. S. Vink, 412, 113 [Google Scholar]
 Owocki, S. P., & Gayley, K. G. 1997, in Luminous Blue Variables: Massive Stars in Transition, eds. A. Nota & H. Lamers, ASP Conf. Ser., 120, 121 [NASA ADS] [Google Scholar]
 Owocki, S. P., & Rybicki, G. B. 1984, ApJ, 284, 337 [NASA ADS] [CrossRef] [Google Scholar]
 Parker, E. N. 1958, ApJ, 128, 664 [Google Scholar]
 Parker, E. N. 1966, ApJ, 143, 32 [NASA ADS] [CrossRef] [Google Scholar]
 Petrovic, J., Pols, O., & Langer, N. 2006, A&A, 450, 219 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Pistinner, S., & Eichler, D. 1995, ApJ, 454, 404 [NASA ADS] [CrossRef] [Google Scholar]
 Poe, C. H., Owocki, S. P., & Castor, J. I. 1990, ApJ, 358, 199 [NASA ADS] [CrossRef] [Google Scholar]
 Puls, J., Vink, J. S., & Najarro, F. 2008, A&A Rev., 16, 209 [Google Scholar]
 Raga, A. C., & Canto, J. 1989, ApJ, 344, 404 [NASA ADS] [CrossRef] [Google Scholar]
 Ro, S., & Matzner, C. D. 2016, ApJ, 821, 109 [NASA ADS] [CrossRef] [Google Scholar]
 Sander, A., Hamann, W.R., & Todt, H. 2012, A&A, 540, A144 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Sander, A., Hamann, W.R., Hainich, R., Shenar, T., & Todt, H. 2015, in WolfRayet Stars: Proceedings of an International Workshop held in Potsdam, eds. W.R. Hamann, A. Sander, & H. Todt, (Potsdam: Potsdam University Press) 139 [Google Scholar]
 Sander, A. A. C., Hamann, W.R., Todt, H., Hainich, R., & Shenar, T. 2017, A&A, 603, A86 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Sanyal, D., Grassitelli, L., Langer, N., & Bestenlehner, J. M. 2015, A&A, 580, A20 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Sanyal, D., Langer, N., Szécsi, D., C Yoon, S., & Grassitelli, L. 2017, A&A, 597, A71 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Schaerer, D. 1996, A&A, 309, 129 [NASA ADS] [Google Scholar]
 Schmutz, W., Leitherer, C., & Gruenwald, R. 1992, PASP, 104, 1164 [NASA ADS] [CrossRef] [Google Scholar]
 Schootemeijer,A., & Langer, N. 2017, ArXiv eprints [arXiv: 1709.08727] [Google Scholar]
 Shaviv, N. J. 1998, ApJ, 494, L193 [NASA ADS] [CrossRef] [Google Scholar]
 Shore, S. N. 2007, Astrophysical Hydrodynamics: An Introduction (Hoboken: Wiley) [CrossRef] [Google Scholar]
 StLouis, N., Chené, A.N., Schnurr, O., & Nicol, M.H. 2009, ApJ, 698, 1951 [NASA ADS] [CrossRef] [Google Scholar]
 Sundqvist, J. O., & Owocki, S. P. 2013, MNRAS, 428, 1837 [NASA ADS] [CrossRef] [Google Scholar]
 Sundqvist, J. O., Petit, V., Owocki, S. P., et al. 2013, MNRAS, 433, 2497 [NASA ADS] [CrossRef] [Google Scholar]
 Tamazawa, S., Toyama, K., Kaneko, N., & Ono, Y. 1975, Ap&SS, 32, 403 [NASA ADS] [CrossRef] [Google Scholar]
 Thorne, K. S. 1981, MNRAS, 194, 439 [NASA ADS] [Google Scholar]
 Thorne, K. S., Flammang, R. A., & Zytkow, A. N. 1981, MNRAS, 194, 475 [NASA ADS] [CrossRef] [Google Scholar]
 Vink, J. S. 2015, in Very Massive Stars in the Local Universe (Basel: Springer International Publishing), Astrophys. Space Sci. Lib., 412 [Google Scholar]
 Vink, J. S., & de Koter A. 2005, A&A, 442, 587 [NASA ADS] [CrossRef] [EDP Sciences] [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]
 Visser, M. 1998, Classical and Quantum Gravity, 15, 1767 [Google Scholar]
 Yoon, S.C., Langer, N., & Norman, C. 2006, A&A, 460, 199 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Similar to the spectroscopic HR diagram (Langer & Kudritzki 2014).
Abbott (1980) introduced the concept of effective sound speed and radiativeacoustic waves while treating the CAK critical point.
The second sound speed is the speed at which a perturbation in the heat flux travels, i.e. the speed of heat. It arises when the Fourier equation for diffusive transport gets a hyperbolic form via a relativistic description of the phenomena (Ali & Zhang 2005).
All Tables
All Figures
Fig. 1 Velocity as a function of radius for a set of 15 M_{⊙} stellar models with different massloss rates and boundary conditions. The black continuous lines (from right to left) indicate models with Ṁ equal to 5, 8,10, 20, 30, 50 × 10^{−6} M_{⊙} yr^{−1} and the pink dashed line indicates the model with Ṁ = 4 × 10^{−6} M_{⊙} yr^{−1}, all having sonic point boundary conditions. The orange dotdashed lines indicates the model with Ṁ = 1 × 10^{−5} M_{⊙} yr^{−1} computed with plane parallel grey atmosphere boundary conditions. The green dashed line indicates the local isothermal sound speedfrom the pink dashed model. 

In the text 
Fig. 2 Density profile (top, in units of g cm^{−3}) and opacity profile (bottom, in units of cm^{2} g^{−1}) as a functionof radius for the same 15 M_{⊙} stellar models with different massloss rates and boundary conditions as in Fig. 1. The colours are as in Fig. 1, and the red circles indicate the location of the sonic point and the green dashed line the Eddington opacity (as defined in Eq. (24) in Sect. 4) of the pink dashed model. 

In the text 
Fig. 3 Velocity profile as a function of radius for a set of 10 M_{⊙} stellar models with different massloss rates (as in Fig. 1). The black continuous lines indicate models with Ṁ ≥ 2 × 10^{−5} M_{⊙} yr^{−1} (2, 3, 5 × 10^{−5}), while the pink dashed lines indicate models with Ṁ ≤ 1 × 10^{−5} M_{⊙} yr^{−1} (5, 8, 10 × 10^{−6}). The green dashed line and the orange dotdashed line indicate the local isothermal sound speed of the model with Ṁ = 5 × 10^{−6} M_{⊙} yr^{−1} and the plane parallel grey atmosphere model with Ṁ = 8 × 10^{−6} M_{⊙} yr^{−1}, respectively. 

In the text 
Fig. 4 Velocity profile as a function of radius for a set of 20 M_{⊙} stellar models with different massloss rates (as in Fig. 1). The black continuous lines indicate models with Ṁ ≥ 5 × 10^{−6} M_{⊙} yr^{−1} (5, 10, 20, 30, 50 × 10^{−6}), while the orange dotdashed lines indicate models with Ṁ = 8 × 10^{−6} M_{⊙} yr^{−1} computed with plane parallel grey atmosphere. The green dashed line indicates the local isothermal sound speed. 

In the text 
Fig. 5 Contour plot showing the opacity (in cm^{2} g^{−2} ) from OPAL tables (colourcoded, see scale at right) as a function of temperature (in K) and of density (in g cm^{−3} ) for a chemical composition of [X, Y, Z] = [0, 0.98, 0.02], with X hydrogen mass fraction, Y helium mass fraction, and Z metallicity. The black lines correspond to the outer layers of the hot, more compact 15 M_{⊙} stellar models with Ṁ ≥ 8 × 10^{−6} M_{⊙} yr^{−1} from Fig. 1 (see also Table 1), while the pink line indicates the cooler model with Ṁ = 4 × 10^{−6} M_{⊙} yr^{−1}, and the orange dotdashed line the plane parallel grey atmosphere model with Ṁ = ×10^{−5} M_{⊙} yr^{−1}. The yellow circles indicate the locations of the sonic points of this set of models, while the white lines indicate the isocontours of the Eddington opacity for the 10 M_{⊙} (dotted), 15 M_{⊙} (continuous), and 20 M_{⊙} (dashed) stellar models. Grey areas indicate regions in the log(ρ) − log(T) diagramwhere the sonic point of a radiationdriven wind cannot be located (d κ∕d r ≤ 0). 

In the text 
Fig. 6 Contour plot showing the massloss rate (colourcoded in units of M_{⊙} yr^{−1}) as a function of sonic point temperature (in K) and of the luminositytomass ratio (in L_{⊙} / M_{⊙} ) of a star. Contours are derived based on the OPAL tables for a chemical composition [X,Y,Z]=[0,0.98,0.02], with X hydrogen mass fraction, Y helium mass fraction, and Z metallicity. 

In the text 
Fig. 7 Contour plot showing the luminositytomass ratio (colourcoded in units of L_{⊙} ∕ M_{⊙} ) as a function of sonic point temperature (in K) and of the massloss rate (in M_{⊙} yr ^{−1}) of a star. Contours are derived based on the OPAL tables for a chemical composition [X,Y,Z]=[0,0.98,0.02], with X hydrogen mass fraction, Y helium mass fraction, and Z metallicity. 

In the text 
Fig. 8 Minimum luminositytomass ratio (in units of L_{⊙} / M_{⊙} , left vertical axis) and corresponding luminosity according to the mass–luminosity relation from Langer (1989, right vertical axis), as a function of massloss rate (in M_{⊙} yr ^{−1}) for an outflow to be driven to supersonic velocities by the radiative acceleration from the iron opacity bump. The green symbols indicate Hfree WNE stars from Hamann et al. (2006), with the different shapes indicating different spectral subclasses. 

In the text 
Fig. 9 Sonic HR diagram relating the luminosity (in L_{⊙} ) to the sonic point temperature (in K). The background is colourcoded according to the massloss rates by stellar wind (legend on the right, in units of M_{⊙} yr ^{−1}) predicted via the use of the OPAL opacity tables and the proximity of the sonic point to the Eddington limit (see Fig. 6). The black symbols are observed Galactic WNE stars from Hamann et al. (2006), located in this diagram viatheir observed luminosity and massloss rates and with symbols indicating their spectral subclass, as in Fig. 8. The white dashed line is the best fitting linear relation for the observed WNE stars. 

In the text 
Fig. A.1 Density profiles in units of g cm^{−3} of a 15 M_{⊙} helium star model with different iron bump opacities. The blue dotted model is computed with the opacity table corresponding to a zero metallicity composition (i.e. without the iron opacity bump), the orange continuous model is computed with the opacity table corresponding to Z = 0.02, while the black dashed model is computed with the opacity table corresponding to Z = 0.05. 

In the text 
Fig. D.1 Temperature (in K) and velocity profile (in km s^{−1}) as a function of radius for the WR111 stellar structure and wind model. The orange profile indicates the subsonic structure derived with our newly developed models, while the blue profile is the supersonic wind profile from Gräfener & Hamann (2005). 

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.