© ESO, 2017

1. Introduction

The strong winds of emission-line stars like Wolf-Rayet (WR) stars or luminous blue variables (LBVs) are of key importance for the mechanical and chemical feedback of massive stars on the interstellar medium. They decide how much mass and angular momentum massive stars lose before their final collapse, and whether the most massive stars end their lives as neutron stars or black holes.

The occurrence of emission lines in the spectra of WR stars and LBVs is related to the high optical thickness of their winds. First of all, this means that at short wavelengths, ionising photons are efficiently absorbed by optically thick continua within the extended wind. The subsequent recombination then leads to the observed emission-line cascades at longer wavelengths. Secondly, this means that photons can be absorbed and re-emitted several times before they escape the wind (Lucy & Abbott 1993). In this way, the radiation field can transfer a multiple of its photon momentum to the wind material, leading to wind-efficiency numbers η = Ṁν_∞/ (L/c) > 1, where η denotes the ratio between the wind momentum rate (the product of the mass-loss rate Ṁ and the terminal wind velocity ν_∞) and the momentum of the radiation field (the stellar luminosity L divided by the speed of light c). As we discuss later in this work, this also means that the sonic point of optically thick winds is located at flux-mean optical depths τ_s > 1 (cf. Vink & Gräfener 2012).

In the limit of a high optical depth τ_s ≫ 1 and for sufficiently small velocity gradients (cf. Eq. (1)), the physics of radiation-driven winds simplifies because the radiative transfer can be described in the diffusion limit, and the flux-mean opacity κ_F equals the Rosseland-mean opacity κ_R(ρ,T), which is a function of density ρ and temperature T only. In this limit, the equation of motion has a critical point at the sonic point that allows inferring some properties of winds in this regime. The critical conditions arising in this limit for the winds of WR stars have been investigated by Nugis & Lamers (2002).

For WR stars, Lamers & Nugis (2002) distinguished between hot and cool winds depending on the temperature at their sonic point. They identified a hot regime with sonic-point temperatures ≳160 kK slightly above the temperature of the Fe-opacity peak, and a cool regime with temperatures of 40–70 kK related to a cooler opacity bump. Detailed hydrodynamic atmosphere/wind models for hot WR stars of the carbon sequence (WC stars, Gräfener & Hamann 2005) and cool H-rich WR stars of the nitrogen sequence (WNh stars, Gräfener & Hamann 2008) confirmed these temperature ranges, although there may be an ambiguity whether the opacities in the cool temperature range originate from He, C, or Fe. While Nugis & Lamers (2000), Ro & Matzner (2016) identified He and C (the latter in carbon-rich mixtures) as the main opacity source, Lamers & Nugis (2002), Gräfener & Hamann (2008) identified Fe as the main contributor. In the non-LTE models of Gräfener & Hamann, these differences are most likely caused by the moderate wind optical depths of WNh stars, and the consequent break-down of the diffusion approximation in this regime.

These studies further showed that it is possible to drive the outer part of optically thick, WR-type winds by radiation pressure alone, and that optically thick winds are formed when stars are approaching the Eddington limit. The resulting dependence of the mass-loss rate on the classical Eddington factor Γ_e, where Γ_e = κ_eL/ (4πcGM) denotes the Eddington factor defined with respect to the free-electron opacity κ_e, has been confirmed by means of Monte Carlo models by Vink et al. (2011), who identified a kink in the Γ_e-dependence between the optically thick and optically thin regime.

WR-type winds are indeed found for a variety of objects with high Eddington factors mainly because of their high L/M ratios. Evolved massive stars, such as H-free WN and WC stars, reach L/M ratios of the order of 10⁴L_⊙/M_⊙ because of the enhanced mean molecular weight in their He-burning core. [WC]-type central stars of planetary nebulae reach high L/M ratios as a consequence of shell burning directly above their electron-degenerate core (cf. Gräfener et al. 2008). WNh stars, on the other hand, are main-sequence stars with masses ≳100 M_⊙. They reach high L/M ratios because L increases much more steeply than M for increasing stellar masses.

A related effect that occurs near the Eddington limit is the envelope inflation effect (Ishii et al. 1999; Petrovic et al. 2006; Gräfener et al. 2012; Sanyal et al. 2015). This effect leads to a radial extension of the outer stellar envelope with a density inversion (Joss et al. 1973), and may be responsible for the large radii that are deduced in spectroscopic analyses of WR stars (e.g. Hamann et al. 2006; Sander et al. 2012; Hainich et al. 2014). The discrepancy between the empirical radii of H-free WR stars and those that are expected for stars on the He-main sequence can amount to up to a factor 10, and is known as the WR radius problem. Gräfener et al. (2012) found that an enhanced mean opacity that is due to an inhomogeneous (clumped) structure within the inflated envelope can possibly account for the large observed radii. The presence of inhomogeneities in similar radiation-dominated envelopes is supported by three-dimensional models of Jiang et al. (2015), however, these authors argued that the effects of porosity in an inhomogeneous medium can reduce the mean opacity (cf. Shaviv 1998; Owocki et al. 2004; Oskinova et al. 2007; Owocki 2008).

Based on a semi-empirical analysis, Gräfener & Vink (2013) investigated the sonic-point conditions that follow from the observed properties of the outer winds of WC stars. They found sonic-point temperatures and densities that would imply enhanced opacities for the majority of stars, which may be the consequence of a clumped envelope structure. Their results suggested that a small number of stars may have smooth envelopes with a compact structure and a hot sonic-point temperature, while the majority has clumped envelopes with an inflated structure and cool sonic-point temperatures.

It is currently not entirely clear how such a clumped and inflated structure can be realised in nature. In many cases, an inflated envelope structure would suggest densities that are so low that the velocities within the envelope would reach values comparable to those in a stellar wind because of the equation of continuity (cf. Petrovic et al. 2006; Ro & Matzner 2016). Except for very low mass-loss rates, the models in both works have small sonic-point radii, suggesting a compact rather than an inflated envelope structure.

To reconcile these results with the radii derived from spectroscopic analyses, it may be necessary to invoke more exotic scenarios for the layers above the sonic point. Petrovic et al. (2006), Ro & Matzner (2016) both found that the wind acceleration ceases beyond the Fe-opacity bump. Gräfener & Hamann (2005) also reported problems with the wind acceleration in the same region. It may thus be that inflated envelopes have more in common with a failed wind, or that there are multiple components of material that are partly falling back on the star.

In the present work, we intend to test how the assumption of a smooth, hot/compact stellar wind complies with the observed properties of H-free WN stars. This type of star has the advantage that it is well described by pure He models that contain trace elements beyond He in solar-like abundances. Even if their cores are enriched by the products of He-burning, the mean molecular weight of these stars stays almost the same. This means that their masses and in particular their L/M ratios are well determined by a given (observed) luminosity (Langer 1989).

In the following, we construct models for the winds of hot/compact He stars, assuming that they are optically thick and radiation-driven. To this purpose, we combine the critical conditions arising at the sonic point with additional conditions imposed by the outer wind, assuming a fixed outer wind structure. In Sect. 2 we explain our method and discuss the dependence of our results on the adopted input parameters. In Sect. 3 we apply our models to a comprehensive sample of H-free WN stars in the Galaxy and the Large Magellanic Cloud (LMC). In Sect. 4 we discuss the results.

2. Optically thick winds

As indicated in the previous section, radiation-driven winds can be divided into optically thin and optically thick winds, depending on the location of their critical point. For the thin winds of OB stars, Castor et al. (1975, CAK) formulated a closed theoretical description that allows predicting wind parameters, such as the mass-loss rate Ṁ and terminal wind velocity ν_∞, for a given set of stellar parameters, such as the stellar mass M, luminosity L, and radius R. Within the CAK theory, the wind parameters are determined using so-called force multipliers, that is, parameters that characterise the strength and the distribution of the spectral lines that drive the wind. Line driving is particularly efficient in optically thin winds because the Doppler shifts that arise in the accelerated wind expose spectral lines to an unattenuated radiation field, and reduce the effects of line self-shadowing.

In optically thin winds, these Doppler shifts introduce a dependency of the radiative line acceleration on the velocity gradient ∂ν /∂r. This term alters the equation of motion in a way that a critical point arises at the velocity of a radiative-acoustic wave mode that is faster than the sonic speed (the so-called Abbott speed, cf. Abbott 1980). Moreover, the dependence on ∂ν /∂r enables a simultaneous determination of the mass-loss rate and wind velocity from the conditions arising at the critical point.

