© ESO 2021

1 Introduction

The evolution of stars with initial masses higher than eight times that of the Sun plays an essential role in the chemistry and dynamics of galaxies similar to our Milky Way (Crowther 2007; Doran et al. 2013; Ramachandran et al. 2018; Prantzos et al. 2018). These massive stars are a vital source of heavy elements and UV radiation, enriching their surroundings through strong radiation-driven winds (Lucy & Solomon 1970; Castor et al. 1975; Puls et al. 2008; Vink et al. 2001; Björklund et al. 2021). In the final evolutionary stages, some of those massive stars become hydrogen-depleted, usually core He-burning progenitors of neutron stars or black holes (Yoon et al. 2012; Groh et al. 2013). Among these are so called classical Wolf-Rayet (WR) stars (Wolf & Rayet 1867), which are characterised by strong spectral emission lines and a high luminosity to mass ratio, L∕M_⋆ ~ 10⁴L_⊙∕M_⊙, (Crowther 2007). Classical WR stars are distinct from very massive, main-sequence WR stars, which are core H-burning (de Koter et al. 1997; Crowther et al. 2010). Classical WR stars are also different in their composition and evolutionary state from hydrogen-deficient WR central stars of planetary nebula (e.g. Todt et al. 2010).

As first suggested by Beals (1929), the prominent emission lines visible in WR spectra indicate strong stellar winds with high terminal speeds (v_∞~ 2000−3000 km s⁻¹) and mass-loss rates (Ṁ ~ 10⁻⁵−10⁻³  M_⊙ yr⁻¹) (Hamann et al. 2019; Sander et al. 2019). However, while the overall wind properties of massive main-sequence OB-stars (see Puls et al. 2008 for a review) are rather well reproduced by the line-driven wind theory formulated first by Castor et al. (1975; CAK), this standard theory typically fails to explain the order of magnitude higher mass-loss rates of classical WR stars (Cassinelli 1991; Lamers & Leitherer 1993).

Nonetheless, WR winds are still thought to be radiation driven, as their high L∕M_⋆ ratio brings them close to the limit at which the acceleration due to radiation g_r balances that of gravity g (i.e. close to Γ = 1, for Eddington factor Γ = g_r∕g). In fact, consulting the opacities used in stellar structure calculations, which are generally obtained from tabulations computed assuming a static medium (e.g. Iglesias & Rogers 1996), the Eddington limit for classical WR stars is already reached in sub-surface layers with temperatures of about 150−250 kK. Since convection is highly inefficient in these layers (Gräfener et al. 2012), static models then typically display density inversions and highly inflated stellar envelopes (Ishii et al. 1999; Petrovic et al. 2006; Gräfener et al. 2012; Sanyal et al. 2015). For corresponding stellar evolution calculations, such envelope inflation is computationally difficult to treat and various numerical tricks are thus typically required to make computations tractable (Paxton et al. 2013; Ekström et al. 2012). This then leads to very hot (T_eff ≳100 kK) and compact (R ~ R_⊙) WR stellar surfaces. On the other hand, spectroscopic studies aiming to constrain the WR surface from observationally inferred T_eff typically deduce hydrostatic radii that are factors of ~2−3 higher than predicted by such evolution models (e.g. Crowther 2007; Sander et al. 2019; Hamann et al. 2019); this mismatch is sometimes referred to as the core-radius problem of classical WR stars.

Rather than retaining a static envelope, breaching the Eddington limit in sub-photospheric layers can initiate an optically thick supersonic outflow. Nugis & Lamers (2002) and Grassitelli et al. (2018) suggested that WR mass-loss rates can be computed simply by considering only the conditions at the point where the sonic speed is reached (at Γ ≈ 1, since gas pressure terms typically are very small in comparison). Using the OPAL tables of Rosseland mean opacity (Iglesias & Rogers 1996), Ro & Matzner (2016) showed that supersonic velocities are indeed found in deep sub-photospheric layers, but that a successful wind solution that could bring the initiated mass flux to infinity could not be found.

However, these models neglect the strong enhancement of the line opacity expected in a supersonic outflow. Previous modelling attempts including the Doppler effect in the line opacity calculations have either relied on a pre-assumedfixed velocity field to solve for the mass loss (e.g. see Lucy & Abbott 1993; Springmann 1994; Springmann & Puls 1998; de Koter et al. 1997), or attempted an iterative (assuming time independence) solution towards a self-consistent velocity field and mass loss (Gräfener & Hamann 2005; Sander et al. 2020; Sander & Vink 2020).The latter models have used comoving-frame (CMF) radiative transfer for the calculation of g_r. Such CMF transfer is a computationally intensive numerical technique that requires the velocity field to remain smooth and monotonic (see, e.g. discussion in Sander et al. 2020). Moreover, all these previous studies of classical WR outflows have been performed in the steady-state limit.

In this paper we present a first attempt to build a model that addresses both (dynamic) envelope inflation and line driving, using a hybrid opacity approach based on combining the opacities used for static stellar structure calculations with a simple variant of the standard parameterisation for line opacities in supersonic flows. This formalism then allows for computation of both steady-state and time-dependent WR wind structures. The organisation of the paper is as follows: in Sect. 2, we describe our basic physical set-up. We present dynamical models in the steady-state limit and compare these to corresponding static calculations in Sect. 3. Then in Sect. 4 we compare these steady-state models to full time-dependent radiation-hydrodynamical simulations of dynamically inflated WR outflows. In Sect. 5 we discuss our results and some open questions regarding WR stellar outflows, and finally we summarise our results and provide an outlook for future work in Sect. 6.

2 Physical model

We describe the WR wind outflow by the appropriate hydrodynamical equations of mass and momentum conservation assuming spherical symmetry: (1) (2)

Here ρ, v and P_g = k_BTρ∕(μm_H) are mass density, velocity and gas pressure, where T is gas temperature, k_B is the Boltzmann constant, μ is the mean molecular weight, and m_H is the mass of the hydrogen atom, g = GM_⋆∕r² is the gravitation acceleration for a constant stellar mass M_⋆ and (3)

is the acceleration due to stellar radiation, for radiation flux F and flux-weighted opacity (mass-absorption coefficient) κ_F in cm² g⁻¹. Using the Eddington factor we can also write the equation of motion (e.o.m) as: (4)

The Eddington ratio can be expressed as: (5)

for a stellar luminosity L = 4πr²F. In this paper, we assume that the radiative luminosity remains constant throughout the outflow. This means that we neglect a term corresponding to the work of the radiation field against gravity. We can estimate the corresponding expected luminosity variation by computing the photon-tiring parameter m = Ṁ∕Ṁ_max (Owocki et al. 2017), which is the ratio of the stellar mass loss Ṁ to the maximum amount of mass loss that the stellar luminosity can drive, Ṁ_max = LR_c∕(M_⋆G), where R_c is the core radius defined at the fixed lower boundary. Taking the stellar parameters used throughout this paper, M_⋆ = 10 M_⊙, R_c = 1R_⊙, log₁₀(L∕L_⊙) = 5.416, and a typical order of magnitude mass loss, Ṁ ~ 5.0 × 10⁻⁵ M_⊙ yr⁻¹, one gets (6)

which demonstrates that the luminosity variation is only a marginal effect for the cases considered in this paper.

In practice, especially for time-dependent dynamical computations with possibly non-monotonic velocity fields, it is not computationally feasible to derive the temperature structure, for example, from radiative equilibrium by means of full solutions to the frequency-dependent radiative transfer equations. To simplify,we therefore follow the common approach of replacing the full energy equation by the Lucy (1971) analytic radiative equilibrium model for a grey, spherically symmetric, diluted atmosphere, however, replacing the grey opacity by the actual flux weighted κ_F (e.g. see Lucy & Abbott 1993). This allows us to write the temperature structure as: (7)

where is a core effective temperature, defined here from the lower boundary set at R_c, σ_SB is the Stefan-Boltzmann constant, is the dilution factor, and τ_sp is the spherically modified optical depth (8)

Here again the radius of our fixed lower boundary R_c is used as a scaling radius for the spherically modified optical depth scale. This is similar to the WR models by Sander et al. (2020) who also define their stellar radius from the hydrostatic lower boundary. It is different, however, than, for example, Sundqvist et al. (2019) who set the reference radius for the temperature structure in their O-star models at a photospheric radius R_ph. An argument against using R_ph as a scaling radius for optically thick WR winds is, however, that R_ph ≫ R_c, such that the geometric dilution of the radiation field would be neglected for large parts of the outflow (see also discussion in Nugis & Lamers 2002). In principle, the explicit choice of the reference radius will affect the temperature structure throughout our models. However, test calculations have shown that this is a quite marginal effect for the resulting mass loss and velocities of the simulations considered in this paper. In any case, to facilitate comparison with other models we also introduce a stellar photospheric effective temperature T_{ph, eff} at this photospheric radius R_ph. Since the definition of a stellar photosphere in a spherically diluted stellar envelope is non-trivial¹ and depends on the specific computation of the temperature structure, we here simply approximate (9)

such that the photospheric and core effective temperatures and radii are related through (10)

For given stellar parameters L, M_⋆, and R_c, and a known variation of κ_F, Eqs. (1)–(8) form the basic system of coupled differential equations under investigation in this paper. In general, κ_F is a complicated function that depends on chemical composition, density, temperature, velocity and redial position of the absorbing and emitting gas. In this study, we focus on two cases. First, we consider static opacities such as those used for stellar structure and evolution computations (see Ro & Matzner 2016). Due to the high continuum optical depth of WR outflows, conditions in the wind launching regions resemble those of radiation diffusion, allowing us to replace the flux-weighted opacity with a Rosseland mean, and so to utilise standard OPAL tables for the opacity mapping as a function of T and ρ (Iglesias & Rogers 1996). However, since these tabulations ignore the influence of Doppler shifts on the line opacity (which are critical for WR wind driving, e.g. Sander et al. 2020), we also consider a hybrid model accounting at least approximately for this effect. More specifically, we use a CAK-like parametrisation to add the cumulative force from an ensemble of spectral lines to the OPAL tables (see, Sect. 3.3).

3 Steady-state approximation

To build physical insight we first examine a simplified steady-state case. Equations (1) and (4) are then written as: (11) (12)

for an isothermal sound speed a² = k_BT∕(μm_H). Equation (11) defines the spherically symmetric mass-loss rate Ṁ. Here we assume μ = 4∕3 corresponding to a fully ionised helium plasma.

It is commonly believed that the winds of WR stars are initiated in the deep sub-photospheric layers where the gas temperature is T ~ 150−200 kK. In this temperature region, due to iron recombination, a significant increase of the Rosseland mean opacity is observed, sometimes called the iron opacity bump. A simple comparison of scales reveals that throughout the supersonic wind outflow the first two terms of ther.h.s. in Eq. (12) are much smaller than the radiation and gravitation terms, so that we can approximate the e.o.m with (13)

Here, the sound-speed term on the l.h.s. is retained in order to enable a mapping of the supersonic solution onto the subsonic part.

3.1 Boundary conditions

Equations (7)–(13) impose a two-point boundary value problem, with boundaries at R_c and at² r → ∞. The inner boundary at r = R_c is defined to be a sonic point, with , where a_s is defined to be the sound speed at this sonic point. As can be seen from Eq. (13), this implies Γ(R_c) = 1. Below the sonic point as the velocity rapidly diminishes from its sonic point value a_s inwards, the density structure is well approximated by the hydrostatic solution (e.g. see Ro & Matzner 2016; Grassitelli et al. 2018). In contrast, in the supersonic part, the density structure diverges from the structure given by the hydrostatic solution. As such, we initiate the calculation at the sonic point and only consider the supersonic part of the outflow, assuming that the sub-sonic region still can be described by hydrostatic equilibrium. The tabulated OPAL opacities can then be used to locate points for which Γ = 1, constraining the lower boundary density as a function of temperature. In this way, the range of possible mass-loss rates is set as a function of the lower boundary temperature.

The outer boundary is set by requiring that the wind temperature decreases sufficiently at large radii. Replacing density with the mass-loss rate using Eq. (11), the spherically modified optical depth (Eq. (8)) is used to compute the temperature structure according to Eq. (7). We compute the optical depth at the outer boundary τ_out by radially integrating Eq. (8) in r ∈ [R_max,  ∞), assuming a terminal velocity v(r ≥ R_max) = v(R_max) = v_∞ and κ_F(r ≥ R_max) = κ_e, where κ_e is the Thompson scattering mass absorption coefficient. The outer boundary temperature is then set by the maximum between applying τ_out in Eq. (8) and a floor value, which is typically set to a few tenths of the stellar effective temperature (see Puls et al. 2005; Sundqvist et al. 2019). By matching the inner and outer boundary constraints, we uniquely constrain the wind structure. In turn this then sets the unique mass-loss rate of the model from the range of possible solutions at the lower boundary.

3.2 Energy requirement

Additionally, for wind material to escape the stellar potential, the global energy requirement must be satisfied, which means that the integrated mechanical energy of the wind has to be positive. Taking the supersonic limit a² ∕v² ≪ 1 in Eq. (13) yields: (14)

Here v₀ ~ a_s is the initial velocity and is the escape velocity from R_c. The first term on the r.h.s. corresponds to the work done by the stellar radiation and the second term is the work done by gravity. We can then introduce as the net energy change in units of : (15)

For the wind to escape to r →∞, this net energy change has to be positive, that is . Alternatively, we can also introduce the ratio of kinetic to core potential energy as a square of the velocity to core escape velocity ratio: (16)

which goes to as r→∞ so that by construction we require w(r →∞) > 0 in order to escape.

3.3 Failed Winds in the static opacity limit

We first investigate the possibility of driving the stellar outflow in the static opacity limit, i.e. the limit in which we only consider OPAL opacities, thus neglecting the Doppler effect; κ_F ≈ κ_OPAL (see also Ro & Matzner 2016). To set up the problem, we choose the stellar mass and the core radius (as given in Table 1) such that they roughly represent the average value of the mass-to-radius ratio of the observationally inferred WR stellar population from Hamann et al. (2019) and Sander et al. (2019). The luminosity is then set such that the resulting Thompson scattering Eddington ratio Γ_e = 0.40, which also corresponds approximately to the observationally inferred average (Hamann et al. 2019; Sander et al. 2019). All models presented inSect. 3 and 4 of this paper use these stellar parameters, and also the OPAL tables for solar metallicity Z = 0.02 and composition by Grevesse & Noels (1993). This set-up then gives a range of possible sonic-point values of such that the condition is fulfilled (see Fig. 1a). These points are used as a single point boundary condition, after which forward integration from the sonic point is carried out. The results of these integrations for various fixed mass-loss rates Ṁ = 4πR_cρ_sa_s are shown in Fig. 1a. The figure displays a colour map of Γ with the temperature on the abscissa and the ratio between the radiation pressure (17)