When winds become optically thick, the Doppler shifts become less efficient because the photon mean-free path decreases. Whether Doppler shifts are important depends on the CAK optical-depth parameter (Castor et al. 1975) (1)which denotes the ratio between the Sobolev length (ν_Dop/ (∂ν /∂r)) and the photon mean-free path (1 /κ_eρ). Here κ_e denotes the reference continuum opacity due to free electrons, and ν_Dop the line-broadening velocity due to thermal and turbulent gas motions. For WR stars, Nugis & Lamers (2002) estimated values up to t_CAK ≈ 50–150 at the sonic point, which would suggest that Doppler shifts may indeed not affect the sonic-point conditions (but see our discussion in Sect. 4.1). In the remainder of this work, we assume that the wind optical depth is large enough for Doppler shifts to be neglected.

In the optically thick limit, the radiation field is better described by photon diffusion. In this limit, photons tend to avoid strong lines and choose the path of lowest resistance through gaps between lines in frequency space. For this reason, it is much more difficult to launch and accelerate optically thick winds, such as the winds of WR stars or LBVs.

In the diffusion limit, and assuming atomic-level populations in local thermodynamic equilibrium (LTE), the radiative acceleration g_rad can be expressed much more easily by the product of the Rosseland-mean opacity κ_R and the radiative flux F divided by the speed of light c, (2)In this equation, the first equality denotes the general case with the flux-mean opacity κ_F, and the second equality the diffusion limit where κ_F can be expressed as a function of density and temperature via the Rosseland-mean opacity κ_R. In the following, we discuss the sonic-point conditions arising in this limit in Sect. 2.1 and the mass-loss relations that follow from a connection between stellar wind and envelope at the sonic point in Sect. 2.2.

2.1. Sonic-point conditions

The sonic-point conditions for radiation-driven winds in the diffusion limit have been discussed in detail by Nugis & Lamers (2002), Ro & Matzner (2016). Most importantly, the equation of motion has a critical point at the sonic radius R_s where the velocity ν(r) equals the isothermal sound speed a.

According to Nugis & Lamers (2002), this critical condition can be expressed in very good approximation as (3)where ρ_s and T_s denote the density and temperature at the sonic point. This means that the wind solution has to cross the Eddington limit very close to the sonic radius R_s. Furthermore, κ_R needs to increase with radius to ensure that the wind can be accelerated through the sonic point.

For their wind models, Nugis & Lamers (2002) adopted Γ(R_s) = 1, or equivalently, κ_R(ρ_s,T_s) = κ_Edd with the Eddington opacity κ_Edd = 4πcGM/L. This means that in principle, the sonic-point conditions can be extracted directly from opacity tables, where κ_R is given as a function of ρ and T.

Based on dynamical wind models, Ro & Matzner (2016) have shown that the subsonic density and temperature structure of their WR models almost precisely matches the structure of hydrostatic models from Gräfener et al. (2012). For the parameter range that we are interested in, these models also approach Γ = 1 near the Fe-opacity peak, meaning that dynamical terms do not significantly affect the subsonic structure. Therefore we take advantage of this situation and use the density and temperature from hydrostatic models to determine the sonic-point conditions instead of employing more complex hydrodynamical models such as those by Petrovic et al. (2006), Ro & Matzner (2016).

We used the gas pressure P_gas and radiation pressure P_rad as independent variables (instead of density and temperature) for our analysis to describe the conditions at the sonic point. As by definition, ν= a at the sonic point and thus we can express P_gas and P_rad at the sonic point for a given mass-loss rate Ṁ and sonic radius R_s, by (4)where we used the definition a² = ℛT/μ with the gas constant ℛ and mean molecular weight μ, and (5)with the Stefan-Boltzmann constant σ and speed of light c.

This leads to an expression for P_gas as a function of Ṁ, R_s, and P_rad at the sonic point (6)or (7)with P_rad and P_gas in cgs units and μ = 4/3 for an ionised He plasma.

Fig. 1 Critical conditions in Prad–Pgas plane. The hydrostatic envelope solution for a 20 M⊙ He star from Gräfener et al. (2012) is indicated by the black/red curve (the axes in this diagram are inverted so that the stellar interior is located in the lower left and the stellar surface in the upper right). On this curve, κ(Prad,Pgas) ≈ κEdd. Areas in the Prad–Pgas plane with higher/lower opacities are indicated in grey/white and by plusandminus signs. The blue line indicates possible locations of the sonic point for an adopted mass-loss rate of 10-4.8M⊙ yr-1 following Eq. (7). The circle indicates the innermost intersection point with an outward-increasing opacity, giving a precise estimate for the location of the sonic point (see text).

We illustrate in Fig. 1 for a 20 M_⊙ He star with solar metallicity (Z = 0.02) and a given mass-loss rate, how the intersection point between the hydrostatic envelope solution and Eq. (7) determines the values of P_rad and P_gas that a dynamical model would have at the sonic point. The envelope models are identical with the stellar structure models for He stars from Gräfener et al. (2012), which are based on OPAL Rosseland-mean opacity tables (Iglesias & Rogers 1996) with solar-like trace element abundances (Grevesse & Noels 1993).

The envelope model in Fig. 1 follows an S-shaped solution with high P_rad and P_gas in the stellar interior and low P_rad and P_gas at the stellar surface. As pointed out before, the stellar envelope is dominated by radiation pressure and very closely follows a contour with Γ = 1 in the P_rad–P_gas plane as long as the ratio P_rad/P_gas is high. The S-shape of the envelope solution is a consequence of the topology of the opacity in the P_rad–P_gas plane, and in particular of the presence of the Fe opacity peak near log (P_rad) ≈ 6 in cgs units.

The envelope solution and Eq. (7) in Fig. 1 have three intersection points that could in principle become the sonic point of an optically thick wind. However, only for the inner- and outermost intersections does the opacity increase towards the stellar surface. These two points represent the hot and cool wind solutions discussed by Nugis & Lamers (2002).

In the following, we concentrate on the hot solution branch, which is indicated in red in Fig. 1. For these solutions the stellar envelope is compact and the sonic point lies at hot temperatures, below the Fe opacity peak with P_rad ≫ P_gas.

2.2. Connection of wind and envelope

To constrain the mass-loss rates of WR stars, we further take the sonic-point conditions into account that are imposed by the outer wind. This connection of wind and envelope has been discussed in a qualitative manner by Gräfener & Vink (2013).

The underlying idea of their approach was to estimate the temperature at the sonic point of a radiation-driven wind based on a prescribed velocity structure ν(r) assuming that it is radiation-driven, that is, using the opacities that are needed to accelerate the wind. Based on numerical computations, they found a relation between P_gas and P_rad (or alternatively ρ and T) at the sonic point that is imposed by the back-warming effect of the outer wind, namely that P_rad/P_gas ≈ const. In Sect. 2.2.1 we describe our approach to compute this new “wind condition”.

In Sect. 2.2.2 we elaborate on an approximate approach to estimate the wind condition analytically. While the applicability of this result may be limited, the derivations in this section help to understand why the wind condition results in relations with an almost constant ratio of P_rad/P_gas ≈ ν_∞/a_s. We further combine the approximate wind condition with the critical condition (Eq. (3)) to estimate mass-loss rates for WR stars.

In Sect. 2.2.3 we evaluate the wind condition from Sect. 2.2.1 numerically and derive improved mass-loss estimates. We further discuss how the numerical relations for P_rad vs. P_gas relate to the results from Sect. 2.2.2. Finally, we discuss in Sect. 2.2.4 the dependence of the derived mass-loss rates on the dominant input parameter of our models, the terminal wind velocity ν_∞.

2.2.1. P_rad and P_gas at the sonic point

For a given outer wind structure, the value of P_rad, or equivalently the temperature at the sonic point, is a consequence of the back-warming effect of the optically thick wind. To compute P_rad, we used the same formalism as Gräfener & Vink (2013, cf. their Sects. 3.1 and 3.2 for more details,) which is based on the original works by Lucy (1971, 1976), Lucy & Abbott (1993).

P_rad is connected to the radiative acceleration g_rad (Eq. (2)) by(8)This equation can be integrated if an initial value P_ref ≡ P_rad(R_ref) at a given reference radius R_ref, the flux-mean opacity κ_F(r), and the radiative luminosity L_rad(r) are known. Our final models in Sect. 2.2.3 are based on a direct numerical integration of this equation, taking into account that L_rad is a function of radius (cf. Eq. (30)).

To describe the solution of Eq. (8) analytically, we made the simplifying assumption that L_rad = L = const. throughout the wind. As demonstrated by Gräfener & Vink (2013), integrating Eq. (8) from R_ref to the sonic radius R_s leads to the expression (9)for the radiation pressure at the sonic radius, where P_ref is the radiation pressure at the reference radius R_ref, and a weighted mean optical depth of the form (10)This expression for is designed to be comparable to the classical flux-mean optical depth τ_s (cf. Eq. (21)). Because the integral in Eq. (10) is weighted towards small radii with r ≳ R_s, we have .

To compute P_ref and R_ref, we followed the approach of Lucy & Abbott (1993), Lucy (1976) and assumed that P_ref is connected to the radiative flux F_ref at R_ref by an effective (flux) temperature T_ref, that is, (11)where R_ref is the radius where . This relation is to some extent analogous to the concept of the effective temperature in plane-parallel atmospheres.

With Eq. (11), Eq. (9) becomes (12)We note that the only approximations in Eq. (12) are our assumptions that L = L_rad, and that P_ref is given by Eq. (11). For WR stars L = L_rad is usually fulfilled within a margin of several percent, up to ≈10% for the most extreme objects with the strongest winds (e.g. Crowther et al. 2002). The validity of Eq. (11) has been confirmed with a high level of accuracy by Gräfener & Vink (2013) through direct comparison with a dynamical non-LTE atmosphere model for an early-type WC star from Gräfener & Hamann (2005), where the temperature structure is computed from the requirement of radiative equilibrium (cf. Gräfener et al. 2002; Hamann & Gräfener 2003).

In the limit of large wind optical depth, we have and R_ref ≫ R_s, so that Eq. (12) simplifies to (13)where φ_rad(R_s) denotes the radiative momentum flux at the sonic radius. Similarly, Eq. (4) can be expressed as (14)where φ_wind(R_s) denotes the mechanical wind momentum flux at the sonic point. The ratio P_rad/P_gas at the sonic radius can thus be expressed as (15)where Φ_rad(R_s) = L/c and Φ_wind(R_s) = Ṁa denote the total radiative and wind momentum at R_s.

2.2.2. Approximate relation for P_rad/P_gas at the sonic point

An analytical expression for the wind condition in Sect. 2.2.1 can be derived based on Eq. (15), using a well-known relation between the wind-efficiency factor η = Ṁν_∞/ (L/c) (the ratio between the wind momentum and the momentum of the radiation field) and the sonic-point optical depth τ_s, namely that η ≈ τ_s (cf. Netzer & Elitzur 1993; Gayley et al. 1995; Lamers & Cassinelli 1999; Vink & Gräfener 2012).

This relation can be understood under the assumption that photons are scattered multiple times while escaping the stellar wind, undergoing a random walk. For a random walk, the effective number of scatterings that a typical photon must undergo to transfer the necessary momentum to the wind is η². From the point of view of radiative transfer, the number of scatterings in such a picture would be τ². As a consequence, η ≈ τ. Because the layers below the sonic point are near hydrostatic equilibrium, the relevant optical depth for the wind acceleration above this point is the sonic-point optical depth τ_s.

To derive the resulting wind condition, we start with the equation of motion for a spherically expanding, radiation-driven wind with a velocity field ν(r). Above the sonic point, this equation is well described by (16)where gas-pressure terms are neglected. In this equation, the wind acceleration ν× ν′ on the left-hand side is balanced by the outward-directed radiative acceleration and the inward-directed gravitational attraction on the right-hand side. With the Eddington factor Γ = κ_FL_rad/ (4πcGM), defined as the ratio between radiative and gravitational acceleration, this equation can be rewritten as (17)When we integrate this equation from the sonic point (with r = R_s, ν= a_s) towards infinity (i.e. r → ∞, ν→ ν_∞), the value of Γ will vary as a function of r, starting with Γ ≃ 1 at r = R_s, and reaching a maximum near the point of the highest wind acceleration dν/dr. The value of Γ in this region will dominate the integral. For the qualitative considerations in the remainder of this section, we therefore adopt Γ = Γ_w = const., where Γ_w is meant to be representative for the value of Γ in the region of highest wind acceleration. Integration of Eq. (17) then results in (18)with the escape velocity defined as . With ν_∞ ≫ a_s , this simplifies to (19)or equivalently, (20)In a second step, Eq. (17) can be rewritten using the equation of continuity Ṁ = 4πρνr²(21)with the flux-mean optical depth dτ = −κ_Fρ dr. Again, integrating from the sonic point towards infinity and adopting Γ = Γ_w = const., we obtain (22)with the sonic-point optical depth τ_s. Finally, we obtain (23)for the wind efficiency factor η, where we made use of Eq. (20) and ν_∞ ≫ a_s. Analogous to Eq. (15) Φ_rad(∞) = L/c and Φ_wind(∞) = Ṁν_∞ are the total radiative and wind momenta for r → ∞.

Given that usually ν_∞ > ν_esc for hot-star winds, this means that the wind efficiency η is of the same order of magnitude as τ_s with a correction factor that reflects the fraction of the wind momentum that is used to overcome the gravitational attraction of the star. We note that this relation is only approximate, primarily because of our assumption that Γ = Γ_w = const., but it will likely provide a good qualitative picture of the involved dependencies.

Combining Eqs. (15) and (23) with (cf. Eq. (10)), we finally obtain (24)or (25)for the ratio P_rad/P_gas at the sonic point. This means that the ratio of P_rad/P_gas is of the same order of magnitude as ν_∞/a_s and thus almost a constant, in line with the numerical results of Gräfener & Vink (2013).

Equation (25) is independent of Ṁ and can be used to estimate Ṁ as a function of ν_∞. To this purpose, we express the sonic speed a_s through the radiation pressure P_rad by (26)so that Eq. (25) can be written as (27)and with μ = 4/3(28)This relation provides an estimate of the sonic-point conditions that are imposed by an optically thick wind with a specified terminal wind velocity ν_∞ in the P_rad–P_gas plane. ν_esc and μ are given by the stellar mass, radius, and chemical composition.

Fig. 2 Mass-loss determination in Prad–Pgas plane. The plot is analogous to Fig. 1, but now the black dashed line indicates the sonic-point conditions of an optically thick wind with an adopted terminal wind velocity of ν∞ = 1800 km s-1 following Eq. (27). Again, the circle indicates the innermost intersection point of Eq. (25) with the static envelope solution. Given the values of Prad and Pgas at the intersection point, Ṁ can be determined from Eq. (7).

The combination of the wind condition (Eq. (28)) with the critical condition Γ ≃ 1 at the sonic point (Eq. (3)) provides an estimate of the mass-loss rate of an optically thick wind. This is illustrated in Fig. 2, where the intersection point between Eq. (28) and the hydrostatic envelope solution (with Γ ≃ 1) specifies the values of P_rad and P_gas for which envelope and wind can be connected at the sonic point. The mass-loss rate that corresponds to these values can then be computed from Eq. (7) in the same way as illustrated in Fig. 1. For the example in Fig. 2, we obtain a value of Ṁ = 10^-4.72M_⊙ yr^-1 for an adopted terminal wind velocity of ν_∞ = 1800 km s^-1, which lies within the observed range of WR mass-loss rates.

Notably, the denominator in the last term of Eqs. (27) and (28) has a minimum for ν_∞ = ν_esc. This means that P_gas, and thus also Ṁ, are reaching a maximum for ν_∞ = ν_esc. If the approximations made in this section were strictly valid, this maximum value would introduce new upper limits for the mass-loss rates of WR stars and LBVs.

It is further important to note that Eqs. (25) to (28) do not depend on the uncertain clumping properties of the wind material, as the underlying Eqs. (17) and (21) only rely on the assumption that the winds are radiation-driven. In case of an inhomogeneous wind structure, τ_s would thus represent an integral over an effective mean opacity, as discussed for instance by Oskinova et al. (2007), Owocki (2008), Sundqvist et al. (2014). An analysis involving the ratio P_rad/P_gas, such as that performed by Gräfener & Vink (2013), thus provides a new important means to diagnose the clumping properties of the wind material in deep layers that are not directly observable.

2.2.3. Numerical wind models