and gas pressure P_g on the ordinate (the ratio P_r∕P_g is a proxy for the density). The integration curves show the density and temperature structure for the different cases, starting from various positions on the Γ = 1 curve at the inner boundary R_c and extending outwards. Figure 1b then shows the corresponding net gain of kinetic energy (Eq. (16)), demonstrating that after an initial acceleration and increasing velocities, all of the curves start to decelerate until eventually reaching zero velocity and so terminating the outward integration (some solutions had to be formally terminated before actually reaching zero velocity, due to numerical difficulties in the steep deceleration region). This indicates that none of the potential solutions starting from the Γ = 1 curve are able to escape the stellar gravitation potential. These results are consistent with Ro & Matzner (2016), who also found that OPAL opacities were not able to sustain a radiation-driven mass loss initiated in the deep layers around the iron bump, and indicates that an additional source of opacity is required, raising the question of a missing force.

Fig. 1 (a) Colour map of ΓOPAL for the stellar parameters from Table 1. Temperature is given on the abscissa and radiation-to-gas pressure ratio on the ordinate. is identified with a white contour. The colour bar on the right corresponds to the colour coding of the plot and gives the numerical value of ΓOPAL for a given temperature and radiation-to-gas pressure ratio. Different line styles show the radiation-to-gas pressure ratio and temperature structure for the different fixed mass-loss rates. (b) Ratio of kinetic to core potential energy of the corresponding solutions from panel a. The abscissa here shows the radius coordinate x.

3.4 Missing force

The natural candidate for the missing opacity source discussed above comes from the cumulative effect of spectral lines due to the Doppler effect. To estimate this we use the CAK parametrisation for a distribution of spectral lines, using here a simple approximation of taking the absolute value of the velocity gradient: (18)

where is an effective strength of the line ensemble (Gayley 1995), and α sets the CAK line distribution power index. In reality, once regions where d v ∕dr < 0 arise, it is possible for a photon proceeding from the stellar core to be Doppler shifted into resonance with the same spectral line at multiple locations in the outflow, thereby requiring a fully non-local computation of the line acceleration. As discussed further below, this could (in principle) be handled in a more self-consistent manner by, for example, Monte-Carlo (MC) line-force calculations.In this initial study, however, we opt for a simpler approach, in order to make time-dependent simulations feasible as well. Two possible approximations that can be taken then are to either take the absolute value of the velocity gradient or use max(0, dv∕dr). The first approach estimates (approximately) an upper limit to the CAK-like force in regions with negative velocity gradients, whereas the second approach would correspond to a lower limit. As most previous studies of line-driven dynamical flows have used (for example, in investigations of line-driven winds from rapidly rotating stars, Petrenz & Puls 2000, and around discs, Proga et al. 1998; Kee et al. 2016), we follow this standard approach in this paper. However, test computations using instead the lower limit max (0, dv∕dr) show that although some details of the radiation force profile then are changed, none of our conclusions are affected. With this approximation, using Eq. (18) for Γ_CAK then accounts for the force from spectral lines in supersonic regions, while Γ_OPAL handles the force from the continuum and spectral lines in the static limit. As such, the total radiative acceleration can be simply estimated by the sum of the CAK-like contribution and the static OPAL opacities, Γ_tot = Γ_OPAL + Γ_CAK.

Let us note here directly that the CAK line force is originally derived for application in winds that are optically thin for continuum radiation. In the lower layers of the dense WR outflows considered here, this will typically not be the case. A more rigorous treatment would then also need to take into account the change in continuum intensity by solving the transfer equation for line+continuum opacities. While approximations for the diffusion limit exist (Gayley et al. 1995) no general formalism has been developed. As such, in this first study of WR wind dynamics we opt for the simple radial-streaming CAK form above.

Furthermore, due to the explicit dependency of the CAK line force on d v∕dr the sonic point is now formally not a critical point (see Castor et al. 1975). However, our iteration scheme circumvents this issue by always applying the velocity gradient from the previous iteration, meaning that Γ_tot(r) technically is no longer an explicit function of dv∕dr. Nevertheless, since we require dv∕dr to be converged between iterations the corresponding feedback upon Γ_tot(r) is still (implicitly) accounted for. The agreement between the converged structures of these steady-state models and the time-dependent simulations presented in the next section brings further support to this method.

A key question then becomes whether Γ_CAK also has a significant impact on the conditions at the optically thick lower boundary at the sonic point. Assuming typical O-star values α =0.66, (e.g. Puls et al. 2000), approximating the velocity gradient³ as d v∕dr ~ Δv_max∕ΔR(v = v_max), and using the mass density and velocity computed in the static opacity limit for Ṁ = 5.0 × 10⁻⁵ M_⊙ yr⁻¹, we find the Eddington ratio for the CAK force at the lower boundary to be Γ_CAK = g_CAK∕g ~ 0.1. Such a small contribution from the CAK force is not surprising due to the high density at the lower boundary⁴ (Γ_CAK ~ 1∕ρ^α, with α >0). However, from test calculations it turns out that although the CAK contribution indeed increases as we move outwards in the wind, a simple model with constant α = 0.66, is still not able to provide the force necessary to achieve a positive total wind energy and so drive the wind to infinity (see previous section).

There are, however, several indications that the line force in the outer winds of classical WR stars indeed might be significantly enhanced as compared to O-type stars. For example, the independent MC line-force calculations by Lucy & Abbott (1993) and Springmann & Puls (1998) both find global momentum rates for WR stars that exceed those of O-stars by an order of magnitude or more. In these MC computations, the outer-wind enhancement of the line force stems from the decreasing level of ionisation with increasing distance from the star in the optically thick WR outflows. Due to this shift from high-to-low ionisation stages as photons move outwards in the wind, new spectral lines become available for photons to interact with. As these new lines have different interaction frequencies than those previously available in the inner wind parts, this means that the line-frequency gaps typically observed for the nearly constant (frozen-in) ionisation conditions of optically thin O-star winds are effectively closed. This then leads to photon trapping (Lucy & Abbott 1993) and efficient multi-line scattering, which significantly increases the line force.

To mimic the effect of an outer wind line-force enhancement, we introduce here a simple spatial variation α =α(r), which increases the CAK force away from the star by a modest decrease of α. To this end, we assume α(r) as a piecewise linear function with a constant α_max from the lower boundary to r₁ and a constant α_min from r₂ to the outer boundary, interpolating between the two values in the range [r₁, r₂]. This then allows the wind outflow initiated at the lower boundary to achieve positive total energy, and so be sustained until the upper boundary. While the specific treatment of α is done here in an ad hoc manner, the overall assumed form is quite consistent with the multi-scattering line-force enhancement found in the MC line-force computations for a CAK-like spectral line distribution by Springmann (1994). The approximate values of α(x) that correspond to the line force computed by Springmann (1994) are given in Fig. 2 as black asterisks, with the α(x) as adapted in the standard model of this paper over-plotted as a solid line.

The figure shows that the order of magnitude outer-wind line-force enhancement applied here indeed is on the same order as that found in these MC calculations.

We note, however, that the Springmann (1994) calculations were performed only for parameterised (monotonic) β-type velocity laws. As such, there is not a complete one-to-one correspondence between the spatial and velocity scales in these MC calculations and those of the present paper. Nonetheless, the overall similarity with the very simple radial dependence of the α parameter adopted here is encouraging. In a follow-up work, we plan to couple the hydrodynamic simulations presented here directly to such MC computations, building on the method developed by Puls et al. (2000) for proper estimation of occupation numbers, and accounting fully for multi-scattering effects as well as non-monotonic velocity fields.

As outlined in detail in the following sections (for our standard values listed in Table 2) the formalism outlined above gives rise to dynamically inflated wind solutions, where the structure in the optically thick deep layers near the iron bump is primarily governed by Γ_OPAL but where Γ_CAK takes over the driving in the diluted outer regions. This then causes a typical wind structure that is characterised by a high mass-loss rate ignited from the lower boundary, a slow (end even partially negative) acceleration of the deeper layers, and a rapid re-acceleration of the outer wind. As such, the typical velocity laws we find not only drastically deviate from the standard β-type laws assumed by most atmospheric models, but they also contain a significant region of negative acceleration (see also Fig. 14 in Pauldrach et al. 1993).