In the following, we use numerical wind models to derive sonic-point conditions and mass-loss rates in an analogous way as in Sect. 2.2.2. To this purpose, we followed the approach of Gräfener & Vink (2013) and solved Eq. (8) numerically, in line with the weighted mean optical depth and reference radius R_ref. To do this, we adopted a β-type wind velocity law of the form (29)with β = 1 and a prescribed value of ν_∞, and compute κ_F(r) under the assumption of radiative driving, that is, in line with the wind acceleration ν× ν′ resulting from the prescribed velocity law (cf. Eq. (16)). Furthermore, we computed L_rad(r) by taking into account that the local radiative luminosity in a radiation-driven wind is reduced by the mechanical wind luminosity, that is, (30)Here we neglected enthalpy terms, which we found to be negligible. The inclusion of L_wind is important for cases where L_wind → L. In such cases, L_rad is significantly reduced by the mechanical work that is needed to accelerate the wind, and radiative wind driving becomes increasingly difficult. No wind acceleration is possible beyond the so-called photon-tiring limit with L = L_wind (Owocki & Gayley 1997; cf. also Heger & Langer 1996; Gräfener et al. 1998; Owocki et al. 2004).

The use of a prescribed velocity law means that we did not solve the equation of motion explicitly. Instead, we derived the sonic-point conditions under the assumption that the winds are radiation-driven. Because the sonic-point values of P_rad and P_gas are the result of the integrated radiative acceleration throughout the outer wind, our results depend predominantly on the adopted terminal wind velocity ν_∞, and to a lesser extent on the detailed velocity structure. Gräfener & Vink (2013) performed a detailed comparison of the sonic-point conditions arising from a β-type velocity law as we used here, and the dynamically consistent WC star model from Gräfener & Hamann (2005). Although the model of Gräfener & Hamann displayed a qualitatively different velocity structure with two acceleration zones, the resulting sonic-point conditions were remarkably similar, with differences of 0.06 dex in P_rad and 0.01 dex in P_gas.

The mass-loss rates based on our numerical models were obtained in analogy to Sect. 2.2.2, that is, we computed a static stellar structure model for a given stellar mass M and used the resulting stellar luminosity L and radius R as an input for our wind models. As stellar envelopes often inflate near the region of the Fe-opacity peak and we are interested in compact solutions where the sonic point is located just below the inflated region, we did not extract the formal surface radius from our models, but used the core radius R_c at the bottom of the inflated layer instead (for more details see Gräfener et al. 2012). R_c serves as an inner boundary for our wind models and is usually located very close to the sonic radius R_s.

The precise value of R_s can only be determined if Ṁ is known. In our final models in Sect. 3.2, we resolve this problem through an iterative procedure. In the present section, we accept a small inconsistency (up to ~5% in R and 10% in Ṁ), which arises because R_c deviates slightly from R_s in Eq. (7).

To derive the sonic-point conditions imposed by the outer wind, we computed a series of wind models (typically for ten different values of Ṁ) with the stellar parameters from above, and a given terminal wind velocity ν_∞. The resulting sonic-point values of P_rad and P_gas are indicated by the blue curve in Fig. 3. The mass-loss rate is computed from the sonic-point values of P_rad and P_gas at the intersection point between the envelope and wind solutions, using Eq. (7). In our example, the resulting mass-loss rate of Ṁ = 10^-4.65M_⊙ yr^-1 compares very well with the approximate value derived in Sect. 2.2.2.

Fig. 3 Mass-loss determination in Prad–Pgas plane. The plot is analogous to Fig. 2, but now the sonic-point conditions following from our numerical models (indicated by the solid blue curve) are used to determine Ṁ.

Fig. 4 Sonic-point ratio of Prad/Pgas. Different contributions to the ratio of Prad/Pgas at the sonic point, as imposed by the outer wind, are indicated on a logarithmic scale. The left panel shows the same sequence of wind models as in Fig. 3, i.e., models with a fixed terminal wind velocity of ν∞ = 1800 km s-1 as a function of the adopted mass-loss rate Ṁ. The right panel shows the same model sequence as in Fig. 5, i.e., models for which Ṁ is derived consistently with the critical condition (Eq. (3)) as a function of ν∞. We plot the ratio of Prad/Pgas that follows from our numerical models (solid blue line), from the approximate relation Eq. (25) (black dashed line), and Eq. (31) (red circles). The grey dotted and dash-dotted lines indicate the (weighted) sonic-point optical depths τs and respectively. The red solid line indicates the ratio τs/η, and the red dashed line the factor . The black solid line indicates the offset factor in Eq. (31).

Figure 3 demonstrates that there is a very good agreement between our approximate wind relation (Eqs. (27) and (28), indicated by the black dashed line) and our numerical wind models (blue) in the region near the actual sonic point. Deviations occur for very high wind densities, when L_wind becomes comparable to L, and for low wind densities, when the wind becomes optically thin.

To better understand the deviations from Eqs. (27) and (28), we investigate in Fig. 4 (left panel) the different contributions to the ratio P_rad/P_gas for the model sequence from Fig. 3. To this purpose, we rewrote Eq. (12) using the definition of the wind efficiency η = Ṁν_∞/ (L/c), and divided it by Eq. (4), which led to(31)In Fig. 4 we extract the different quantities on the right-hand side (rhs) of this equation from our numerical models. The approximations that we made when deriving Eqs. (25)–(28) were that the second term on the rhs of Eq. (31) is well described by Eq. (23), and that terms three and four are of order unity. We note that we expect Eq. (31) to represent our numerical results very accurately, except for the fact that we assumed L_rad = const. instead of taking Eq. (30) into account.

Our numerical results in Fig. 4 (indicated by the blue solid line) deviate from Eq. (25) (black dashed line) for low and high Ṁ. For the case of low Ṁ, it is clear that the optical depths τ_s and (indicated in grey) become too small to justify our assumption of large optical depth. As a consequence, the last term in Eq. (31) (black solid line) exceeds unity for low Ṁ. For high the situation is different because the sum of the terms on the rhs of Eq. (31) (indicated by red circles) displays a very good agreement with Eq. (25). Our numerical results deviate from this value for high Ṁ because the mechanical wind luminosity L_wind becomes comparable to L in Eq. (30). As a consequence, L_rad is reduced in the outer wind, leading to a weaker back-warming effect, and thus lower P_rad.

Furthermore, Fig. 4 demonstrates that the ratio on the rhs of Eq. (31) deviates from unity, which is expected given the definition of (cf. Eq. (10) and the related comments). In our example, this deviation is canceled by the fact that the second term τ_s/η (red solid line) is larger than what is expected from Eq. (23) and exceeds (red dashed line) by a similar amount.

Altogether, our analysis in Fig. 4 indicates that Eq. (25) may indeed be useful to understand the qualitative dependence of the sonic-point conditions on the outer wind, but its application in quantitative models is probably hampered by a too low accuracy. Instead, it is necessary to employ numerical models for quantitative predictions.

2.2.4. Mass-loss rates as a function of ν_∞

In Sect. 2.2.3 we established a new method to compute the mass-loss rates of hot stars with optically thick winds with a given terminal wind velocity ν_∞. For practical applications, it will be possible to estimate ν_∞, for instance, for a certain class of objects, either observationally (based on spectroscopic results) or from common theoretical arguments (e.g. that ν_∞ should be related to ν_esc). In any case, it will be of primary interest how sensitively our results depend on the choice of ν_∞.

For this reason, we performed the same computations as described in Sects. 2.2.2 and 2.2.3 for a 20 M_⊙ He star model for various values of ν_∞. The results are shown in Fig. 5. As expected, the mass-loss rates obtained from our approximate relation Eq. (27) (black dashed curve in Fig. 5) show a maximum for ν_∞ = ν_esc. This maximum does not occur for our numerical models (blue curve). The reason for this discrepancy is most likely our simplifying assumption that Γ = Γ_w = const. in Eqs. (20) and (23), leading to the given dependence of Eq. (25). Nevertheless, for ν_∞ < ν_esc the qualitative behaviour of Ṁ vs. ν_∞ is well described by our approximate relation. Near the maximum of Ṁ at ν_∞ = ν_esc , the numerical results also start to saturate, and the dependence on ν_∞ becomes relatively weak. For the lowest values of ν_∞ , the discrepancy increases again, most likely because our assumption that ν_∞ ≫ a does not hold anymore.

In the right panel of Fig. 4 we show the different contributions to the ratio P_rad/P_gas as a function of ν_∞. First of all, the figure shows that τ_s as well as are large enough to justify our assumption of an optically thick wind. As a consequence, the last term in Eq. (31) is indeed very close to one. As in our example from Sect. 2.2.3, the second and third term in Eq. (31) deviate from unity. However, in the present case, the respective errors do not compensate each other. In particular, the difference between τ_s and increases for large ν_∞, leading to a discrepancy between our approximate and numerical results described above.