Fig. 2 Plot that demonstrates values of α (see text) required to reproduce the line force computed from MC radiative transfer calculations by Springmann (1994; asterisks). The solid line then over-plots the piecewise linear function adopted in this paper to describe the corresponding variation of α. The vertical dashed lines correspond to r1 (x1) and r2 (x2).

3.5 Dynamically inflated steady-state model

Based on the hybrid formalism introduced above we can now numerically solve Eqs. (7)–(13), iterating towards convergence using a simple scheme based on Runge-Kutta integrations. As discussed in Sect. 3.1, to constrain the mass-loss rate it is critical to here take into account the two-point boundary value nature of the problem. At a basic level, every iteration is performed using two steps: (i) an inside-out integration R_c → R_max, where the velocity, density, and temperature structure are computed, followed by (ii) an outside-in integration R_max → R_c updating the temperature structure using the correct boundary condition from the optical depth found in the first step. The mass-loss rate is then updated for the next iteration by finding the correct location on the Γ =1 curve at thelower boundary, using the updated temperature structure from the outside-in integration.

Using this scheme, Fig. 3a shows the final, converged temperature profile for our standard stellar parameters (Table 1) and choice of line-force parameters (Table 2). The final (also converged) mass-loss rate for the model is Ṁ = 1.47 × 10⁻⁵ M_⊙ yr⁻¹.

In addition to a high predicted mass-loss rate, the solution displays a non-monotonic velocity field. The resulting velocity profile shown in Fig. 3b can be split into three basic regions: (1) an initial acceleration from the core followed by (2) stagnation and even a deceleration part (and so a negative velocity gradient, marked with crosses in the figure), and (3) an outer region of fast re-acceleration to outflow velocities > 1000 km s⁻¹. The distinct behaviour of these three regions can be understood via the different contributing parts of the total radiation force.

Figure 3c illustrates this, plotting the individual contributions of OPAL and CAK-like opacities to the total radiative acceleration. From the figure, on one hand, it is clear that in the inner regions the radiation force is dominated by OPAL opacities and the CAK-like contribution is small; indeed, the sonic point conditions at the lower boundary can here be quite well approximated using only OPAL opacities, i.e. neglecting any influence from the Doppler-shift enhancement of line opacities. On the other hand, in the outer regions Γ_CAK becomes thedominant part, causing a re-acceleration of the flow and setting a high terminal wind velocity. This re-acceleration ensures thatthe energy requirement of Eq. (15) is satisfied. As such, it is this part that ultimately allows the wind to be sustained all the way to the outer boundary, in contrast to the failed outflow and fallback seen in the previous section.

To explore the influence of different opacity contributions on the mass loss of the model, we explore the logarithm of the spherically modified optical depth (Eq. (8)) as a function of radius in Fig. 3d. The dashed vertical line marks the location of the stellar photosphere defined by Eq. (9). We note, however, that at the photosphere τ_sp (R_ph)≠2∕3 as R_c is used as ascaling radius in Eq. (8). For most stars, this photospheric radius R_ph serves as the hydrostatic stellar radius; by contrast, for the WR stars examined in this paper, the photosphere is located well in the outflowing regions. Inspection of Fig. 3d shows that the majority of τ_sp is accumulated within few stellar radii in the regions dominated by OPAL opacities, thus the converged lower boundary temperature, and thus the mass-loss rate, is primarily controlled by Γ_OPAL.

Inspection of the density distribution in Fig. 3e further clarifies the terminology dynamic inflation. In this figure, we have again identified the location R_ph with a vertical dashed line. The figure shows a relatively slow decline of the density in the sub-photospheric parts, which is somewhat reminiscent of the profiles seen in hydrostatic stellar models near the Eddington limit that undergo envelope inflation (see further below). But unlike such static inflation mod0els, the simulations here allow for a velocity field to develop; the low-density envelope now represents a dynamically inflated star.

3.6 Comparison to static inflation

For the same lower boundary conditions as adopted above, we can also compute static models. This can then serve as a direct comparison between the dynamic inflation found in the models above, and the static inflation that often occurs in stellar structure models neglecting the v dv∕dr term in the e.o.m. (and so assuming a hydrostatic stellar envelope throughout, e.g. Petrovic et al. 2006; Gräfener et al. 2012).

To construct such static models, we simply set v ~dv∕dr ~ 0 in our steady-state e.o.m. above. Introducing a column mass (19)

we write the resulting equation of hydrostatic equilibrium as: (20)

We further simplify the temperature calculation by assuming radiative diffusion in these optically thick hydrostatic layers: (21)

for a radiation pressure according to Eq. (17). Combining the equations for gas and radiation pressure, one obtains (22)

where now Γ = Γ_OPAL throughout the complete model.

This set-up imposes a similar two-point boundary value problem as previously analysed. To provide a fair comparison with those dynamical simulations, we assume the same lower boundary radius r = R_⊙ as before and re-use the lower boundary temperature T ≈ 230 kK found from the converged hydrodynamic steady-state model presented above. However, at this inner boundary we now impose Γ < 1 as Γ ≥ 1 would lead to non-zero gas pressure at r →∞ (see also Sect. 2.2 in Gräfener et al. 2012). The outer boundary conditions for gas and radiation pressure (i.e., density and temperature) are also set up in the same way as before. Assuming Γ = const = Γ_e for the region r ∈ [R_max,  ∞), we may integrate the hydrostatic equations to estimate the temperature and the density at the outer boundary r = R_max. This fixes the effective temperature and the spherically modified optical depth scale of these hydrostatic (HS) models, requiring that the outer boundary always satisfies . Thus in this model, the stellar photospheric radius is defined as (23)

for the corresponding hydrostatic effective temperature (24)

quite analogous to how these quantities typically are defined in static stellar structure and evolution models.

These boundary conditions are used to numerically solve the new set of equations from the outer boundary inwards. A consistent solution is obtained by shooting for the fixed inner boundary temperature and radius, varying the stellar photospheric radius and assuming inner-boundary values as discussed above. In this way we obtain static envelopes corresponding directly to our previous dynamic models. The right panel of Fig. 4 shows the predicted densities for static models with different assumed⁵ Γ_e ∈ [0.37, 0.40, 0.41] (for comparison see Fig. 1 in Gräfener et al. 2012). These densities are then compared to those of the steady-state computation in Sect. 3, for which Γ_e = 0.40.

We note directly that in the static case an extended region of low density is formed. The photospheric radius of this inflated envelope is extremely sensitive to Γ_e, with growing notably as Γ_e is increased. For the static model with highest assumed Γ_e = 0.41, we already find a large accompanied by a relatively cool effective temperature T_{HS, eff} = 21 kK; further increasing Γ_e causes the inferred photospheric radius to inflate to increasingly unrealistic, extremely large values. For the Γ_e = 0.40 that we use in our dynamic models we find and T_{HS, eff} = 56 kK. The reason for the inflation here is the same as in the models by Petrovic et al. (2006) and Gräfener et al. (2012). Basically, it is related to the fact that for high values of Γ_e the integration moves along a Γ ≈ 1 curve, forcing the density to remain low and almost constant over a large radial extent. Since enough column mass (or equivalently optical depth) still has to be accumulated in order to meet the lower boundary condition, the only way the star can react is by expanding in radius (see also discussion in Owocki 2014). Comparing this now to our dynamic solutions shows that, for the same Γ_e, inclusion of the dynamical terms reduces the inflation somewhat. Nonetheless, R_ph still lies very far away from the core in such outflow models of dynamically inflated envelopes; indeed, for our standard Γ_e = 0.40 we obtain R_ph = 4.48 R_⊙ and T_{ph, eff} = 62 kK. Such a close agreement of photospheric radii in dynamic and static models is however not generally found. Namely, a very small increase to Γ_e = 0.41 in the static model leads to an order of magnitude increase of the star’s photospheric radius (see Fig. 4), in contrast to the dynamic model where this barely affects the location of R_ph.

The analysis above assumes that all luminosity is carried by radiation. In stellar models, the strong increase in opacity around the iron bump will typically lead to an onset of convection, such that L = L_rad + L_conv. In principle, this may alleviate the inflation by allowing a portion of the luminosity to be carried by convection, lowering the effective Γ_e ~ L_rad. However, for the conditions prevailing in WR stars convection in these near-surface layers should be very inefficient. Following Gräfener et al. (2012), an upper limit to the fraction of the total stellar flux F that can be carried by convection is: (25)

This means that for the T_{c, eff} > 10⁵ kK considered in this paper, less than a percent of the total luminosity will be carried by convection.

This convective inefficiency is confirmed by the computation of WR stars using the stellar structure and evolution code package MESA (Paxton et al. 2019). Namely, using the standard mixing-length theory (MLT) prescription for convection shows that such stellar structure and evolution computations⁶ indeed undergo similar inflation in the WR stage as seen here in Fig. 4. However, in these stellar evolution simulations envelope inflation and the associated density inversions often imply prohibitively short time steps within the model calculations. To prevent this from happening, such stellar models often invoke an additional energy transport that is far greater than that implied by the upper-limit estimate of L_conv above (see, for example, the so-called MLT++ prescription in MESA, Sect. 7.2 in Paxton et al. 2013). This then forces L_rad to always be low enough such that the calculation can proceed.

To mimic these tricks, in the next step we artificially reduce Γ_e ~ L_rad by introducing L = L_conv + L_rad. Then using L_rad in Eq. (20), while computing the radiation pressure from the total luminosity in Eq. (21), we solve the differential equations exactly as above. By then assuming different ratios L∕L_rad for the same fixed total luminosity, we can directly compare dynamic and static inflation models with and without such a reduction in radiative luminosity. Analogous to the left panel, the right panel of Fig. 4 shows the radial density structure of a case with reduced radiative luminosity L∕L_rad = 4 (solid line), comparing this to the original case L∕L_rad = 1 (solid, marked line). The ratio of L∕L_rad = 4 here is chosen such that the stellar envelope no longer experiences a density inversion. The figure demonstrates clearly that the two solutions have vastly different radial scales; while the original model had and T_{HS, eff} = 56 kK, the model with reduced radiative luminosity has a very hot T_{HS, eff} = 121 kK and a photospheric radius that lies less than a percent above the core R_c = R_⊙.

The reason for these drastically different stellar photospheric radii is the same as described earlier; including an additional (artificial) energy transport reduces the effective Γ_e ~ L_rad, so that a dense exponential atmosphere with small scale height is formed. This then avoids inflation and density inversions, and the upshot is a hot and compact WR stellar surface. This is a good example of how envelope inflation and density inversions are often avoided in stellar structure computations by reducing L_rad well beyond the limits implied by standard convective energy transport.

A key result of our analysis here is thus that such tricks of the near-surface regions in WR stars are necessary only because of enforcing a hydrostatic solution. As shown above, when including the dynamical terms the stellar envelope will instead quite naturally develop a dynamically inflated outflow region.

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

Fig. 4 Density profiles of different static and steady-state solutions. On both plots the solid lines show the different cases of hydrostatic solutions and dashed lines show the sub-photospheric part of the steady-state outflow (see Fig. 3e). The different markers then compare hydrostatic solutions for different assumed Thompson scattering Eddington factor Γe (left panel) and for different assumed L∕Lrad ratios for thefixed total luminosity log10(L∕L⊙) = 5.416 (right panel). The numbers at the right end of each density profile correspond to the photospheric effective temperature of the corresponding model in kilo-Kelvins, as defined in the text. To better visualise the innermost parts of the simulations, the right panel displays the logarithm of the scaled radius x = 1 − Rc∕r on the abscissa.

4 Time-dependent numerical hydrodynamic simulations

The previous section presented calculations in the steady-state limit, using various degrees of approximations for solving the dynamical (and static) equations. Building on the same basic formalism, this section now presents full time-dependent numerical hydrodynamics simulations for spherically symmetric WR wind outflows. A key objective is to examine whether such time-dependent simulations relax to similar solutions as the simplified steady models of Sect. 3.

The hydrodynamic simulations are performed using the finite volume code MPI-AMRVAC (Xia et al. 2018) to solve Eqs. (1), (4), and (7). We use the HLL solver (Harten 1983) with a MINMOD flux limiter (van Leer 1979). The radiative acceleration is computed as in Sect. 3.3 by summing up the contributions from OPAL and CAK-like opacities, Γ_tot = Γ_OPAL + Γ_CAK. Here the OPAL contribution is trivially obtained from the local density and temperature, and the CAK contribution from a simple central finite difference scheme applied to Eq. (18). To avoid information propagating over multiple cells in a single integration we further limit the integration time step d t by:

where dt_CFL is the time step set by the standard Courant-Friedrichs-Lewy condition (Courant et al. 1928), and Δr the width of one numerical cell.

As Eq. (7) is not the standard way of computing the energy balance in MPI-AMRVAC we have implemented an additional routine to handle this. In this subroutine the optical depth Eq. (8) is first computed by trapezoid integration starting from the outer boundary i = i_max, setting τ(i_max) = τ_out as before, inwards to i = i_min. The radial temperaturestructure T_i is then updated at each hydrodynamical time step according to Eq. (7). For simplicity, we here define R_c to always be at the lower boundary of the simulation.

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

4.1 Initial and boundary conditions

Analogous to our steady-state simulations, we define the outer boundary to be located at some radius R_max and use the same optical depth and temperature conditions as described in the previous section, Eq. (7) and the pre-specifiedminimum temperature (see above). As in, for instance, Driessen et al. (2019), the numerical outer boundary conditions on the velocity and density are set by simple extrapolation from the outer most points of the simulation grid. To properly resolve the sonic point we set the inner boundary to a lower radius than the expected sonic point R_a, i.e. now R_min ≡ R_c = 0.99 R_⊙ < R_a. At this R_min we then fix the lower boundary density ρ_⋆ and allow the velocity to adjust to the overlying wind conditions (see, e.g. Owocki et al. 1988). ρ_⋆ is initially estimated from our corresponding steady-state models by simply assuming hydrostatic equilibrium below the sonic point. The stellar mass and luminosity are still as in Table 1.

The initial wind condition is an analytic outflow structure, in the form of a simple β-law: (26)

where the terminal velocity , and β = 0.5 are derived from a simple assumed Γ = const > 1. The parameter b sets the minimum initial velocity v_min at r = R_c, here chosen such that v_min ≈ 50 km s⁻¹. The equation of continuity then directly provides the initial density structure. Finally, the initial temperature structure is computed from Eq. (7) with the opacity set according to the assumed input Γ from above. Using these initial and boundary conditions, numerical time integration is carried out until the wind relaxes to a (quasi) steady state. We, however, note that the results of our simulations are not sensitive to these assumed initial conditions.

4.2 Time-dependent dynamic inflation model