To summarise, our approximate relation Eq. (25) is not applicable in all cases. However, the flattening of the relation Ṁ(ν_∞) near the maximum at ν_∞ = ν_esc in Eq. (25) also leads to a flattening of this relation in our numerical models, and thus to a weak dependence of Ṁ on ν_∞. This means that our numerical mass-loss estimates should be reliable as long as reasonable estimates of ν_∞ can be provided.

3. Mass-loss relations for Wolf-Rayet stars

In this section we use our method from Sect. 2 to derive mass-loss relations for H-free WN stars. To this purpose, we investigate the terminal wind velocities of these objects in the Galaxy and the LMC in Sect. 3.1, and compare the resulting mass-loss relations with empirical results in Sect. 3.2.

3.1. Terminal wind velocities of H-free WN stars

In this section we investigate the terminal wind velocities ν_∞ of H-free WN stars from previous empirical studies of WR stars in the Galaxy and LMC. The main purpose is to determine which values of ν_∞ are appropriate as input parameters for our models. An important question in this context is the dependence of ν_∞ on stellar parameters.

Fig. 5 Mass-loss predictions. The derived mass-loss rates for a 20 M⊙ He star are indicated as a function of the adopted terminal wind velocity ν∞. The solid blue curve indicates numerical wind models (cf. Sect. 2.2.3) and the dashed curve our approximate estimates (Sect. 2.2.1).

A common way to implement such a dependency is to assume a constant ratio of ν_∞/ ν_esc. The idea behind this approach is a (near) equipartition between the gravitational and kinetic contributions to L_wind in Eq. (30). Such an equipartition is predicted in the theory of optically thin, radiation-driven winds of Castor et al. (1975), and would also follow from our Eq. (19) assuming that WR stars of similar spectral type display similar Γ_w. Empirically, relations of this form have been confirmed for the winds of OB stars, which tend to show constant ratios of ν_∞/ ν_esc well above one (e.g. Lamers et al. 1995). In addition, for the optically thick winds of WC stars, Gräfener & Vink (2013) found a relatively clear relation with ν_∞/ ν_esc ≈ 1.6.

In Fig. 6 we show the terminal wind velocities ν_∞ for putatively single, H-free WN stars in the Galaxy and LMC as a function of ν_esc. The data are compiled from the comprehensive studies of WN stars in the Galaxy by Hamann et al. (2006) and the LMC by Hainich et al. (2014). While Hamann et al. concentrated on single WN stars in their work, the study of Hainich et al. also included binaries, and we removed all stars with indications for multiplicity from the sample. The main parameters of our sample stars are compiled in Tables 1 and 2.

Fig. 6 Terminal wind velocities ν∞ from Hamann et al. (2006), Hainich et al. (2014). Empirical values for single, H-free WN stars in the Galaxy (blue) and LMC (grey) are indicated as a function of the escape velocity νesc. Dashed lines indicate relations with ν∞/ νesc = 1.3 and ν∞ = 1800 km s-1.

The ν_esc in Fig. 6 are estimated on the basis of the spectroscopically determined radii R_⋆ resulting from the analyses of Hamann et al. (2006), Hainich et al. (2014), and stellar masses M obtained from the mass-luminosity relation for pure He-stars from Gräfener et al. (2011). The radii R_⋆ denote the inner boundary of the atmosphere models used to analyse the stars. They are located at large optical depth (typically of the order of τ = 20) in the nearly hydrostatic layers of the atmosphere models and are almost identical to the sonic radius R_s.

It is important to note that because of the WR radius problem, the spectroscopic radii used here do not necessarily resemble the radii that we compute in our present (compact) stellar structure models (cf. also the discussion in Sect. 4.3). As a consequence, only part of the stars in the observed sample may be representative for our models. Our present models display escape velocities in the range ν_esc ≈ 2200–2600 km s^-1, meaning that they hardly vary.

In Fig. 6 only the coolest stars, with ν_esc ≲ 1200 km s^-1, seem to follow a relation with ν_∞/ ν_esc ≈ 1.3. This trend vanishes at higher escape velocities and the terminal wind velocities seem to saturate. For the Galactic stars, this saturation seems to occur near 1800 km s^-1 and only few objects lie in the saturated regime. For the LMC stars, almost all H-free WN stars seem to lie in the saturated regime. According to Hainich et al. (2014), their spectra are in agreement with a constant value of ν_∞ = 1600 ± 200 km s^-1, with only two exceptions.

For comparison, we extracted terminal wind velocities from the morphological ananlyses of UV line-profiles of Niedzielski & Skorzynski (2002), Niedzielski et al. (2004), which are based on spectroscopic data taken with the International Ultraviolet Explorer (IUE) satellite. The results for the sub-sample overlapping with our sample are plotted in Fig. 7 (cf. also Tables 1 and 2). There are significant discrepancies with respect to the values derived by Hamann et al. (2006), Hainich et al. (2014), and the data of Niedzielski & Skorzynski (2002), Niedzielski et al. (2004) have a tendency to show higher values of ν_∞. In Fig. 7 now the complete Galactic sample seems to follow a relation with ν_∞/ ν_esc ≈ 1.3, and the LMC data are, with few exceptions, in much better agreement with this relation. These discrepancies have also been discussed by Hainich et al. (2014), who noted that adopting the higher ν_∞ from Niedzielski et al. leads to a mismatch of the line widths in the optical range.

The reason for the discrepancies is unclear. Notably, Niedzielski et al. (2004) also reported peculiarities and inconsistencies between ν_∞ determined from different ionic species. In particular for the LMC, the quality of the IUE data is partly very poor. Moreover, the pointing accuracy and large aperture size of the IUE may not be sufficient to separate stars in dense fields in the LMC. One could also speculate that the UV profiles do not always stem from the same cool wind component as the optical lines, similar to what has been discussed for clumped, multi-component winds of O stars (Sundqvist et al. 2014).

In the following analysis we acknowledge that the terminal wind velocities of WR stars are uncertain and compute different model sequences with ν_∞ = const. and ν_∞/ ν_esc = const.

Fig. 7 Terminal wind velocities ν∞ from Niedzielski & Skorzynski (2002), Niedzielski et al. (2004). Empirical values for a subset of the stars in Fig. 6 are indicated as a function of the escape velocity νesc. Dashed lines indicate relations with ν∞/ νesc = 1.3 and ν∞ = 1800 km s-1.

3.2. WR mass-loss rates

In this section we use our models from Sect. 2.2.3 in combination with terminal wind velocities from Sect. 3.1 to predict mass-loss rates of H-free WN stars. We computed two model sequences for Galactic metallicity (Z = 0.02), one with a fixed terminal wind velocity of ν_∞ = 1800 km s^-1, and one with ν_∞/ ν_esc = 1.3. For the compact stellar structure models used in this work, the latter relation results in values of ν_∞ = 2800–3400 km s^-1. For the LMC metallicity (Z = 0.008), we computed two sequences with ν_∞ = 1600 km s^-1 and ν_∞/ ν_esc = 1.3.

In contrast to our previous computations, we estimated the sonic radius R_s as precisely as possible by matching the sonic-point density following from Eq. (4) iteratively for each given value of Ṁ with the corresponding density in our static stellar structure models. The resulting mass-loss rates Ṁ are shown as a function of luminosity L in Figs. 8 and 9.

In these two figures, we compare our results with our samples of putatively single H-free WN stars in the Galaxy and LMC from Tables 1 and 2 extracted from Hamann et al. (2006), Hainich et al. (2014). Because different wind clumping factors D were adopted in these two studies, we scaled the mass-loss rates from Hamann et al. (which were originally obtained with D = 4) with to match those from Hainich et al. with D = 10. The same inconsistency has recently been pointed out by Yoon (2017), who obtained a clumping-corrected relation for H-free WN stars by fitting the corrected data from Hamann et al. (2006) using the luminosity dependence with Ṁ ∝ L^1.18 derived by Hainich et al. (2014) for WN stars in the LMC.

In Figs. 8 and 9 we compare our results with the commonly used empirical mass-loss relation from Nugis & Lamers (2000) scaled with the Z-dependence from Vink & de Koter (2005), that is, (32)and the relation from Yoon (2017), that is, (33)The theoretical relations from our Galactic models in Fig. 8 display a nice agreement with relation from Yoon (2017). In particular, the slope of our predicted mass-loss relations (with roughly Ṁ ∝ L^1.3) compares well with the slope proposed by Yoon (Ṁ ∝ L^1.18) and is significantly shallower than the relation from Nugis & Lamers (Ṁ ∝ L^1.63).