Starting from the initially assumed structure the wind outflow gradually relaxes to a hydrodynamically consistent structure. The velocity relaxation is displayed in Fig. 5 as a surface plot, where along one horizontal axis the physical time t in ks of the simulation is shown and along the other the radius in x = 1 − R_c∕r. The vertical axis then measures velocity (in km s⁻¹). At t = 0 the initial velocity structure, v_β, is visible, however as time increases this β-type structure gradually evolves towards the three-phased velocity profile discussed in Sect. 3.5. Around at t ≈ 40 ks the simulation has relaxed to a state very similar to that found in our corresponding steady-state model (Sect. 3.5).

Running the simulation for in total 600 ks, Figs. 6a and b show the Eddington ratio Γ and the velocity profile of the final time step. From Fig. 6a we observe that, just as in the steady-state model, the CAK-like force near the lower boundary on average only has a small contribution to the total force; contrarily, near the outer boundary, itis dominant. The radiation force from the OPAL contribution has the opposite behaviour, dominating at the sonic point and inthe inner regions, but becoming secondary in the outer parts. This then again leads to a velocity profile (solid line in Fig. 6b) with a non-monotonic behaviour, which clearly can not be fit with the simple β-type velocity law assumed in the initial conditions (dotted line).

This velocity profile can be directly compared to that found in our corresponding steady-state model, illustrated by the dashed line in Fig. 6b, which shows close agreement between the two. The mass loss of the time-dependent model is Ṁ = 1.64 × 10⁻⁵ M_⊙ yr⁻¹, which also agrees well with the mass-loss rate of Ṁ = 1.47 × 10⁻⁵ M_⊙ yr⁻¹ from our steady-state model (see Sect. 3). In this way, our time-dependent models demonstrate that the dynamicallyinflated solution discussed in the previous section indeed works as a stable attractor for the stellar wind envelope, confirming our basic conclusions regarding a wind outflow controlled primarily by OPAL opacities in the inner regions, and by the CAK-like force in the outer. The small differences in velocity magnitude and mass loss between steady-state and time-dependent models are mainly related to differences in setting the lower boundary. In steady-state computations we define the sonic point to be exactly at the lower boundary, whereas in time-dependent computations it is allowed to self adjust so that the sub-sonic part is also resolved.

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

4.3 Ṁ dependency on CAK α

A general finding of our baseline simulations in Sects. 3 and 4 is that (as long as the wind can escape), the mass-loss rate is primarily controlled by the conditions in the lower wind; in contrast, the terminal wind speed is controlled mainly by the outer wind conditions. Moreover, while it has turned out that a steady-state model, such as the one presented inSect. 3, is quite difficult to converge, the time-dependent simulations of this section are computationally somewhat easier to handle.

Therefore, we have used our time-dependent set-up to further explore the impact of the CAK-like force parameters on the predictions of Ṁ and v_∞. We investigate the impact of these parameters by varying only α_min and α_max, however we point out that comparable effects could also be achieved by instead introducing corresponding radial variations of the second line-force parameter . To this end, we test a range of values α_max ∈ [0.60, 0.66] and α_min ∈ [0.48, 0.53], corresponding roughly to the scatter of α observed in the MC simulations displayed in Fig. 2. The results of these tests are summarised in Tables 3 and 4, along with the mass-loss rates and terminal velocities found for the different model runs. In addition, we also tested variations in the spatial position of the α_max to α_min transition for the standard CAK force parameters in the range x ∈ [0.3, 0.7]. However, for these tests, no significant influence on the mass-loss rates of the models have been observed. Therefore, the following subsections focus exclusively on varying the values of α_max and α_min.

4.3.1 Dependency on α_min

As Γ_CAK is the prevalent driver in the outer wind, one might initially expect a similar dependence of the mass-loss rate on α_min as given by the standard CAK model. The CAK expression for the mass-loss rate Ṁ_CAK, typically derived from critical-point analyses for Γ_tot = Γ_e + Γ_CAK in the radial streaming limit, is: (27)

From this expression, it directly follows that Ṁ_CAK increases steeply if α is decreased while keeping all other parameters constant. However, it can be seen from Table 3 that in our models Ṁ varies only marginally across the range of considered α_min values. This distinguishes these optically thick WR simulations from standard CAK-type models applied to the optically thin winds of O-type stars. We further confirm the finding from the steady models that if α_min is increased too much, the line force becomes insufficient to sustain the outflow. This happens because the line acceleration in the outer wind must be sufficient to prevent fallback of the mass loss initiated at the Fe-opacity bump. That is, the line acceleration must provide enough additional energy such that the previously introduced condition is fulfilled. As alluded to above, the CAK mass-loss rate formula cannot be directly used to predict the mass loss of this hybrid model. Nevertheless, the limiting value of α_min that can sustain the wind can be estimated by inverting the above equation, solving for the α at which the CAK mass-loss rate equals that of our actual model. For the given model with Ṁ = 1.64 × 10⁻⁵ M_⊙ yr⁻¹, one finds α_min ≈ 0.52. The results of our numerical simulations, summarised in Table 3, indeed suggest that this provides a reasonable estimate for the limiting α_min value.

This is also supported by closer inspection of the simulations with α_min = 0.52, which results in a flattening of the wind velocity profile as . With a further increase to α_min = 0.53, we find thatthe CAK-like force is unable to keep up with the mass flux initiated at the lower boundary. In the 1D simulations presented here, this then leads to a similar failed wind situation as discussed in the previous section; we discuss potential implications of this behaviour for multi-dimensional versions of our simulations at the end of the next section.

Finally, as summarised in Table 3 the terminal velocity indeed is very sensitive to small changes of α_min. Overall, this thus suggests that within this formalism, Ṁ and v_∞ are determined quite independently. As discussed further in Sect. 5.2, such a basic decoupling of Ṁ and v_∞ indeed seems to be supported by observational data of WR winds (e.g. see Hamann et al. 2019; Sander et al. 2019). On the other hand, this result also seems to contradict that of the most recent study from Sander & Vink (2020), who find a correlation between Ṁ and v_∞. We further discuss possible reasons for this in Sect. 5.4.

4.3.2 Dependency on α_max

In contrast to α_min, α_max has little to no effect on determining the terminal velocity of the outflow. Nevertheless, it can affect Ṁ. This is illustrated in Table 4, where a modest increase in Ṁ is found for decreasing values of α_max. We note, however, that again the magnitude of the mass loss variation is small in comparison to that seen in CAK-type wind models of optically thin O-stars.

This rather weak dependence on the CAK force stems from the fact that the dominant contribution to the total radiation force comes from the OPAL opacities in the regions near the stellar core (see also Sect. 3). At the sonic point for α_max = 0.6, we find Γ_CAK ≈ 0.14 and for the α_max = 0.66 we find Γ = 0.02, which is much lower even than our initial order of magnitude estimate Γ ≈ 0.1 in Sect. 3.4. This further supports that the CAK force has a negligible contribution at the sonic point. Moreover, since those inner parts accumulate the majority of the integrated wind optical depth, they also primarily control the lower-boundary temperature to which the overlying wind adjusts, and thereby also the mass-loss rate.

5 Discussion

5.1 Implication for core radius problem

The dynamic models presented in Sects. 3 and 4 provide a natural explanation to the so-called core radius problem of classical WR stars (see Sect. 1). To understand this, we introduce the standard Rosseland mean optical depth scale (28)

where κ_Ros is computed using the Rosseland mean opacities from the OPAL tabulation, such that κ_Ros = κ_OPAL throughout this work.

Then, following the standard spectroscopic model approach of fitting an assumed v_β velocity profile to the entire supersonic region, Fig. 7 directly illustrates how the hydrostatic core (where v_β reaches sub-sonic velocities) in that case would be located at a radius . Indeed, the core radius recovered in this way is a factor ≈2.6 times larger than the true value, demonstrating the basic challenge of determining WR hydrostatic radii from spectroscopic data (e.g. Crowther 2007). The discrepancy essentially stems from the fact that τ_Ros = 2∕3 is reached already at very high radii (see the grey area in Fig. 7), so that only the outermost velocity law is visible to us. Observational estimates of the hydrostatic stellar core radius must thus be obtained by some extrapolation of the velocity law assumed in the spectroscopic model, without any knowledge of the true velocity field in the invisible inner regions.

The dynamic models presented in this paper explain this long-standing discrepancy simply by their extended (dynamically inflated) region of low-velocity outflowing material above the hydrostatic stellar core.

To further illustrate this effect, let us inspect the observational positioning of classical WR stars on the Hertzsprung–Russell diagram (HRD), where the stars are placed according to their observationally inferred luminosities and hydrostatic core effective temperatures. Figure 8 here shows the HRD where the observed galactic sample of WR stars by Hamann et al. (2019) and Sander et al. (2019) is placed. Positions of these stars on the HRD can then be compared to the location of the so-called helium zero age main sequence (He-ZAMS; dashed line in figure), computed here with MESA using the standard mixing-length theory for convection with solar metallicity and relative metal fractions as given by Asplund et al. (2009)⁷. On average, observationally inferred positions of WR stars are then expected to lay to the left of this He-ZAMS line. This is because typical stellar evolution models, where an additional energy transport is typically also invoked in order to prevent the WR cores from inflating (see Sect. 3.5), occupy the region to the hotter side of He-ZAMS on the HRD. Figure 8, however, shows that the average position of these WR stars (indicated by + in the figure) is located on the cooler side of the He-ZAMS instead. To improve the agreement between spectroscopy and evolution predictions a simple correction factor ≈ 3 can be applied to the spectroscopic radii, which then shifts the core effective temperature location (now identified on Fig. 8 by x). And indeed, such a correction factor matches quite well with the correction factor implied by our dynamic inflation models above.

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

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

Fig. 9 Distribution of mass-loss rates and terminal velocities derived from the spectroscopic fitting by Sander et al. (2019) and Hamann et al. (2019) for an observed stellar population of galactic WR stars of nitrogen (WN), carbon (WC), and oxygen sequences (WO). Mass-loss rates and terminal velocities derived in our steady-state (log10 (Ṁ∕ M⊙ yr−1) = −4.83, v∞ = 2200 km s−1) and time-dependent (log10(Ṁ∕ M⊙ yr−1) = −4.78, v∞ = 2063 km s−1) simulations are consistent with the range of mass loss shown on this figure.

5.2 Empirical constraints and model dependencies of mass loss and terminal velocity

The basic model presented in the previous sections shows an overall quite good agreement with empirical constraints on terminal velocities and mass-loss rates of WR stars. Figure 9 plots such empirically inferred values of Ṁ and v_∞, again from Hamann et al. (2019) and Sander et al. (2019). Compared to the mass-loss rates and terminal velocities predicted by our hydrodynamic and steady-state models (summarised in Table 5), the latter is located at the high end of the inferred values. However, since v_∞ is primarily determined by the line force in the outer wind regions, which here is parametrised by the simple radial streaming CAK-approximation outlined in previous sections, a more detailed approach will be necessary to make more quantitative comparisons to empirical v_∞. Inspection of Fig. 9 also reveals that the empirical mass-loss rates and terminal wind speeds show no prominent anti-correlation. This further supports another basic characteristic of our models, namely that Ṁ and v_∞ are quite decoupled (see Sect. 4.3.1).

5.3 Dependence on metallicity

The model of dynamic inflation presented in this paper relies on the fact that the Γ =1 condition is reached already at the Fe-opacity bump. As such it is sensitive to the variation of opacity at the Fe-opacity bump, which is directly controlled by metallicity (e.g. Cantiello et al. 2009). Qualitatively, the model of dynamic inflation should respond to a decrease or increase in metallicity by reducing or enlarging the scale of inflation and the associated mass loss. Moreover, at low enough metallicity the Γ = 1 conditionwill not be fulfilled at the Fe bump, in which case dynamic inflation will no longer take place from the hot Fe bump, and a standard optically thin wind develops instead. This behaviour is also in line with the recently published models by Sander & Vink (2020).

5.4 Comparison to co-moving frame models

As mentioned above, recently Sander et al. (2020) and Sander & Vink (2020) published time-independent wind simulations of classical WR stars. These models solve the same steady-state e.o.m. as here⁸, however using a different method to calculate the radiative acceleration. There Γ(r) is derived from full solutions of the frequency-dependent radiative transfer equation in the co-moving frame (CMF), without the use of any parameterisation to estimate the line-force contribution (see also Sundqvist et al. 2019; Björklund et al. 2021, for similar approaches applied to O-stars). Although these models thus compute Γ(r) for a given velocity and density structure in a more detailed way than here, this CMF method (as currently implemented) has one major disadvantage. Because of the way the CMF transfer equation is numerically solved, it can only be applied to monotonically increasing (or decreasing) velocity fields v(r). As such, this method cannot be used to derive Γ(r) for the type of non-monotonic velocity structures found here (see Figs. 3b, 6b).

The way Sander et al. (2020) overcome this issue is by assuming an ad hoc very high degree of clumping in the WR outflow (see also Gräfener & Hamann 2005). They parameterise this using a simple model where all wind mass is assumed to be concentrated within overdense clumps, such that the densities of these clumps are ρ_cl = ⟨ρ⟩∕f_vol = D⟨ρ⟩, where ⟨ρ⟩ is the mean density (assumed to be preserved with respect to a smooth model), f_vol the volume fraction of the total wind contained by clumps, and D the clump-overdensity factor.

By adopting a very high D = 50 (f_vol = 0.02), the ionisation balance of the wind is shifted so that more effective driving lines become available, and the line force is thereby increased enough to avoid any deceleration regions. When assuming instead a more modest D = 10 in their models, they do experience the same type ofΓ < 1 regions as here, which leads then to a non-monotonic v(r) that cannot be handled by the CMF radiative transfer (see Fig. 1 in Sander et al. 2020, and their corresponding discussion). Since observations of electron scattering wings in WR outflows suggest quite modest values D ≈ 4−10 (see overviews by Crowther 2007; Puls et al. 2008), and theoretical wind-clumping models following the line-driven instability (Owocki et al. 1988) have so far only been developed for OB-stars (with again typical values D ≈ 10; Sundqvist et al. 2018; Driessen et al. 2019), it is at present not clear how the high clump densities needed to keep v(r) monotonic in CMF-based WR models might be physically justified.

The models presented here are thus very complementary to such CMF-based simulations. Although our more simplified treatment of Γ involves certain approximations (as discussed in previous sections) it allows us to model regions of stagnation and deceleration (see Sects. 3 and 4), whereas such non-monotonic flows are incompatible with the methods currently used in CMF-based calculations. Moreover, the method presented here allows for time-dependent simulations that otherwise would not be computationally feasible. In this respect, we also recall our finding that if the line force in the outer wind is decreased too much, the wind not only experiences regions of deceleration but may eventually reach zero velocity. In the 1D simulations presented here, this leads to failed wind solutions and model termination. In a multi-dimensional simulation, one may speculate that it might instead result in a very complex wind structure of co-existing regions of upflows and inflows. In turn, this may then cause shocks and perhaps increased levels of clump formation. These could be simulated in multi-dimensional hydrodynamics simulations using the method we present here.