Fig. 8 Mass-loss relations for Galactic metallicity. Mass-loss rates for H-free WN stars in the Galaxy as a function of luminosity. Black/grey stars indicate empirical mass-loss rates for Galactic WN stars with/without known distance from cluster membership (Hamann et al. 2006). The empirical mass-loss rates have been scaled down to those expected for a wind clumping factor of D = 10. The red curves indicate mass-loss rates obtained from our models with Z = 0.02 for ν∞/ νesc = 1.3 (top) and ν∞ = 1800 km s-1 (bottom), the labels indicate stellar masses in M⊙. The dotted/dashed lines indicate the empirical mass-loss prescriptions for H-free WN stars from Nugis & Lamers (2000) and Yoon (2017).

Fig. 9 Mass-loss relations for LMC metallicity. Mass-loss rates for H-free WN stars in the LMC as a function of luminosity. Black stars indicate empirical mass-loss rates from (Hainich et al. 2014). The blue curves indicate mass-loss rates obtained from our models with Z = 0.008 for ν∞/ νesc = 1.3 (top) and ν∞ = 1600 km s-1 (bottom), the labels indicate stellar masses in M⊙. The dotted/dashed lines indicate the empirical mass-loss prescriptions for H-free WN stars from Nugis & Lamers (2000) scaled with the metallicity dependence from Vink & de Koter (2005), and from Yoon (2017), respectively. Notably, our models do not predict solutions below ~25M⊙, which contradicts the observed range of luminosities for H-free WN stars in the LMC.

Fig. 10 Limiting masses. Wind and envelope solutions for a 14 M⊙ He star analogous to Fig. 3. For Z = 0.02 the envelope solution (red/black curve) just matches the wind solution with ν∞ = 1800 km s-1 (blue curve). For Z = 0.008 (black dashed curve), no match can be obtained.

Another important result is that we encounter stellar mass limits below which no wind solutions exist. For our Galactic models, this situation occurs for models ≲14 M_⊙. To illustrate the situation, we show the solution topology for our lowest-mass model with 14 M_⊙ in Fig. 10. For this mass, the wind solution with ν_∞ = 1800 km s^-1 just touches the tip of the envelope solution, so that there still exists an intersection point. For lower masses, the L/M ratio would decrease, so that the Eddington opacity κ_Edd = 4πcGM/L, which needs to be matched at the sonic point due to the critical condition (Eq. (3)), would increase. As a consequence, higher densities would be required at the sonic point, and the envelope solution in Fig. 10 would move downward, towards higher values of P_gas. In such a case, there would be no intersection point between envelope and wind solution, which means that a connection of wind and envelope at the sonic point would be physically impossible.

The same situation occurs for lower metallicities. The envelope solution for a 14 M_⊙ model with LMC metallicity is indicated as a dashed curve in Fig. 10. In this case, higher densities are also needed to produce the required value of κ_Edd. Clearly, no intersection point exists between the dashed envelope solution and the wind solution in Fig. 10.

To shift the wind solution accordingly, it would be necessary to reduce the ratio of P_rad/P_gas in our example by a factor of about five. Since P_rad is confined within a small range to match the temperature of the Fe opacity peak, this would mean that P_gas, and thus Ṁ (according to Eq. (4)), would need to increase by a similar factor. Following our discussion in Sect. 3.1, such an increase would demand a substantial increase of ν_∞ , which is not observed. Because the mechanical part of the wind luminosity L_wind in Eq. (30) increases with , such a wind will approach the photon-tiring limit (with L = L_wind, cf. Eq. (30)) before the desired ratio of P_rad/P_gas is reached.

As a consequence of the metallicity-dependence of the mass limits discussed above, our LMC models with Z = 0.008 in Fig. 9 only produce solutions for masses ≳ 25 M_⊙. This limit corresponds to a minimum luminosity of log (L/L ⊙ ) ≳ 5.8, a value above which hardly any stars are observed. Clearly, our models do not reproduce the optically thick winds of the bulk of the H-free WN stars in the LMC.

Our high-mass models with M ≳ 25 M_⊙, suggest a metallicity dependence Ṁ ∝ Z^γ with an exponent of γ = 0.8–1.0, which is comparable to the value of γ = 0.86 from Vink & de Koter (2005) and higher than the value of γ = 0.6 that has been derived by Yoon (2017) based on observations of H-free WN stars from Hamann et al. (2006) and the Ṁ(L) relation for LMC stars from Hainich et al. (2014).

4. Discussion

We investigated whether WR-type stellar winds can be described by the “classical” assumptions of a compact stellar structure and a radiation-driven, optically thick, smooth stellar wind. As it became clear in the previous sections, these classical assumptions can describe many of the observed properties of H-free WN stars well, but do also fail in several respects. In the following, we summarise and discuss our results. In Sect. 4.1 we start with a discussion of the wind physics implied by our models as well as their possible shortcomings, followed by a discussion of their dependence on ν_∞ in Sect. 4.2. In Sect. 4.3 we discuss the implications of the observed terminal wind velocities, and in Sect. 4.4 the relevance of the derived WR mass-loss relations and mass limits. In Sect. 4.5 we compare our results with previous theoretical studies.

4.1. Wind physics

In Sect. 3.2 we provided the first theoretical Ṁ(L) relations for WR stars that are obtained through a combination of the critical condition that occurs at the sonic point of optically thick winds (Eq. (3)), and the sonic-point condition arising from the back-warming effect of the outer wind (the “wind condition”, Eq. (12)). While the critical condition only depends on the local conditions at the sonic point, the wind condition depends on the integrated properties of the outer wind above the sonic point, in particular, on its optical depth. The mass-loss rates obtained from our models are thus a consequence of both the local conditions at the sonic point and the physics of the outer wind.

The temperature structure resulting from Eq. (12) is a consequence of radiative equilibrium. The interplay between the sonic point and the outer wind is thus moderated by radiative heat transfer. This means that information on the state of the outer wind is communicated downstream, towards the sonic point, through radiation.

In our current models we assume that the outer wind is radiation-driven. Notably, we do not explain how the wind acceleration is sustained, especially beyond the Fe-opacity peak. In previous wind models by Gräfener & Hamann (2005), this was achieved by non-LTE effects, namely through the population of meta-stable energy levels that typically occur directly above the ground state of Fe-group ions. Nevertheless, Gräfener & Hamann reported that WR wind driving is still problematic in the general case because of the lack of opacities in the gap between the Fe-opacity peak and the opacity peak at cooler temperatures.

A possible way to increase the flux-mean opacity in stellar winds is the line-deshadowing effect in the presence of velocity gradients. As discussed in Sect. 2, the occurrence of this effect is governed by the CAK optical depth parameter t_CAK (Eq. (1)). A critical parameter in Eq. (1) is the line-broadening velocity ν_Dop. Identifying ν_Dop with the thermal velocity of protons, Nugis & Lamers (2002), Ro & Matzner (2016) argued that Doppler shifts are not important in the sonic-point region of WR-type winds. In the models of Gräfener & Hamann (2005), high values of ν_Dop were also adopted and velocity gradients played a subordinate role.

Given the extremely narrow thermal width of iron lines, however, ν_Dop may be much smaller. This could mean that the effects of velocity gradients have been underestimated in previous studies. In the presence of high line densities, this could be further facilitated by multi-line scattering effects (Friend & Castor 1983; Gayley et al. 1995), potentially leading to a violation of the diffusion approximation even for near-sonic velocities.

In this context, it is crucial whether the line-broadening velocity ν_Dop is dominated by the slow thermal motions of Fe ions or by turbulent motions that may be orders of magnitude higher. While Grassitelli et al. (2016), for example, argued for a turbulent structure in the sub-sonic part of WR envelopes, the low values of ν_Dop may be supported by the suppression of convection in the presence of velocity gradients as discussed by Ro & Matzner (2016). We can therefore not exclude at the current stage that velocity gradients play an important role in WR wind driving, in contrast to the underlying model assumptions in this work, but our present models serve as an important test case for WR wind theories.

4.2. Dependence on wind velocities

An interesting new aspect of our wind models is their dependence on the terminal wind velocity ν_∞, which is used as an input parameter in our models. In nature, ν_∞ is set by the complex physics of the outer wind. There are, however, situations where small changes in stellar parameters lead to drastic changes in the outer wind structure and ν_∞. A well-known example is the bi-stability jump for B-supergiants (Vink et al. 1999; Petrov et al. 2016). Differences in the chemical composition will also be able to cause variations of ν_∞. Because we did not compute the velocity structure self-consistently, such situations can only be considered in our models by varying ν_∞ as input parameter.