The models presented here also demonstrate that even with multiple extended regions where Γ <1, the accumulated net energy of the wind can still be positive (i.e., we can still have ). This suggests that the WR outflow can be driven without invoking such a high degree of clumping, but instead allowing for non-monotonic velocities. For example, let us consider the lower panel of Fig. 1 in Sander et al. (2020), where the authors provide Γ(r) not only for their final D = 50 model (their solid line) but also for their collapsed D = 10 solution (their dashed line). From these profiles, we estimated the net energy corresponding to both models and found that, indeed, both of the force configurations provide . As such, both solutions may actually escape the stellar potential, but since the D = 10 model would have contained regions of negative acceleration it could not be treated by the assumed CMF radiative transfer framework.

Note that Sander et al. (2020) make the argument that although the assumed very high degree of clumping influences v_∞ significantly, it does not have a major impact on Ṁ, since clumping is assumed to be present only at supersonic velocities. This is quite consistent with our finding here. On the other hand, Sander & Vink (2020) also find a correlation between Ṁ and v_∞, which we do not necessarily observe in our models. This result, however, might simply be related to them forcing a monotonic velocity field in their simulations, since the same correlation was also found in the β-law WR models by Gräfener et al. (2017).

6 Summary and future work

This paper has presented a model for the wind outflows of classical Wolf-Rayet stars, based on a hybrid opacity method where OPAL opacities are used in combination with a CAK-like line force parametrisation. The latter is used here to account for the effects of Doppler shifts upon spectral line opacity in an approximate way. Our hybrid opacity model essentially results in a two stage WR stellar wind: a supersonic outflow is first initiated in optically thick layers (at the so-called Fe bump) by the OPAL opacities. However, these opacities are not able to prevent fallback upon the core, and so it is the CAK-like force that takes over the driving in the outer regions and ensures that the initiated outflow can also escape the stellar potential.

Our approach leads to deep-seated wind initiation with high mass-loss rates and terminal velocities (see Table 5), in the range of values typically inferred for classical WR stars by spectroscopic analyses. The simulations display dense and slow – dynamically inflated – sub-photospheric layers, offering a natural explanation to the so-called core-radius problem of classical WR stars. Direct comparison with hydrostatic computations further shows that neglecting the dynamical terms in the e.o.m. instead leads to drastic inflation of the then assumed static envelope. In this way, we demonstrateexplicitly that the need to invoke huge amounts of ad hoc energy transport in these convectively inefficient regions of stellar structure and evolution calculations is simply an artefact of neglecting the dynamic outer boundary.

The paper here uses a line-force parameterisation that introduces a very simple radial variation to the CAK power-law index α in order to sustain the outer wind outflow. The quite simple nature of this basic set-up means that it should be rather straightforward to extend these first simulations towards exploration of a more complete stellar parameter space for classical WR stars. Since the assumed CAK line-force parameterisation proved to have little effect on the resulting mass loss (as long as we ensure a positive wind energy), such model grids might then offer good alternatives for direct inclusion of mass-loss rates into stellar evolution models.

The assumed treatment of α does influence the outer wind velocity profile and terminal speed. As such, constraining the morphology of the re-accelerating part of the outflow will require more detailed studies. This is a very challenging task, however, since a non-parametrised calculation of the line force in these regions must be able to deal with non-monotonic velocity fields. As such, current CMF radiative transfer methods, like those applied in some alternative WR- (Gräfener & Hamann 2005; Sander et al. 2020; Sander & Vink 2020) and O-type star wind models (Sundqvist et al. 2019; Björklund et al. 2021), will not be sufficient.

The assumed variation of α can also be varied in multi-dimensional extensions of the time-dependent simulations presented here. In such a study we may allow parts of the wind material to fall back into the deep regions of the outflow. While in 1D simulations such fallback eventually always leads to model termination, in a multi-dimensional setting it may instead create a highly structured, turbulent flow. This might then still deposit enough momentum into the plasma that a wind can be sustained, but now with clumps of material perhaps mechanically transporting some of the momentum through collisions.

As the wind velocity field also affects spectral line formation, another important follow up regards analysis of spectral features associated with the non-monotonic velocities found in this paper. In particular, a key question becomes if there are any fundamental differences in the predicted strengths and shapes of strategic spectral lines calculated from our models as compared to those predicted by standard spectroscopic models assuming a monotonic β-type velocity law (e.g. cmfgen, Hillier & Miller 1998). Using the newly developed 3D radiative transfer code by Hennicker et al. (2020), studies of this are already underway and will be presented in an upcoming paper. Finally, a very interesting problem to address concerns stellar structure and evolution, and the manifestation of dynamic inflation in different evolutionary stages. Such a study might be especially important since quite many stars (some already at earlier evolution stages) reach the Eddington limit at opacity bumps in regions that might have insufficient convective energy transport. As demonstrated in Sect. 3, dynamically inflated envelopes should then develop, which in turn could have a significant feedback on the stellar structure as well as on the further evolution of the star.

Acknowledgements

The authors would like to thank J. Puls, T. Shenar and all members of the KUL EQUATION group for fruitful discussion, comments, and suggestions. L.P., J.S., N.K., L.D., A.d.K., L.M., and H.S. acknowledge support from the KU Leuven C1 grant MAESTRO C16/17/007. L.P. further acknowledges the additional support from Shota Rustaveli Georgian National Foundation Grant Project No. FR17-391. J.S. further acknowledges additional support by the Belgian Research Foundation Flanders (FWO) Odysseus program under grant number G0H9218N.

References

All Tables

All Figures

Fig. 1 (a) Colour map of ΓOPAL for the stellar parameters from Table 1. Temperature is given on the abscissa and radiation-to-gas pressure ratio on the ordinate. is identified with a white contour. The colour bar on the right corresponds to the colour coding of the plot and gives the numerical value of ΓOPAL for a given temperature and radiation-to-gas pressure ratio. Different line styles show the radiation-to-gas pressure ratio and temperature structure for the different fixed mass-loss rates. (b) Ratio of kinetic to core potential energy of the corresponding solutions from panel a. The abscissa here shows the radius coordinate x.

In the text

	Fig. 2 Plot that demonstrates values of α (see text) required to reproduce the line force computed from MC radiative transfer calculations by Springmann (1994; asterisks). The solid line then over-plots the piecewise linear function adopted in this paper to describe the corresponding variation of α. The vertical dashed lines correspond to r1 (x1) and r2 (x2).
In the text

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

In the text

Fig. 4 Density profiles of different static and steady-state solutions. On both plots the solid lines show the different cases of hydrostatic solutions and dashed lines show the sub-photospheric part of the steady-state outflow (see Fig. 3e). The different markers then compare hydrostatic solutions for different assumed Thompson scattering Eddington factor Γe (left panel) and for different assumed L∕Lrad ratios for thefixed total luminosity log10(L∕L⊙) = 5.416 (right panel). The numbers at the right end of each density profile correspond to the photospheric effective temperature of the corresponding model in kilo-Kelvins, as defined in the text. To better visualise the innermost parts of the simulations, the right panel displays the logarithm of the scaled radius x = 1 − Rc∕r on the abscissa.

In the text

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

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

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

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

In the text

Fig. 9 Distribution of mass-loss rates and terminal velocities derived from the spectroscopic fitting by Sander et al. (2019) and Hamann et al. (2019) for an observed stellar population of galactic WR stars of nitrogen (WN), carbon (WC), and oxygen sequences (WO). Mass-loss rates and terminal velocities derived in our steady-state (log10 (Ṁ∕ M⊙ yr−1) = −4.83, v∞ = 2200 km s−1) and time-dependent (log10(Ṁ∕ M⊙ yr−1) = −4.78, v∞ = 2063 km s−1) simulations are consistent with the range of mass loss shown on this figure.

In the text