In Sect. 2.2.4 we found that Ṁ increases with increasing ν_∞ for our optically thick wind models. This is in stark contrast to what is found for the optically thin winds of OB stars for which the wind momentum Ṁν_∞ is expected to be almost invariant for fixed stellar parameters. As a consequence, the wind models for OB stars by Vink et al. (2000), Vink & de Koter (2002) display the opposite behaviour to ours, namely that Ṁ decreases for increasing ν_∞.

The reason for the peculiar behaviour of optically thick winds is that because of the effect of the outer wind, the ratio of P_rad/P_gas at the sonic point decreases for increasing values of ν_∞ (corresponding to the case with ν_∞ < ν_esc in Eq. (25)). Because P_rad is constrained through the temperature of the Fe-opacity peak, this means that P_gas, and thus the density ρ_s, at the sonic point increases for increasing ν_∞. As a consequence, the mass-loss rate increases for increasing ν_∞ as well. This may be the reason why the observed WR wind velocities saturate easily, leading to terminal wind velocities ν_∞ of similar order or even lower than the escape velocity ν_esc.

4.3. Observed wind velocities and radii

The terminal wind velocity is the key parameter for the parametrisation of the outer wind in our models. For this reason, we reviewed the observed wind properties of H-free WN stars in Sect. 3.1 using stellar parameters and ν_∞ from Hamann et al. (2006), Hainich et al. (2014), and alternative values determined from UV line profiles by Niedzielski & Skorzynski (2002), Niedzielski et al. (2004).

Based on the spectroscopically determined WR radii and ν_∞ from Hamann et al. (2006), Hainich et al. (2014), we found that predominantly late spectral sub-types in the Galaxy seem to follow a relation with ν_∞/ ν_esc = const., similar to the relation for WC stars from Gräfener & Vink (2013). The rest of our sample stars seems to display almost constant values of ν_∞. Using the from Niedzielski & Skorzynski (2002) and Niedzielski et al. (2004), all Galactic objects and, with few exceptions, also the LMC objects seem to follow a relation with ν_∞/ ν_esc = const.

The existence of such relations is relevant for the WR radius problem, that is, for the question whether the radii of WR stars are compact, as adopted in our present models, or if they are much larger, as suggested by spectroscopic analyses. Whether the spectroscopic radii are reliable depends on the wind density. Hamann & Gräfener (2004) discussed the limit where the visible spectrum is formed so far out in the wind that it is not possible to determine the radius of the wind base. Sander et al. (2012), Hainich et al. (2014) discussed that this limit predominantly affects the earliest WR sub-types, and that the spectroscopic radii of most late-type WR stars are most likely well determined.

Based on our discussion in Sect. 3.1, the existence of a relation with ν_∞/ ν_esc = const. for WC stars and at least part of the WN stars suggests that their wind base may indeed be located at larger radii. On the other hand, a constant value of ν_∞ would suggest that ν_esc is also (almost) constant, as is expected for compact WR stars. Because we cannot clearly distinguish between the two scenarios based on the current data, it is currently not clear whether WR stars are generally compact. However, it is important to keep in mind that our assumption of a compact core may not be correct for all stars in the sample.

For our final model computations we used both approaches, that is, we used prescriptions with ν_∞/ ν_esc = const. and ν_∞ = const. Because the first relation applies to late sub-types with large radii, the resulting values for our compact models (ν_∞ = 2800–3400 km s^-1) are larger than what is typically observed for this type of objects. Nevertheless, it is important to investigate how different assumptions on ν_∞ affect our results.

4.4. Mass-loss rates and mass limits

In Sect. 3.2 we compared the Ṁ(L) relations resulting from our models with observations. We found a good agreement with the new empirical mass-loss relation for H-free WN stars of Yoon (2017). Empirical mass-loss relations for WR stars are commonly characterised by power laws with fixed exponents for the dependence on luminosity and metallicity. Our models display such a relation with Ṁ ∝ L^δZ^γ, δ = 1.3, and γ = 0.8–1.0. The exponents in this relation do not depend significantly on the adopted prescription for ν_∞.

In an absolute sense, our predictions are in good agreement with the mass-loss prescription of Yoon (2017). However, the results depend upon the adopted prescription for ν_∞, in the sense that the models with constant (and lower) ν_∞ predict slightly lower mass-loss rates (cf. Sect. 4.2).

The most notable difference between the observed WR population and our models are the Z-dependent mass limits below which our models do not provide wind solutions. In particular for the low metallicity in the LMC, we find that the Fe peak does not provide enough opacity to support WR-type mass-loss rates in the region around the sonic point.

A possible solution to this problem has been discussed previously by Gräfener & Vink (2013). These authors suggested that clumping could increase the mean opacity at the sonic point. According to Gräfener et al. (2012), clumping could also lead to an enhanced inflation effect and would thus help to explain the large radii that are derived in spectral analyses of WR stars. However, as we discussed earlier, such a configuration cannot be realised with a smooth and continuous wind flow (cf. Ro & Matzner 2016) and would most likely require a very complex wind structure. An alternative way to increase the flux-mean opacity near the sonic point could be the line-deshadowing and multi-line scattering effects discussed in Sect. 4.1.

In this context, it is interesting that according to Sander et al. (2012), the population of WC stars in the Galaxy shows a clear dichotomy between early and late WC sub-types, where the late sub-types occur predominantly at low luminosities ≲ 10^5.3L_⊙. Sander et al. argued that the large radii inferred for the latter group suggest that they are inflated. According to our models, there are indeed no continuous wind solutions for smooth/compact WR-type winds in this luminosity regime, so a clumped/inflated wind structure may be preferred.

As discussed in Sect. 3.2, the encountered mass limits result from the critical condition (Eq. (3)) that the Eddington opacity κ_Edd = 4πcGM/L needs to be matched at the sonic point. For a given chemical surface composition, the truly decisive parameter is thus the L/M ratio rather than the stellar mass. Only for objects with a given L(M) relation, as is the case for the core He-burning objects investigated here, the encountered limits translate into limits for the stellar mass M, or alternatively for the luminosity L. In particular, the derived mass limits do not apply to WR-type central stars of planetary nebulae (PNe), which are in a phase of shell burning. Strictly speaking, our results only apply to pure He stars. In reality, the cores of He-burning stars become enriched in C and O during their evolution. This leads to a slight increase in mean molecular weight, and consequently, to an increase of the L/M ratio. However, this effect is not strong enough to resolve the discrepancies that we encounter at LMC metallicity.

Even if the encountered mass limits need to be modified, the existence of physical mass limits per se is interesting. The WR stars in the Galaxy and LMC, including the samples discussed here, do display lower luminosity/mass limits. In the framework of single-star evolution, such mass limits are expected (although current evolution models struggle to reproduce the luminosities of the observed WR population, e.g. Sander et al. 2012). On the other hand, binary-interaction models predict the existence of a continuous mass distribution of WR stars in binary systems, including low-mass WR stars, which are believed to be supernova (SN) Ib/c progenitors (e.g. Smith et al. 2017; Zapartas et al. 2017). In this framework, the existence of mass limits for optically thick winds could explain why such low-mass WR stars have not yet been observed. He-stars below the corresponding WR luminosity limits would not be able to support optically thick winds, that is to say, they would not display the prototypical WR emission lines.

4.5. Comparison with previous wind models

Previous models for core He-burning WR stars have been computed based on Monte Carlo models (Lucy & Abbott 1993; Vink & de Koter 2005) and based on hydrodynamic non-LTE model atmospheres (Gräfener & Hamann 2005).

The Monte Carlo models focused only on the outer wind using β-type velocity laws with prescribed ν_∞. Mass-loss rates Ṁ were determined through the condition that the global wind luminosity (Eq. (30)) equals the work performed by the radiative force resulting from the Mote Carlo computations. This means that no local dynamical consistency was obtained, and in particular, that the critical and wind conditions at the sonic point were ignored. Our present models are complementary to this approach in the sense that we assume that the outer wind is radiation-driven and that we focus only on the resulting conditions at the sonic point. In comparison with our Ṁ(L) relations in Fig. 8, the Galactic early-type WN model of Lucy & Abbott with [log (L/L_⊙),log (Ṁ/M_⊙ yr^-1)] = [5.45,−4.68] lies 0.2–0.3 dex above our results. The late-type WN and WC models of Vink & de Koter with [5.62,−4.89] and [5.36,−5.34] lie ≈0.1 dex and 0.2–0.3 dex below our relation. However, for these models, much larger radii were adopted and the WN model had a non-zero hydrogen surface mass fraction of X = 0.15.

The WC model by Gräfener & Hamann (2005) was obtained by a simultaneous and locally consistent numerical solution of the wind hydrodynamics and the non-LTE radiative transfer in the co-moving frame. This means that the velocity structure, and in particular ν_∞, were computed within the model, and that the critical-point and wind conditions were fulfilled in the same way as in our present models. Furthermore, the model by Gräfener & Hamann (2005) is equivalent to the models presented here in the sense that it has a smooth/compact wind structure with a sonic point on the hot side of the Fe-opacity peak. With [log (L/L_⊙),log (Ṁ/M_⊙ yr^-1)] = [5.45,−5.14] , the WC model lies near the lower mass limit obtained in the present work, roughly 0.1 dex below the Ṁ(L) relations for WN stars in Fig. 8.

Because of the similarity to our present models, the limitations found in the present work also apply to the models by Gräfener & Hamann (2005). This explains why it was not possible for Gräfener & Hamann to provide models for metallicities below solar, while Vink & de Koter (2005) could provide wind models down to log (Z/Z_⊙) = −4.5 without problems.

5. Conclusions

We examined the physics of optically thick, radiation-driven winds with a smooth wind structure near the sonic point. We focused on H-free WN stars with a compact structure (without envelope inflation) and a high sonic-point temperature. For solar metallicity, our models provide Ṁ(L) relations in good agreement with observations. In particular, their slope agrees well with commonly used mass-loss prescriptions. However, we encountered mass limits that suggest that smooth/compact wind models cannot explain the occurrence of WR-type winds at low masses and/or low metallicities. As suggested by Gräfener & Vink (2013), this could mean that the winds of many WR stars are clumped near the sonic point. Alternatively, the Fe-peak opacities could be enhanced by multi-line effects, or other mechanisms than radiative driving could play a role in WR wind driving.

The presence of physical mass limits for WR-type winds might play a role for the non-detection of low-mass WR stars. These possible SN Ib/c progenitors are believed to be formed through binary interaction, but have not yet been detected spectroscopically. In the absence of strong stellar winds such objects would not display the prototypical WR features, and would become very difficult to detect.

The new sonic-point conditions derived in this work suggest that the ratio P_rad/P_gas at the sonic radius of optically thick winds is of the order of ν_∞/a_s. Because usually ν_∞ ≫ a_s, this means that the formation of optically thick winds is only possible when radiation pressure dominates gas pressure near the sonic radius. This implies that high Eddington factors are a pre-requisite for the formation of optically thick winds.

Otherwise, our models show wind properties that are different from those of OB stars, in particular, they show a qualitatively different dependence on the terminal wind velocity.

The sonic-point conditions derived in the present work can serve as a new means to connect stellar structure models and optically thick winds in a consistent manner at the sonic point. In combination with sonic-point boundary conditions for stellar structure models, as they will be presented in a future work by Grassitelli et al. (2017), such models could help to better understand the interaction between stellar wind and envelope in critical evolutionary phases near the Eddington limit, such as the LBV phase or in phases directly preceding SN explosions.

Acknowledgments

We thank the anonymous referee for the helpful comments and W.-R. Hamann for providing the WRPLOT plotting software. G.G. thanks the Deutsche Forschunsgemeinschaft (DFG) for financial support under grant No. GR 1717/5-1.

References

All Tables

All Figures

Fig. 1 Critical conditions in Prad–Pgas plane. The hydrostatic envelope solution for a 20 M⊙ He star from Gräfener et al. (2012) is indicated by the black/red curve (the axes in this diagram are inverted so that the stellar interior is located in the lower left and the stellar surface in the upper right). On this curve, κ(Prad,Pgas) ≈ κEdd. Areas in the Prad–Pgas plane with higher/lower opacities are indicated in grey/white and by plusandminus signs. The blue line indicates possible locations of the sonic point for an adopted mass-loss rate of 10-4.8M⊙ yr-1 following Eq. (7). The circle indicates the innermost intersection point with an outward-increasing opacity, giving a precise estimate for the location of the sonic point (see text).

In the text

Fig. 2 Mass-loss determination in Prad–Pgas plane. The plot is analogous to Fig. 1, but now the black dashed line indicates the sonic-point conditions of an optically thick wind with an adopted terminal wind velocity of ν∞ = 1800 km s-1 following Eq. (27). Again, the circle indicates the innermost intersection point of Eq. (25) with the static envelope solution. Given the values of Prad and Pgas at the intersection point, Ṁ can be determined from Eq. (7).

In the text

	Fig. 3 Mass-loss determination in Prad–Pgas plane. The plot is analogous to Fig. 2, but now the sonic-point conditions following from our numerical models (indicated by the solid blue curve) are used to determine Ṁ.
In the text

Fig. 4 Sonic-point ratio of Prad/Pgas. Different contributions to the ratio of Prad/Pgas at the sonic point, as imposed by the outer wind, are indicated on a logarithmic scale. The left panel shows the same sequence of wind models as in Fig. 3, i.e., models with a fixed terminal wind velocity of ν∞ = 1800 km s-1 as a function of the adopted mass-loss rate Ṁ. The right panel shows the same model sequence as in Fig. 5, i.e., models for which Ṁ is derived consistently with the critical condition (Eq. (3)) as a function of ν∞. We plot the ratio of Prad/Pgas that follows from our numerical models (solid blue line), from the approximate relation Eq. (25) (black dashed line), and Eq. (31) (red circles). The grey dotted and dash-dotted lines indicate the (weighted) sonic-point optical depths τs and respectively. The red solid line indicates the ratio τs/η, and the red dashed line the factor . The black solid line indicates the offset factor in Eq. (31).

In the text

	Fig. 5 Mass-loss predictions. The derived mass-loss rates for a 20 M⊙ He star are indicated as a function of the adopted terminal wind velocity ν∞. The solid blue curve indicates numerical wind models (cf. Sect. 2.2.3) and the dashed curve our approximate estimates (Sect. 2.2.1).
In the text

	Fig. 6 Terminal wind velocities ν∞ from Hamann et al. (2006), Hainich et al. (2014). Empirical values for single, H-free WN stars in the Galaxy (blue) and LMC (grey) are indicated as a function of the escape velocity νesc. Dashed lines indicate relations with ν∞/ νesc = 1.3 and ν∞ = 1800 km s-1.
In the text

	Fig. 7 Terminal wind velocities ν∞ from Niedzielski & Skorzynski (2002), Niedzielski et al. (2004). Empirical values for a subset of the stars in Fig. 6 are indicated as a function of the escape velocity νesc. Dashed lines indicate relations with ν∞/ νesc = 1.3 and ν∞ = 1800 km s-1.
In the text

Fig. 8 Mass-loss relations for Galactic metallicity. Mass-loss rates for H-free WN stars in the Galaxy as a function of luminosity. Black/grey stars indicate empirical mass-loss rates for Galactic WN stars with/without known distance from cluster membership (Hamann et al. 2006). The empirical mass-loss rates have been scaled down to those expected for a wind clumping factor of D = 10. The red curves indicate mass-loss rates obtained from our models with Z = 0.02 for ν∞/ νesc = 1.3 (top) and ν∞ = 1800 km s-1 (bottom), the labels indicate stellar masses in M⊙. The dotted/dashed lines indicate the empirical mass-loss prescriptions for H-free WN stars from Nugis & Lamers (2000) and Yoon (2017).

In the text

Fig. 9 Mass-loss relations for LMC metallicity. Mass-loss rates for H-free WN stars in the LMC as a function of luminosity. Black stars indicate empirical mass-loss rates from (Hainich et al. 2014). The blue curves indicate mass-loss rates obtained from our models with Z = 0.008 for ν∞/ νesc = 1.3 (top) and ν∞ = 1600 km s-1 (bottom), the labels indicate stellar masses in M⊙. The dotted/dashed lines indicate the empirical mass-loss prescriptions for H-free WN stars from Nugis & Lamers (2000) scaled with the metallicity dependence from Vink & de Koter (2005), and from Yoon (2017), respectively. Notably, our models do not predict solutions below ~25M⊙, which contradicts the observed range of luminosities for H-free WN stars in the LMC.

In the text

	Fig. 10 Limiting masses. Wind and envelope solutions for a 14 M⊙ He star analogous to Fig. 3. For Z = 0.02 the envelope solution (red/black curve) just matches the wind solution with ν∞ = 1800 km s-1 (blue curve). For Z = 0.008 (black dashed curve), no match can be obtained.
In the text