© The Authors 2023

1 Introduction

The mass loss of hot, evolved, massive stars plays a critical role on multiple astrophysical scales: strong stellar winds affect the individual appearance of the stars (e.g., Hamann 1985; de Koter et al. 1997; Hillier et al. 2001; Shenar et al. 2020) and consequently also their ionizing and energetic feedback to the environment (e.g., Smith et al. 2002; Crowther & Hadfield 2006; Hainich et al. 2015; Sander & Vink 2020). Although the timescales for all evolutionary stages beyond the main sequence are comparably short, mass loss in these stages still considerably affects the stellar fates (e.g., Langer et al. 1994; Chieffi & Limongi 2013; Yusof et al. 2022; Moriya & Yoon 2022). In particular for hydrogen-depleted, classical Wolf–Rayet (WR) stars, strong stellar winds provide a major channel for the chemical enrichment of their host environment (e.g., Maeder 1983; Dray et al. 2003; Farmer et al. 2021; Martinet et al. 2022). The strong winds give rise to the emission-line-dominated spectra of WR stars, leaving their imprint even in integrated spectra of whole stellar populations and galaxies (e.g., Conti 1991; Leitherer et al. 1996; Schaerer et al. 1999; Plat et al. 2019). Since the detectability of gravitational waves (Abbott et al. 2016) and along with it the significant amount of black holes (BHs) above 20 M_⊙ (e.g., The LIGO Scientific Collaboration 2021), the interest in a better understanding of the BH-mass limiting WR mass loss as a function of metallicity (Z) has increased even further (e.g., Woosley et al. 2020; Higgins et al. 2021; Vink et al. 2021).

Contrary to their impact, the theoretical understanding of WR-type winds is still rather limited. Despite earlier doubts, partially exacerbated by the high mass-loss rates determined before clumping was incorporated into wind models, the considerations and model efforts of Nugis & Lamers (2002), Gräfener & Hamann (2005) and Vink & de Koter (2005) demonstrated that the winds of WR stars are mainly radiatively driven with iron opacities playing a critical role for the acceleration of the wind and the scaling of the mass-loss rate. Since then, it required a new generation of computers and a considerable update to the modeling techniques to extend these fundamental efforts to a larger parameter space. Only recently did Sander et al. (2020) and Sander & Vink (2020) manage to calculate a larger set of dynamically consistent 1D atmosphere models that were able to predict the winds of classical WR (cWR) stars over a wider parameter space, though still covering two dimensions only. A key ingredient of these models is the detailed calculation of the flux-weighted mean opacity ℌ_F and thus the radiative acceleration a_rad in an expanding environment without requiring the assumption of local thermodynamic equilibrium (LTE). The new generation of computational capabilities has also opened the path toward multidimensional simulations for WR winds (e.g., Poniatowski et al. 2021; Moens et al. 2022). In contrast to the 1D models, these 3D calculations are time-dependent, but so far limited to LTE and very few test cases, making the current insights from 1D and 3D modeling quite complementary.

In this work, we mainly follow up on the work of Sander et al. (2020) and Sander & Vink (2020), using 1D stellar atmosphere models to explore an additional dimension that is very important to determine the properties and strength of WR winds. Given the high computational costs of dynamically consistent atmosphere models, all sequences presented in Sander & Vink (2020) were calculated using a fixed stellar temperature (T_*) defined at a Rosseland continuum optical depth of τ_R,cont = 20. The corresponding radii R_* for the models were thereby given via the Stefan-Boltzmann law $$L=4\pi R_{*}^2{\sigma }_{\text{SB}}{T}_{*}^4.$$(1)

While the value of T_* = 141 kK in Sander & Vink (2020) was well motivated by the prototypical solution for a classical WC star (Gräfener & Hamann 2005), there is a priori no reason to assume that this choice of T_* is valid for all He-burning WR stars. In fact, stellar structure models predict a curvature in the zero age main sequence (ZAMS) for He stars (e.g., Langer 1989; Köhler et al. 2015) with lower temperatures obtained for lower masses. Since we could not take this effect into account in Sander & Vink (2020), we had to limit the applicability of the derived Ṁ recipe to He stars of about 10 M_⊙ and higher.

The radii of WR stars and – as a consequence – also their temperatures are a long-standing topic of active research (e.g., Hillier 1991; Hamann & Gräfener 2004; Grassitelli et al. 2018; Sander et al. 2020). Beside the curvature in the He ZAMS, we are facing a particular challenge for stars with dense winds by the photosphere shifting to highly supersonic velocities. Thereby, the spectral appearance is completely determined in the wind, providing no direct observable (e.g. log g) which could be used to determine the (hydrostatic) stellar radius. This raises the problem of connecting the effective temperatures for a Rosseland optical depth of τ_R = 2/3 (T_2/3), which can be obtained via quantitative spectroscopy, to the (much) deeper subsonic regime represented by T_* for stars with extended envelopes (cf. the discussion in Sander et al. 2020). In principle, the actual hydrostatic radii of the stars could be much larger. This is referred to as (hydrostatic) inflation. Alternatively, a relatively compact star can be cloaked in a wind that is optically thick out to significant radii. Both solutions might actually occur in nature with the realized branch depending on the particular stellar parameters.

Stellar structure calculations can help to get a handle on R_* and T_* for helium stars, but their inherent (and computationally necessary) limitation to gray opacities usually prevents a proper estimation of T_2/3. For a selected evolutionary track of a 60 M_⊙ star, Groh et al. (2014) calculated stellar atmosphere models adopting stellar parameters derived from an evolutionary track. This method led to important revisions on the predicted spectral appearances and their duration during the later evolutionary stages of massive stars. However, despite the more sophisticated method to obtain the improved effective temperatures, their underlying mass-loss rates were taken from a simplified recipe inherent to the evolutionary calculations.

With the calculation of hydrodynamically-consistent atmosphere models, we can now obtain consistent mass-loss rates Ṁ and effective temperatures for He-burning stars without requiring any prescription of Ṁ. In this work, we use this technique to investigate the behavior of T_2/3 and other temperature scales for multiple sequences of models with extended atmospheres. The paper is structured as follows: In Sect. 2, we briefly introduce the model atmosphere code including its recent updates necessary for our study as well as the calculated model sequences. In Sect. 3, we present the resulting temperature trends, starting with an exemplary discussion of one sequence before exploring the full sample. Afterwards, we take a closer look at the obtained trends in the terminal wind velocities in Sect. 4. Evolutionary implications and resulting scaling relations for the imprint of the WR effective temperatures are discussed in Sect. 5. The insights on He II ionizing fluxes are introduced in Sect. 6 before drawing the conclusions in Sect. 7.

2 Stellar atmosphere models

In this work, we employ the PoWR model atmosphere code (Gräfener et al. 2002; Hamann & Gräfener 2003; Sander et al. 2015) in its hydrodynamical branch (PoWR^HD) to calculate stationary, hydrodynamically-consistent atmosphere models. The implementation concepts for coupling hydrodynamics and radiative transfer are described in Sander et al. (2017, 2018, 2020). Our hydrodynamic solutions are calculated in a similar manner as described in Sander & Vink (2020), that is we keep the stellar parameters L and M_* fixed and iteratively adjust Ṁ and υ(r) until a consistent solution is obtained.

In our previous studies, the solutions obtained for the velocity field were always monotonic. In this work, we demonstrate that apart from the high clumping factor (D_∞ = 50), this was mainly a result from choosing T_* = 141 kK as the anchor point of our model sequences.

When extending our modeling efforts to lower T_*, we now approach a regime where Γ_rad = a_rad g⁻¹ can drop below unity after launching the wind. Such a deficiency in the available radiative acceleration is regularly seen in WR atmosphere models including the opacities of the hot iron bump (e.g., Gräfener & Hamann 2005; Aadland et al. 2022). When assuming a prescribed velocity field, as traditionally done in spectral analysis, these deficiency regions have no immediate impact on the model calculations. This is different in our case where we solve hydro-dynamic equation of motion. Here, a region with Γ_rad < 1 in the supersonic regime implies a negative velocity gradient until Γ_rad eventually raises above unity again, resulting in a nonmonotonic velocity υ(r) (see, e.g., the recent calculations and discussions in Poniatowski et al. 2021). Given that our atmosphere modeling technique needs to perform the radiative transfer in a co-moving frame, which cannot handle nonmonotonic velocity fields, we therefore have to modify the υ(r) obtained from the solution of the hydrodynamic equation of motion.

Two approaches are used in this work, which are illustrated in Fig. 1. In the simple, numerically more robust method, we ignore any negative gradients already during the solution of the equation of motion. With this method, we obtain a locally consistent solution at all depth points except for any supersonic regions where Γ_rad < 1. However, this method leads to an over-prediction of υ_∞ as any reduction in υ(r) due to parts with negative gradients is ignored. The example with the dashed-dotted curve in Fig. 1 illustrates that this approach can remove all further structure from the outer velocity field. We thus calculate a second type of models where the integration of the hydrodynamic equation of motion is not perturbed and only the resulting velocity field is interpolated afterwards. For the latter interpolation, we use an “outside-in” approach, starting at the outer boundary of our model (R_max) and cutting away any parts where υ(r) increases inward. While the such modified solution actually leads to more local violations of the hydrodynamic equation of motion (cf. Fig. 2), we preserve not only the correct υ_∞ but usually also the whole υ(r) in the optically thin regime as illustrated with the red-dashed curve in Fig. 1. We therefore use these types of models as the main anchor-point for discussing our results and drawing conclusions.

To be able to compare our results to the large set of calculations performed in Sander & Vink (2020), we keep most of our original model input, including the clumping description (D_∞ = 50 and υ_cl = 100 km s⁻¹, using the “Hillier law” from Hillier & Miller 1999) and the set of considered elements (cf. Table 1). However, we have calculated some additional models, including two complete T_* sequences for 12.9 M_⊙ and 20 M_⊙, with D_∞ = 10 as well as one sequence with D_∞ = 4 to have a comparison sets which turns out to be quite insightful.

In Fig. 3, we display the resulting contributions to the radiative acceleration from two models which only differ in D_∞. The higher D_∞ changes the wind stratification, most notably by stronger recombination from He III to He II (indicated by the higher He II bound-free opacity bump) and an earlier ionization change from Fe VI to Fe V in the outer wind, enabling additional line driving from Fe IV. (A more detailed breakdown with all ionic contributions is provided in Appendix C with Figs. C.1 and C.2 and brief discussion about the comparison.)

Furthermore, we calculate a few additional sequences where we add surface hydrogen (X_H = 0.2) and one sequence where we switch to a WC-type (X_C = 0.4, X_N = 4 × 10⁻⁵, X_O = 0.05) metal composition. A full overview of the model sequences is given in Table 2. Similar to Sander & Vink (2020), we again align L and M_* such that they follow the relations for hydrogen-free stars given in Gräfener et al. (2011). This means that we ignore any potential extra mass (and luminosity) due to surface hydrogen as well as any differences in the L-M relation between WN and WC stars. We do so in order to isolate the effects of different chemical compositions on the resulting wind predictions rather than trying to emulate an observed star or a particular evolutionary model. A more detailed investigation of the impact of (surface) hydrogen on the mass loss of WR stars is currently underway and will follow in a separate paper.

While not tailored to mimic any particular WR star, the spectra resulting from our model sequences display a typical WR-type appearance. As an example, we plot four spectra from the 20 M_⊙ WN sequence with X_H = 0.2, Z = Z_⊙ and D_∞ = 10 in Fig. 4. The hottest model shown mimics typical features of a relative weak-lined, early-type WN with intrinsic absorption lines, e.g., seen in WN3ha stars. Along the sequence, the lines tend to get stronger, but narrower with additional lines from cooler ionization stages appearing in the cooler models that would be classified as later WN types. Fine-tuned model efforts for particular stars will be necessary to further constrain choices of currently free parameters such as microturbulence and clumping. Ideally, one would want to eliminate the necessity for a dedicated input of these parameters completely to get full dynamical consistency, but this would require significant code updates, which is beyond the scope of the present paper.

Fig. 1 Illustrating example of the possible velocity treatments in case that the integration of the hydrodynamic equation of motions yields a nonmonotonic υ(r) (solid blue curve): in the simple treatment, negative velocity gradients are suppressed during the integration, yielding the blue, dash-dotted curve. In some cases, such as illustrated in this plot, this can spoil the terminal velocity υ∞. In the more sophisticated treatment, negative gradients are therefore taken into account and υ(r) is modified such that the interpolated solution keeps the obtained υ∞.

Fig. 2 Resulting acceleration stratification of the converged models with two different treatments of nonmonotonic velocity fields: The dashed red and black lines illustrate the two sides of the hydrodynamic equation of motion for a model where negative velocity gradients are ignored in the solution. The corresponding solid lines show the result for the alternative method where the nonmonotonic υ(r) is interpolated afterwards. The small inlet shows the resulting velocity fields for both models. In this example we show models with T* = 125 kK, log L/L⊙ = 5.475 and M = 15 M⊙.

Fig. 3 Major contributions to the radiative acceleration for two hydrodynamically consistent, hydrogen-free WN models with T* = 130 kK, log L/L⊙ = 5.7, M = 20 M⊙, and D∞ = 10 (upper panel) or D∞ = 50 (lower panel): for the line contributions, all elemental contributions except Fe are summed over all ions. The total radiative acceleration (arad), the Thomson acceleration from free electrons (aThom = Γe · g), and the contribution from gas (and turbulence) pressure (apress) are also shown for comparison. The loosely dashed horizontal line denotes the total Eddington limit that needs to be overcome to launch a wind.

Fig. 4 Example spectra from our WN model sequence with log L/L⊙ = 5.7 and M = 20 M⊙, XH = 0.2, Z = Z⊙, and D∞ = 10 for the optical wavelength regime with different colors corresponding to the different models.

3 Temperatures and radii

For our sequences listed in Table 2, we study the mass-loss rate as a function of the effective temperature at the critical (≈ sonic) point T_eff(τ_crit). The latter value is very close to T_* defined at τ_R,cont = 20, which usually describes the inner boundary of our models, except for models with very high Ṁ where τ_R,cont = 100 needs to be chosen as the inner boundary. The results depicted in Fig. 5 reveal a steep, monotonie decrease of Ṁ with increasing value of T_eff(τ_crit) (corresponding to decreasing radii R_crit). This is not unexpected given that larger radii lower the local gravitational acceleration and thus enable an easier escape of material.

Our sequences with different chemical compositions give us first qualitative insights on the impact of hydrogen on the one hand and carbon and oxygen on the other hand: The direct comparison of the two curves for 20 M_⊙ at Z_⊙ in Fig. 5 show a systematic shift to slightly higher mass-loss rates in the presence of surface hydrogen (as long as it is negligible for the total stellar mass). Interestingly, a hydrogen surface mass fraction of X_H = 0.2 seems to be sufficient to counter the lower metal abundances in the 0.5 Z_⊙ sequence. This result should not yet be generalized given the limited number of sequences, but will be followed up in our dedicated study focusing on surface hydrogen. Contrary to hydrogen, the inclusion of more carbon and oxygen is not beneficial to Ṁ as illustrated in Fig. 5 by the sequence with the WC surface composition. In all cases, the local (electron) temperatures at the launching point of the wind (R_crit) are too high to generate any additional line opacity from С or O. Instead, the slightly lower amount of free electrons compared to a WN surface composition decreases the resulting Ṁ (cf. Sander et al. 2020). The opposite effect instead happens in the case of WN stars with X_H > 0.

For comparison, Fig. 5 also contains empirical results obtained with standard models using a β-law or double-β-law to describe υ(r) from Hamann et al. (2006, 2019); Hainich et al. (2015); Shenar et al. (2016, 2019). Illustrating the well-known “Wolf–Rayet radius problem” (e.g., Grassitelli et al. 2018), most empirically derived temperatures are located at T_*, values that seem to be too cool for their mass-loss rate. As demonstrated by Gräfener & Hamann (2005) and Sander et al. (2020), the critical radii inferred from dynamically-consistent atmosphere models are much smaller than those obtained by a standard β-law due to the opacities of the “hot iron bump” that enable the launch of a supersonic wind already at deeper layers. Moreover, the comparison in the Ṁ-T_eff(τ_crit) plane is not ideal as the empirically derived values of Ṁ depend on distances which are still uncertain for some Galactic targets, but as we see below when discussing the transformed mass-loss rate, this is not a major issue here. An inspection of the model sequences indicates a clear shift for different mass regimes. Thus, one could argue that the whole plane in Fig. 5 could be covered if we would calculate further sequences for lower masses. However, lower masses are most likely not the solution and discrepancies remain even when adjusting the mass-loss rates for stars with different luminosities in Appendix D.

Fig. 5 Mass-loss rate Ṁ as a function of Teff(τcrit) for our model sequences. For comparison, a set of empirically inferred temperatures T* for different types of WR stars is shown as well (gray symbols), illustrating the well-known “WR radius problem.” Since the empirical values are inferred from models without dynamical consistency, the values of T* defined at a Rosseland optical depth of τR,cont = 20, are shown. For our model sequences, the values of T* and Teff(τcrit) align very closely except for the highest mass-loss rates. A direct comparison with T* from the dynamically-consistent models is provided in Fig. E.1.

Fig. 6 Mass-loss rates versus different temperature scales for a series of dynamically consistent atmosphere models with log L/L⊙ = 5.7, M = 20 M⊙, and ХH = 0: the thick red dashed line denotes the effective temperatures defined at a Rosseland optical depth of τR = 2/3, while the green solid line and the blue dashed-dotted line denote the effective temperatures referring to τcrit and τR,cort = 20, respectively. The green dotted line on the right denotes the (electron) temperature at the critical point. Curves in lighter colors reflect models using the simple integration treatment suppressing negative velocity gradients (cf. Sect. 2). The black dashed line shows the hydrodynamic structure solutions by Grassitelli et al. (2018).

3.1 Mass loss versus different temperature scales

To discuss the different temperature scales in WR winds and their scaling with Ṁ, we take a closer look at an individual model sequence. In Fig. 6, we plot the mass-loss rate Ṁ for the sequence of 20 M_⊙ models without hydrogen as functions of different temperature definitions, namely (i) the effective temperature T_2/з at a Rosseland optical depth of τ_R = 2/3 (thick red dashed line), often simply denoted as T_eff in the literature; (ii) the effective temperature T_* commonly used in the model setup for PoWR models, defined at a Rosseland continuum optical depth of τ_R,cont = 20; (iii) the effective temperature at the critical point, denoted T_eff(τ_crit), T_eff(R_crìt), or simply T_eff,crit; and (iv) the electron temperature T_e at the critical point (thin green dotted line). In contrast to the first three temperatures, T_e is not an effective temperature. In the deeper layers of the atmosphere where the deviations from LTE become negligible, T_e aligns with the general temperature T(r) defined in (1D) stellar structure models. To further illustrate the effect of the two different nonmonotonic υ(r)-treatments, we plot the more accurate method with posterior interpolation of the velocity field in strong colors while the simple method ignoring negative gradients is drawn in lighter shades of the same line style. Beyond a certain temperature, there are no more deceleration regions and thus both curves agree.

Except for the regimes with highest mass-loss (Ṁ > 10⁻⁴ M_⊙ yr⁻¹ in Fig. 6), the values of T_* and T_eff,crit closely align. This does not imply that τ_crit has to correspond to τ_R,cont ≈ 20, but that the locations of the radii corresponding to τ_R,cont = 20 and τ_crit are close enough to yield similar effective temperatures. The alignment between T_* and T_eff,crit at lower Ṁ is fulfilled in all of our model sequences (see Figs. 7 and E.3 for further examples) and allows us to discuss the physically more meaningful temperatures at the critical point (i.e., the launch of the wind) instead of the slightly more technical T_*. For higher Ṁ, τ_crit moves inward and eventually surpasses τ_R,cont = 20, explaining the deviation of the curves for the highest mass-loss rates. However, these cases require models very close to the Eddington limit.

The growing difference between the effective temperature at the launch of the wind (T_eff,crit) and T_2/з with increasing Ṁ shows the “extended atmosphere” of a WR star. For low mass-loss rates, the atmosphere is optically thin and the two temperatures align. Although at usually much lower temperatures, this is similar to what we see for most OB-star winds. With increasing Ṁ, we get a more and more extended optically thick layer. Albeit leading to a much cooler appearance of the star, this kind of layer should not be mixed up with the inflated envelope obtained in various hydrostatic structure models (e.g., Petrovic et al. 2006; Gräfener et al. 2012; Ro & Matzner 2016). Instead of a subsonic, but still loosely bound extended layer, our models show supersonic (i.e., unbound) layers moving out with hundreds of km s⁻¹. As illustrated in Sander et al. (2020), the winds often reach more than 0.5 υ_∞ before the atmosphere becomes optically thin. In more recent work (e.g., Poniatowski et al. 2021), this form of an extended photosphere is termed “dynamical inflation” to distinguish it from the hydrostatic inflation. Hydrostatic inflation is not likely to occur if a wind can be launched and maintained, but it could in situations where the latter is not given. In Sect. 3.2, we discuss the limits of launching a wind from the hot iron bump. While a detailed exploration of wind solutions beyond this limit is not feasible in this study, the numerical results of our “failed” models indicate a tendency toward larger sonic radii, potentially indicating some form of hydrostatic inflation.

Finally, we also plot the electron temperature at the critical (≈ sonic) point as a dotted green curve in Fig. 6. The values reflect the expected temperature range of the hot iron bump (around 200 kK), although the particular values are slightly higher than predicted in structural studies employing OPAL opacity tables (Grassitelli et al. 2018; Nakauchi & Saio 2018). The value of T_е(R_crit) appears relatively constant at first sight, but aside from numerical scatter affecting the results a bit, a subtle trend can be noticed: In the regime of optically thick winds, there is a tendency toward increasing T_е(R_crit) with higher mass-loss rates. This is qualitatively in line with the predictions by Grassitelli et al. (2018), who found that in hydrodynamic stellar structure calculations higher mass-loss rates correspond to higher temperatures at the sonic point. However, this trend is interrupted when the winds become optically more thin and eventually T_е(R_crit) increases mildly with lower Ṁ until T_eff,crit surpasses T_е(R_crit).

All of the temperature trends described for the exemplary Fig. 6 are observed for the other model sequences as well. For comparison, we show similar temperature scale plots for the 15 M_⊙ sequence (Fig. 7) and the 12.9 M_⊙ sequence with surface hydrogen (Fig. E.3). While there are shifts in the absolute values, the same general trends are clearly identified.

To study whether our findings are more general or limited to our new sample, we check the behavior of the different temperature scales also for the whole set of model sequences from Sander & Vink (2020). The resulting curves are presented in Fig. 8. Due to the fixed value of T_*(τ_R,cont = 20) = 141 kK in Sander & Vink (2020) and the launching of the winds at high optical depths, the effective temperature referring to the critical point (i.e., the launch of the wind) hardly varies over the whole sample. Still, the resulting T_2/з-temperatures look very similar to those obtained in our new models with varying T_*. On the other hand, T_е(R_crit) varies much more than in any of our new model sequences. As we also see a shift in the T_е(R_crit)-curves between different mass sequences in our new work, we can conclude that Г_е – defined by the chemical composition and L/M – plays a major role in setting the temperature regime of the sonic point. The ratio between the flux and the radius – which is mapped in T_* – instead only has a minor effect. We do not see the interruption of the T_е(R_crit)-trend in Fig. 8 that was apparent in Fig. 6 and the other new model sequences. This is likely due to the different dimensionality of the sequences in Sander & Vink (2020) (fixed T_*, variable L/M per sequence) and this work (fixed L/M, variable T_*, per sequence). In general, we can conclude that for stars further away from the Eddington Limit, the same Ṁ can only be reached by shifting the critical point to lower electron temperatures. For the same L/M-ratio, however, we see a much lower amplitude of changes in T_е(R_crit). In a zeroth-order approximation, one could state that T_е(R_crit) is constant for a given L/M and chemical composition.

Fig. 7 Same as Fig. 6, but for a model sequence with log L/L⊙ = 5.475, M = 15 M⊙, and ХH = 0. Contrary to the situation in Fig. 6 for a 20 M⊙ star, there is an abrupt breakdown of solutions beyond a minimum temperature and a maximum Ṁ. For the T2/з- and Teff(τcrit)-scales these points are marked with a vertical line attached to a gray-hatched area.

Fig. 8 Mass-loss rates as a function of different temperature scales for the L/M model sequences presented in Sander & Vink (2020). Higher mass-loss rates correspond to higher L/M-ratios in this plot.

3.2 WR-type mass loss and its breakdown

The trend of increasing Ṁ with lower T_eff,crit does not automatically continue beyond the plotted values. The comparison between the simpler and the more sophisticated treatment in Fig. 7 already suggests that the effect of deceleration regions has to be taken into account for computing a more realistic Ṁ. In some situations, such as the one illustrated in Fig. 7, the deceleration region can become large enough to reduce the wind to subsonic or even negative velocities, making it impossible to launch a wind from the deeper layers of the “hot iron bump.” This regime occurs right next to the (theoretical) maximum of Ṁ along the T_eff,crit-axis which is reached when the deceleration region is just not strong enough to put υ(r) below the local sound speed in the wind.

The situation of a failed wind launch is illustrated in Fig. 9, where the radiative acceleration barely reaches Γ_rad = 1 in a model for 10 M_⊙. This example also illustrates that for lower L/M values, the regime where no wind can be launched from the hot iron bump gets larger and larger. A hydrogen-rich surface can compensate this to some degree as it helps to get the star closer to the Eddington limit. Still, when getting to lower and lower L/M-ratios the regime of WR winds driven by the hot iron bump eventually vanishes. In our models with a fixed L-M-relation this corresponds to a limit in both luminosity and mass. However, objects of lower masses and luminosity can potentially launch a wind if they have a considerably higher L/M-ratios than homogeneous He stars. This might e.g. be the case in WR-type central stars of planetary nebulae. Gräfener et al. (2017) and Ro (2019) also pointed out that most of the H-free WN population in the LMC presents a challenge as these stars should not be able to launch a wind from the hot iron bump if their masses would adhere to a typical L–M relation for He-burning stars. However, this discrepancy is already reduced if we use the Sander & Vink (2020) models, likely due to the computed flux-weighted opacities exceeding the OPAL Rosseland opacities assumed as a proxy for ℌ_F in previous studies. The discrepancy could potentially be reduced even further slightly lower temperatures are considered as well (cf. Table 3). Nonetheless, the LMC sample remains an interesting test-bed for detailed comparisons with individual objects and the limits of radiation-driven winds from the hot iron bump.

Summarizing the limits of Ṁ along the temperature axis, we see two very different behaviors: Toward cooler temperatures, we have an abrupt breakdown of the thick wind regime when the effect of the deceleration region outweighs the initial acceleration by the hot iron bump. This endpoint is reached close to the highest possible mass-loss rate (for the given stellar parameters) in this whole wind regime. On the hot temperature end, we instead proceed rather smoothly into the regime of optically thin winds with lower and lower values of Ṁ. This drop along the T-axis is significantly shallower than the strong breakdown of Ṁ along the L/M-axis we obtained in Sander & Vink (2020).

Fig. 9 Illustration of the radiative acceleration for a model with 10 M⊙ and T* ≈ 115 kK which is not capable of launching a wind from the “hot iron bump” and thus is dynamically not converged: In the inner part, the total radiative acceleration (red dashed line) approaches Γrad = 1, but does not surpass it sufficiently to launch a wind that could be maintained in the following deceleration region.

3.3 Scaling with the transformed mass-loss rate

Since the different calculated model sequences show very similar slopes for Ṁ(T_2/з), we investigate whether there is a common scaling behind these curves. Given the empirical scaling relations for WR spectra and our findings from Sander & Vink (2020), we study the “transformed mass-loss rate” $${\dot{M}}_{\text{t}}=\dot{M}\sqrt{D}\cdot \left(\frac{1000\text{km\hspace{0.17em}}{\text{s}}^{-1}}{{\upsilon }_{\infty }}\right){\left(\frac{{10}^6{L}_{\odot }}{L}\right)}^{3/4},$$(2)

originally introduced by Gräfener & Vink (2013), for our new model sequences as a function of T_2/з. As depicted in Fig. 10, the resulting curves align extremely well when plotting Ṁ_t instead of Ṁ. Major offsets are only introduced when assuming different (maximum) clumping factors D_∞. Given our findings in Sander et al. (2020) and Sander & Vink (2020), the latter is not much of a surprise. While the clumping does not directly affect the radiative transfer, the solution of the statistical equations are solved for a higher density D · ρ (Hamann & Koesterke 1998). This affects the ionization stratification and usually leads to a larger opacity and thus larger terminal velocity υ_∞. While the mass-loss rate Ṁ is typically not much affected, the resulting Ṁ_t is changed due the increase in υ_∞ being smaller than the increase in $\sqrt{D_{\infty }}$. For our 20 M_⊙ models with X_H = 0.2, the typical increase was about 20% in υ_∞ when increasing D_∞ from 4 to 10 and about 40% when increasing D_∞ from 10 to 50.

To quantify our finding, we flip the axes and show a double-logarithmic plot in Fig. 11. A clear transition between two regimes is evident with a “kink” around log Ṁ_t ≈ −4.5 that appears to be independent of D_∞. The more dense wind regime (T_2/з < 130kK and log Ṁ_t > −4.5) can be reasonably well approximated by a linear fit, yielding $$\log \frac{T_{2/3}}{\text{K}}=\left(-0.49\pm 0.01\right)log\frac{{\dot{M}}_{\text{t}}}{M_{\odot }{\text{yr}}^{-1}}+\left(2.91\pm 0.02\right)$$(3)

for the sequences using D_∞ = 50. Given the inherent numerical scatter, in particular in υ_∞ entering Ṁ_t, we can conclude that in the limit of dense winds $T_{2/3}\propto {\dot{M}}_{\text{t}}^{-1/2}$.

When comparing the relations with empirically obtained values of WN (Hamann et al. 2006, 2019; Hainich et al. 2015; Shenar et al. 2016, 2019) and WC stars (2012, 2019; Aadland et al. 2022), it is immediately evident that comparing T_2/з between empirical and theoretical results yields a much better match than the comparison between the empirical T_*, and our theoretical Т_eff(τ_crit) in Fig. 5 (or the direct T_*-comparison in Fig. E.1). While the mismatch in Fig. 5 illustrates the “Wolf-Rayet radius problem” discussed at the beginning of Sect. 3, the better alignment of the T_2/з values underlines that value of the empirical analysis, despite the dynamical concerns. In empirical studies with fixed velocity fields, models are chosen such that they reproduce the observed spectrum. Although τ_2/3 has a significant wavelength dependence in WR winds, the effective temperature corresponding to the Rosseland mean value provides some form of a representative value for the regime that needs to be met when the light eventually escapes from the star. Our dynamically consistent models can generally reproduce these T_2/з values, but employing more compact radii that better align with structural predictions.

Despite the generally better match when comparing T_2/3, it is also evident from Fig. 11 that all symbols are either on or leftward of the derived curve for D_∞ = 50. The most striking discrepancies are obtained for the SMC WN stars. The Aadland et al. (2022) WC and WO results are very close to our obtained relation. As they are the only one assuming D_∞ = 20, some discrepancies are likely rooted in different clumping assumptions and treatments. The mismatch of the empirical SMC positions however, cannot be explained with clumping differences alone with most of the stars showing empirical Ṁ_t values that are about an order of magnitude lower than our model relations. This could be due to various effects including considerable differences in L/M, e.g., due to having significant hydrogen shells and thus not obeying the assumed L-M-relation in our model sequences, or too low T_2/з estimates. Investigating these and other possibilities would add further dimensions to our model sequences and thus we have to postpone a dedicated analysis of individual targets to a separate follow-up paper.

When considering the obtained curves in the T_2/3-Ṁ_t-plane from the current model sequences, it is so far unclear whether the slope and even the underlying scaling is universal. In Fig. 12, we thus plot the same parameters, now using the sequences from Sander & Vink (2020). A first noticeable difference is the upper horizontal cutoff at T_2/3 ≈ 141 kK, but this is expected due to the fixed value of T_* = 141 kK in the Sander & Vink (2020) sample. The behavior at higher values of Ṁ_t, looks similar to Fig. 11 at first – although with considerably more scatter – but an actual fit of the data reveals a non-negligible difference in the slopes, yielding $$\log \frac{T_{2/3}}{\text{K}}=\left(-0.667\pm 0.009\right)\text{log}\frac{{\dot{M}}_{\text{t}}}{M_{\odot }{\text{yr}}^{-1}}+\left(2.21\pm 0.03\right)$$(4)

for the Sander & Vink (2020) sample. (For comparison, the derived trend for the new sequences is shown as well in Fig. 12.) We discuss possible origins later in Appendix D.

Fig. 10 Transformed mass-loss rate as a function of T2/з for our model sequences.

Fig. 11 Effective temperature at a Rosseland optical depth of 2/3 as a function of the transformed mass-loss rate Ṁt for our new calculated model sequences. The sequences connected by solid lines all employ D∞ = 50, while those with dashed and dotted curves indicate sequences using D∞ = 10 and 4, respectively, as indicated in the plot. For comparison, also various empirical results from the literature are depicted by discrete, gray symbols.

Fig. 12 Effective temperature at a Rosseland optical depth of 2/3 as a function of the transformed mass-loss rate Ṁt for the whole set of models from Sander & Vink (2020). For log (Ṁt [M⊙ yr−1]) < −5.5, T2/3 is effectively independent of Ṁt. In Sander & Vink (2020), the value of T* is fixed for all models, but the difference in L/M still yields a wide range of T2/з values. The dashed-dotted line represents a linear fit of the temperature trend for log (Ṁt [M⊙ yr−1]) > −4.5 while the dotted curve represents the fit for the new model sequence illustrated in Fig. 11.

3.4 Quantitative mass-loss radius-dependence

The similarity of the curves in Fig. 5 and the simplicity of the slopes, indicates a common dependence between log Ṁ and log T_eff,crit for our sequences with fixed L/M. Performing a linear fit, we obtain a relation in the form of $$\log \left(\dot{M}\left[M_{\odot }{\text{yr}}^{-1}\right]\right)=-6\cdot \log \left(T_{\text{eff,crit}}\left[\text{K}\right]\right)+\text{offset}$$(5)

with the detailed fit coefficients being presented in Appendix A and Table A.1. The factor −6 in Eq. (5) implies that the obtained temperature dependence essentially reflects the radius change of the stellar models, since $T_{\text{eff,crit}}^4\propto R_{\text{crit}}^{-2}$ and $$\dot{M}\propto R_{\text{crit}}^3$$(6)

for models with L = const, (as in our model sequences). Our corresponding plot (Fig. A.1) further shows that deviations from the purely geometrical trend occur when we reach the limit of radiative driving discussed in Sect. 3.2 (and listed explicitly for each sequence in Table A.1), as e.g. visible at the upper end of the WC sequence (cf. Fig. A.1).

From a dynamical perspective, the obtained R_crit-trend for Ṁ can be understood as a dependence on the gravitational acceleration on the critical point, where the wind is launched. With the straight-forward definition of $$g_{\text{crit}}=g\left(R_{\text{crit}}\right)=\frac{GM}{R_{\text{crit}}^2}$$(7)

we can rewrite Eq. (6) as $$\dot{M}\propto g_{\text{crit}}^{-3/2}$$(8)

since M is a constant among each of the model sequences. Trend curves reflecting Eq. (8) are displayed in Fig. 13 together with the curves from our model sequences. Generally, a decreasing trend of Ṁ with log g_crit is not surprising as an increased gravitational force needs to be overcome. Given the constant mass M along the model sequences, the change in g_crit expected from Eq. (8) is purely geometrical, that is only from the change in R_crit. The model sequences align well with Eq. (8), but there is a notable flattening for the highest mass-loss rate, that is in the case of more dense winds. We cannot rule out that a numerical effect is playing a role here as these high-Ṁ models often operate on the limits of what the code is capable of. Nonetheless, given that the bending occurs in all sequences, a physical origin seems more likely and we continue our efforts on this assumptions. For some sequences, there are also notable deviations from Eq. (8) at the lower Ṁ-end. From the current set of calculations, the apparent kink in the curves approximately coincides with T_e(R_crit) surpassing T_eff,crit, meaning that the electron temperature at the critical point is higher than the effective temperature at this point. However, the low number of models where this trend is clearly observed and the need to include higher ionization stages in these models, which can cause an additional offset in the numerical solutions if not done early enough in the sequence, currently refrain us from concluding whether there is a clear “kink” or a more gradual change that might potentially be emphasized by a switch in the numerical setup. In any case, the model solutions obtained in this thinner wind regime are characterized by (electron) temperature stratification that remain very high, e.g. larger than 50 kK, until infinity. Their leading acceleration is provided by Fe M-shell ions, which are populated throughout the wind, qualitatively similar to the example shown in Fig. 16 of Sander et al. (2020). When considering the transformed mass-loss rate Ṁ_t instead of Ṁ, the scaling of the velocity with g_crit has to be considered as well, which we do in the following section.

Fig. 13 Mass-loss rate Ṁ as a function of the gravitational acceleration at the critical radius $g_{\text{crit}}=GM{R}_{\text{crit}}^{-2}$. Thin, dotted, gray lines indicate curves with $\dot{M}\propto g_{\text{crit}}^{-3/2}$.

4 Terminal velocity trends

With the intrinsic solution of the hydrodynamic equation of motion, our models automatically predict terminal wind velocities together with Ṁ. While already entering the transformed mass-loss rates, we now take a look at the explicit results for υ_∞ as a function of T_2/з in Fig. 14. Although our sequences are not at all adjusted to match any particular observations, we also plot empirical results for WN stars obtained with PoWR for comparison. It is clear that our models match the general regime of the observed sample, but a closer inspection also shows caveats, for example with hydrogen-free WN stars showing values above the 20 M_⊙ H-free sequence. Various possibilities could explain this (e.g., higher L/M and/or higher clumping), but a thorough investigation is beyond the scope of this paper.

4.1 Impact of clumping

In contrast to Ṁ, the values for υ_∞ tend to scatter a bit more due to being evaluated at the outer boundary of the models. The terminal velocity further strongly depends on the included opacity, so including all ions contributing to the acceleration is necessary in order to avoid underestimating υ_∞. As apparent from Fig. 14, υ_∞ also reacts on the choice of the clumping factor D_∞. With a depth-dependent onset of the clumping, the response of the mass-loss rate to a change of D_∞ is usually small as the resulting differences in D(r) are small in the subsonic layers (see also Fig. 1 in Sander et al. 2020, and Appendix С of this work). In the supersonic layers, however, any differences in D_∞ affect the bulk of the opacities being considered in the hydro-dynamic equation of motion. Consequently, the obtained values for υ_∞ are notably higher for higher D_∞, typically on the order of 20% when calculating a model for the same L, M, T_*, and chemical composition with D_∞ = 50 instead of 10. In particular cases, the effects can be much larger, e.g. when a model has a deceleration regime in the case of a lower D_∞, while there is no such regime for higher D_∞. Significant changes of the wind density regime due to a switch of D_∞ can then lead to a stronger change in the derived Ṁ. Moreover, the assumption of little to no clumping in the deeper layers would become invalid in case of a porous, optically thick medium, which can result from a subsonic, super-Eddington situation (e.g., Shaviv 1998, 2000).

Fig. 14 Terminal velocity as a function of T2/3 for the different model sets. For comparison, also the derived trend for OB-star winds in the Milky Way (dashed line) and the LMC (dashed-dotted line) from Hawcroft et al. (in prep.) are shown.

4.2 Scaling with T_2/3

Beside the differences due to clumping, Fig. 14 demonstrates that any significant change of the chemical composition usually affects the derived υ_∞-values. In the figure, we see higher terminal velocities for the 20 M_⊙ model sequence with surface hydrogen (X_H = 0.2) compared to the corresponding hydrogen-free sequence. This result might be counter-intuitive at first as hydrogen does not provide significant line opacity that could be used to increase υ_∞. However, the additional hydrogen is able to boost the mass-loss rate of the star as the hydrogen atoms in the atmosphere provide a higher budget of free electrons compared to a hydrogen-free atmosphere¹. With more acceleration available already in the deeper layers, the critical point of the wind moves inward to higher optical depths. Line opacities which were subsonic in the hydrogen-free case can now be used to further boost the terminal wind speed. Due to the higher mass loss of the hydrogen-containing model, the value of T_2/3 decreases when comparing models with the same T_*. Interestingly, as we saw in Fig. E.1, the value of ${\dot{M}}_{\text{t}}$ remains the same when comparing against T_2/3. In a follow-up study, we will test whether this behavior is universal when considering surface hydrogen or whether the chosen fraction of X_H = 0.2 coincidentally balances out other effects for a 20 M_⊙ He-burning star as we e.g. saw with the metallicity reduction being balanced by the surface hydrogen when considering only $\dot{M}$ in Fig. 5.

To get some insights on the general scaling of WR-type winds with effective temperature (here: T_2/3), we also compare our sequences to the trends for OB-type stars. In Fig. 14, we plot the trends obtained for the ULLYSES OB stars for the SMC and a Galactic comparison sample by Hawcroft et al. (in prep.) as well as the predictions from Vink & Sander (2021). We see that generally the terminal velocity increases much steeper with the temperature in OB-type winds than in WR-type winds. Only at higher temperatures (T_2/3 > 130 kK), when the winds become more optically thin as the mass-loss rates decrease, the steepness of the v_∞(T_2/3)-curves increases to a value more comparable to those obtained for OB-star winds.

4.3 Critical-point dependencies

In addition to studying the behavior of v_∞ as a function of the observable quantity T_2/3, we further investigate the behavior of v_∞ as a function of T_* or the physically more meaningful T_eff,crit. As noted above, the values of v_∞ tend to scatter a bit more, but if we restrict the linear fitting to the optically thick wind regime, we find $$\log \left({\upsilon }_{\infty }\left[\text{km\hspace{0.17em}}{\text{s}}^{-1}\right]\right)=1.2\log \left(T_{\text{eff}}\left({\tau }_{\text{crit}}\right)\left[\text{K}\right]\right)+\text{offset}.$$(9)

Given that the number of models in the optically thick regime is restricted and not all sequences reach far enough to see the flattening of the trend, this coefficient has to be considered as rather uncertain. Nonetheless, the corresponding scaling of ${\upsilon }_{\infty }\propto R_{\text{crit}}^{-0.6}$ can also be obtained from directly fitting v_∞(R_crit) in the optically thick limit. Using the critical radius to define the escape velocity $${\upsilon }_{\text{esc}}=\sqrt{\frac{2GM}{R_{\text{crit}}}}$$(10)

we obtain $${\upsilon }_{\infty }\propto {\upsilon }_{\text{esc}}^{1.2}.$$(11)

This relation also holds for the effective escape velocity $${\upsilon }_{\text{esc}}=\sqrt{\frac{2GM}{R_{\text{crit}}}\left(1-{\Gamma }_{\text{e}}\right)}$$(12)

since all of the model sequences have a constant L/M and Γ_e is approximately constant. The latter is consequence of the unchanged free electron budget below R_crit in the considered temperature range. Thus, in the optically thick wind regime we have a slight difference with ${\upsilon }_{\infty }\propto {\upsilon }_{\text{esc,eff}}^{1.2}$ along the T_* -dimension compared to the well-known v_∞ ∝ v_esc,eff in the well-known CAK theory (named after Castor et al. 1975). For the sequences along the L/M-dimension from Sander & Vink (2020) we instead obtain a negative trend of log v_∞ ≈ −4.6 log v_esc,ff + const. in the optically thick regime with a flattening of the trend for the highest mass-loss rates. In both cases, the scaling of v_∞ with v_esc,eff remains complicated with no straight-forward prediction as offsets remain in all scalings. This is in sharp contrast to the classical (m)CAK result, where v_∞ follows as an offset-free value from v_esc.

For lower mass-loss rates (log ${\dot{M}}_{\text{t}}<-4.5$), corresponding usually to T_eff,crit > 150 kK, we reach the regime where winds are mostly optically thin. Above, we could show that when reaching this regime, there seems to be an alignment of v_∞(T_2/3), with the slopes known from OB-type winds. With considerable scatter in the exponent of up to ±0.5, we find $${\upsilon }_{\infty }\propto T_{\text{eff},crit}^4,$$(13)

corresponding to ${\upsilon }_{\infty }\propto R_{\text{crit}}^{-2}$ or ${\upsilon }_{\infty }\propto {\upsilon }_{\text{esc,eff}}^4$, that is a steeper relation, contrary to the expected flattening of the slope. However, there is growing evidence from both observations (e.g., Garcia et al. 2014) as well as theoretical CMF-based and Monte Carlo calculations (e.g., Björklund et al. 2021; Vink & Sander 2021) that even in the typical OB-type regime the scaling of v_∞ with v_esc,eff is likely more complicated.

4.4 Influence on M_t

The scaling of v_∞ with R_crit introduces an additional dependency when considering the transformed mass-loss rate ${\dot{M}}_{\text{t}}$ as a function of R_crit or g_crit.

From Eqs. (5) and (6), we know that $\dot{M}\propto T_{\text{eff,crit}}^{-6}$ or $\dot{M}\propto R_{\text{crit}}^3$. From the definition of ${\dot{M}}_{\text{t}}$ (Eq. (2)) we get $$\log {\dot{M}}_{\text{t}}\left(R_{\text{crit}}\right)=\log \dot{M}\left(R_{\text{crit}}\right)-\log {\upsilon }_{\infty }\left(R_{\text{crit}}\right)+\text{offset}.$$(14)

With the different trends derived for the optically thick and thin limit in Sect. 4.3, we obtain ${\dot{M}}_{\text{t}}\propto R_{\text{crit}}^{3.6}$ and ${\dot{M}}_{\text{t}}\propto R_{\text{crit}}^5$ respectively. Using g_crit as defined in Eqs. (7) and (8), this yields $$\log {\dot{M}}_{\text{t}}=1.8\log g_{\text{crit}}+\text{offset}$$(15)

for the optically thick limit and $$\log {\dot{M}}_{\text{t}}=2.5\log g_{\text{crit}}+\text{offset}$$(16)

in the optically thin limit. These trends are depicted as sets of gray lines in Fig. 15, where the curves from the model sequences are shown as well. In contrast to the $\dot{M}\left(g_{\text{crit}}\right)$-behavior discussed in Sect. 3.4, the representation of the slope in the optically thin regime is less precise. In the optically thick limit, some curves align well, but others appear to be slightly steeper or shallower. Hence, the overall results for ${\dot{M}}_{\text{t}}\left(g_{\text{crit}}\right)$ should be considered less robust than the $\dot{M}\left(g_{\text{crit}}\right)$-trends. Interestingly, we do not see the “kink” or clear bending for some sequences at lowest (transformed) mass-loss rates that we see in Fig. 13 for $\dot{M}\left(g_{\text{crit}}\right)$. While it is hard to draw strong conclusions, there at least appears to be one continuous slope for ${\dot{M}}_{\text{t}}\left(g_{\text{crit}}\right)$ in the thinner wind regime, regardless of whether the critical point is located at temperatures below or above T_eff,crit. Since we find a change for $\dot{M}$ alone, this would imply that v_∞ outweighs this effect. While the inspection of the corresponding sequences in Fig. 14 is only indicative here, indeed the v_∞-curves of these sequences bend again toward shallower slopes then plotting them as functions of T_eff,crit.

Fig. 15 Transformed mass-loss rate ${\dot{M}}_{\text{t}}$ as a function of the gravitational acceleration at the critical radius $g_{\text{crit}}=GM{R}_{\text{crit}}^{-2}$. To reflect the expected trends from the Rcrit-fits, thin gray lines are plotted in the background. The dotted, gray lines indicate $\dot{M}\propto g_{\text{crit}}^{-2.5}$ (optically thin regime), while the dashed, gray lines correspond to $\dot{M}\propto g_{\text{crit}}^{-1.8}$ (optically thick regime).

5 Potential consequences for stellar evolution

Our study presents the very first sequences of hydrodynamically consistent atmosphere models in the cWR regime, where we vary the input parameter T_* − corresponding roughly to T_eff,crit for most models. WR mass-loss recipes commonly do not incorporate any temperature/radius dependency, which can be seen as a consequence of the optically dense winds of WR stars. When reproducing their spectra with prescribed velocity fields, there is a degeneracy of solutions making it impossible to find a unique value of T_* for more dense winds (e.g., Hillier 1991; Hamann & Gräfener 2004; Lefever et al. 2022).

From an evolutionary standpoint, one could justify the omission of a T_*-dependence arguing that hydrogen-free WN stars – and to some extend also WC stars – may form a 1D sequence as they represent He-burning stars that do not contain any further shell structure which could skew the relation between the luminosity and mass. In reality, effects such as inflation, convection, or rotation augment the physical conditions of wind launching and mass loss, especially when considering their multidimensional nature. However, in the currently typical 1D spherical approach ignoring such issues, no further parameter would be necessary if the He star evolution could be perfectly mapped to one of the fundamental stellar parameters. For He stars above 10 M_⊙, the intrinsic curvature of the HeZAMS indeed gets relatively small (e.g., Langer 1989) and the obtained tracks of WR evolution in different codes yield very similar temperatures around log(T [K]) = 5.1², regardless whether these have been calculated from pure He stars or including all prior evolution from the ZAMS (e.g., Georgy et al. 2012; Limongi & Chieffi 2018; Higgins et al. 2021).

5.1 Mass-loss comparison for a representative model

In the models from Sander & Vink (2020), we thus ignored the width of ≈0.1 dex in log T_* and fixed T_* in order to keep the total amount of models manageable. In this work, we now calculated a number of model sequences where we vary T_* in order to investigate the effect of a wider range of T_*, which probes not only the curvature of the He ZAMS, but also gives a glimpse of how $\dot{M}$ might be affected for stars which are not yet or no longer (exactly) on the He ZAMS. We find that despite the narrow range in temperature, the effect on $\dot{M}$ is quite noticeable. For our 20 M_⊙ model sequence (at Z_⊙) even a narrow range of only 0.05 dex (i.e., T_* = 125 … 140) results in a factor of two in $\dot{M}$. Whether such a significant correction is really necessary depends on the difference between the most realistic choice of T_* and the fixed value (141 kK) in Sander & Vink (2020). Combining the structural constraints by Grassitelli et al. (2018) with our model sequence, we find an ideal value of T_eff,crit ≈ T_* ≈ 130 kK for a 20 M_⊙ at Z_⊙ model without any hydrogen.

In Table 3, we provide a comparison of the resulting mass-loss rates for a 20 M_⊙ star. Beside the values employing the new estimate of T_* and the solutions for T_* = 141 kK from Sander & Vink (2020), we also list the resulting $\dot{M}$-values from Gräfener et al. (2017) and commonly used (semi-)empirical recipes. For Z_⊙ we find a difference of ~0.2 dex in $\dot{M}$, unaffected by the choice of D_∞. Using the values of Table 3 as an average mass-loss rate during the typical He burning lifetime (300 kyr), this corresponds to a difference between 12.5 M_⊙ and 8.1 M_⊙ at the end of core He-burning. This calculation is of course only a rough estimate and does not take any change of the stellar parameters or surface abundances into account. Nonetheless, the value using the $\dot{M}$ from Sander & Vink (2020) is close to what we obtain with actual stellar evolution calculations in Higgins et al. (2021).

At 0.5 Z_⊙, a value roughly corresponding to the LMC, the 0.2 dex shift holds as well when adopting D_∞ = 50. For the hydrogen-free 130 kK model at 0.5 Z_⊙ with D_∞ = 10, however, we only find a solution if we relax the stability criterion on $\dot{M}$ between consecutive updates that we otherwise enforce. For the 141 kK model, we already see a notable difference in M when reducing from D_∞ = 50 down to 10. The reason is that we are already close to the regime where we can no longer obtain a wind solution driven by the hot iron bump (see Sect. 3.2). For the 130 kK we have reached already a meta-stable situation with respect to the solution stability. Thus, the obtained value of $\dot{M}$ for D_∞ = 10 is much lower than expected. The matching of the absolute values for 130 kK and 141 kK is a pure coincidence with higher, also meta-stable solutions up to ≈ −4.7 for 130 kK being possible as well.

5.2 Structural limits and the role of hydrogen

The presence of hydrogen at the surface can considerably change the limits of the wind onset derived above. In contrast to our hydrogen-free results shown in Table 3 and depicted in Fig. 6, our model sequence with X_H = 0.2 and 0.5 Z_⊙ extends to much cooler temperatures (T_eff,crit < 100 kK) as the additional acceleration from free electrons helps to compensate the effect of the deceleration regions. The choice of D_∞ has some impact on the results, but in both cases the effect of surface hydrogen as such is much larger.

While we do not aim at a detailed comparison with observations in this work – which would require new analyses with dynamically-consistent models – it is striking that all hydrogen-free WN stars in the LMC are of the subtype WN4 or earlier (Hainich et al. 2014; Shenar et al. 2019). Moreover, all WC stars in the LMC show early subtypes as well. While one has to be careful drawing absolute conclusions, our temperature study now indicates that beside the metallicity limiting the lower luminosity of the observed WN population, the observed restriction of the subtype regime might be a direct consequence of the inability to launch WR-type winds below a certain (sonic point) temperature. From the perspective of fixed stellar parameters, the lower mass-loss rate reached at a lower metallicity corresponds to a shift to earlier subtypes (cf. Fig. 4), thereby confirming the suggestion by Crowther et al. (2002).

Coming from a different angle, but addressing essentially the same problem, Grassitelli et al. (2018) and Ro (2019) used hydrodynamic stellar structure models and semi-analytic approaches to predict the existence of a “minimum mass-loss rate” for launching a stellar wind from the hot iron bump. For values of $\dot{M}$ below this limit, extended low-density regions were predicted. In this work, we do not aim to obtain the latter type of solutions, but the calculations failing to launch a wind show a tendency toward trying to launch a wind further out with a (much) lower $\dot{M}$. We can thus qualitatively confirm the structural predictions by Grassitelli et al. (2018) and Ro (2019), assuming that we are bound to a choice of T_* following the – ideally hydrodynamical – structure calculations for the HeZAMS. This underlines once more that unifying structural and atmosphere models remains a challenge that requires a new generation of both atmosphere and structure models.

5.3 An approximated handling of the temperature shift

In light of the structural considerations above, it appears likely that any future description of WR-type mass loss needs a temperature or radius-dependency. Simpler treatments might be reasonable in a time-averaged situation, but cannot predict a realistic mass loss for individual points along an evolutionary track. Hence, a detailed update of the Sander & Vink (2020) formula will eventually be necessary, but the current amount of models does not allow a wide-space parameter investigation. The applicability to lower masses also turns out to be nontrivial: One the one hand, the curvature of the HeZAMS toward cooler temperatures should soften the sharp drop obtained in Sander & Vink (2020) of $\dot{M}$ toward lower He star masses. On the other hand, we reach lower limits of T_eff,crit for driving winds by the hot iron bump (cf. Sect. 3.2). In our 12.9 M_⊙ sequence with X_H = 0.2, this limit is at ≈97 kK. Given that this limit seems to increase to slightly higher temperatures for lower L/M values, it appears unlikely that one can find any solution for a wind driven by the hot iron bump for stars with M ≤ 10 M_⊙ fulfilling the L−M relation from Gräfener et al. (2011).

While a full coverage of the driving limit at cooler temperatures will require its own tailored study, we can use Eq. (5) from Sect. 3.4 to derive a decent temperature description up to discontinuity in $\dot{M}$. Since Eq. (5) seems to be valid across both the optically thick and thin regime, we can approximate $\dot{M}$ for WN winds driven by the hot iron bump via $$\log \left(\frac{\dot{M}}{M_{\odot }{\text{yr}}^{-1}}\right)=\log \left(\frac{{\dot{M}}_{\text{SV2020}}}{M_{\odot }{\text{yr}}^{-1}}\right)+3\log \left(\frac{R_{\text{crit}}}{R_{\text{crit,}T141}}\right)$$(17)

with ${\dot{M}}_{\text{SV2020}}$ denoting the mass-loss rate from Sander & Vink (2020) and R_crit,T141 = R_crit (T_* = 141 kK) being the critical radius (in R_⊙) of their corresponding model (e.g., 1.217 R_⊙ for the 20 M_⊙ He star without hydrogen). Although R_crit,T141 could be obtained from L/M or Γ_e via a nonlinear fit of the Sander & Vink (2020) data, it is much more convenient to reformulate Eq. (17) in terms of the effective temperature at the critical (≈ sonic) point T_eff,crit, yielding $$\log \left(\frac{\dot{M}}{M_{\odot }{\text{yr}}^{-1}}\right)=\log \left(\frac{{\dot{M}}_{\text{SV2020}}}{M_{\odot }{\text{yr}}^{-1}}\right)-6\log \left(\frac{T_{\text{cff,crit}}}{141\text{\hspace{0.17em}kK}}\right).$$(18)

Apart from small deviations for the highest mass-loss rates $\left(\dot{M}\gg {10}^{-4}{M}_{\odot }{\text{yr}}^{-1}\right)$, the fixed value of 141 kK accurately represents the value of T_eff,crit in Sander & Vink(2020), as illustrated previously in Fig. 8.

This adjusted $\dot{M}$-recipe requires the knowledge of either T_eff,crit or R_crit. As we did include only a small microturbulent velocity in our modeling efforts (30 km s⁻¹), the quantities can be replaced by the sonic point values without introducing a considerable error. Still, the accurate use of Eqs. (17) and (18) requires models with a meaningful sonic point in a hydrodynamical sense to prevent reintroducing any further radius/temperature discrepancies. Stellar atmosphere analyses typically employ predefined velocity fields (usually β-laws) and thus do not have a sonic point that is hydrodynamically consistent. Purely hydrostatic stellar structure calculations are problematic as well as they do not yield a sonic point by construction. This underlines that in order to obtain a really insight- and meaningful comparison between theory and observation for optically thick winds, a new generation of both atmosphere and stellar structure models will be necessary.

The results obtained in our study could be helpful to eventually obtain realistic predictions for the effective temperatures (T_2/3) of WR stars in stellar evolution models. A route toward such a recipe based on our findings is given in Appendix D.

6 Transparency to He II ionizing photons

Despite having intrinsically quite hot temperatures, classical WR stars do not necessarily emit a significant number of ionizing photons beyond the He II ionization edge, that is below 227 Å or above 54 eV. As first described in Schmutz et al. (1992), the transparency of the wind for photons with energies above 54 eV depends on the mass-loss rate, and thus the density of the wind. In more dense winds, He III recombines to He II, making the atmosphere opaque to He II ionizing photons out to very large radii. The presence of line blanketing further affects the absolute ionizing fluxes significantly (cf. Smith et al. 2002). Beside usually leading to a reduction of the He I ionizing flux, it can also affect the He II ionizing flux transition by a few orders of magnitude as we will see in our model calculations. To cover the region where the (continuum) optical depth drops below unity in this wavelength region, we extended the outer boundary radius R_max of our atmosphere models to extremely large values, often up to 100000 R_*. (Typical atmosphere models for spectral fitting require only R_max = 100 … 1000 R_*.)

In Fig. 16, we plot the rate of ionizing photons per second Q_{He II} as a function of the transformed mass-loss rate ${\dot{M}}_{\text{t}}$. The absolute numbers in the regime with low Q_{He II} are more uncertain as they depend strongly on the precise boundary treatment. If the wind becomes transparent around and below 227 Å only very close to the model boundary or even remains optically thick at R_max, the value for Q_{He II} can be underestimated, but should never exceed 10⁴¹ s⁻¹. Given that this is many orders of magnitude below the actually strong Q_{He II} emitters with rates >10⁴⁷ s⁻¹, the values of the shaded regime in Fig. 16 should have no practical consequences.

Displaying the He II ionizing flux as a function of ${\dot{M}}_{\text{t}}$ in Fig. 16 confirms that similar to other quantities, the switch in transparency is caused by the lower wind density. In fact, there seems to be a critical lower boundary around $\log \left({\dot{M}}_{\text{t}}\left[M_{\odot }{\text{yr}}^{-1}\right]\right)=-4.6$ to −4.4 where all model sequences switch abruptly. To study whether this transition value might be more universal, we created the same plot for the model sequences from Sander & Vink (2020). Their model set, shown in Fig. 17, clearly hints at a Z-dependency for the transition, which we do only sparsely map in our new model set. In fact, our new sequences are quite complementary to the datasets from Sander & Vink (2020), indicating that the transition does not only depend on Z as a total value, but likely on the detailed composition and – notably especially at the lower end of the transition in Fig. 16 – on the choice of the clumping factor. Nonetheless, we can conclude that all stars in our large model sample with $\log \left({\dot{M}}_{\text{t}}\left[M_{\odot }{\text{yr}}^{-1}\right]\right)<-4.6$ are strong emitters of He II ionizing flux. Thus, we propose to take this value as an upper limit of whether to consider WR stars as notable contributors to the He II ionizing photon budget.

Fig. 16 Number of helium ionizing photons per second QHe II as a function of the transformed mass-loss rate ${\dot{M}}_{\text{t}}$ for the new model sequences calculated in this work.

Fig. 17 Number of helium ionizing photons per second QHe II as a function of the transformed mass-loss rate ${\dot{M}}_{\text{t}}$ for the model sequences computed in Sander & Vink (2020).

7 Summary and conclusions

In this work, we presented an exploratory study for the temperature-dependency of radiation-driven winds launched by the so-called hot iron opacity bump. For the first time, we calculated temperature-dependent sequences of hydrodynamically consistent stellar atmosphere models in the cWR regime. To achieve our results, we had to allow for nonmonotonic velocity field solutions when solving the hydrodynamic equation of motion. In order to perform the necessary radiative transfer in the comoving frame, we afterwards interpolated the obtained velocity fields such that the main wind properties ($\dot{M}$, v_∞) as well as the characteristics in the outer wind were maintained. We draw the following conclusions:

• The mass-loss rates $\dot{M}$ depend significantly on the critical radius R_crit and thus also on the assumed model temperature setting T_eff(R_crit). For model sequences with constant luminosity L and stellar mass M, we obtain $\dot{M}\propto R_{\text{crit}}^3$ over a wide range with moderate deviations from this purely geometrical effect occurring at the lower and upper end of our sequences. This finding can also be expressed in the form of $\dot{M}\propto g_{\text{crit}}^{-\text{3/2}}$, reflecting that larger radii for the critical point imply a lower gravitational force. Our findings underline that WR-type mass-loss depends on multiple parameters and the 2D description from Sander & Vink (2020) needs to be extended further to describe all relevant effects.

• Except for very dense winds – corresponding to $\log \left({\dot{M}}_{\text{t}}\left[M_{\odot }{\text{yr}}^{-1}\right]\right)\approx -3.0$ and above – the effective temperature at the critical point T_eff(τ_crit) is close to the effective temperature at a Rosseland continuum optical depth of τ_R,c = 20. For WN-type models τ_R,c = 20 typically corresponds to τ_Thom ≈ 17, albeit with considerable scatter along the model sequences.

• We find a characteristic value of $\left({\dot{M}}_{\text{t}}\left[M_{\odot }{\text{yr}}^{-1}\right]\right)\approx -4.5$ for the transition between the optically thin and thick regime. While there is some scatter between different model sequences, this characteristic value of ${\dot{M}}_{\text{t}}$ (plus some error margin) provides a very convenient tool to distinguish between the regimes as ${\dot{M}}_t$ can also be determined with empirical models. Known WC stars show values well above this (e.g., Gräfener & Vink 2013; Sander et al. 2019) while WO stars might be found on both sides of the transition. Whether the characteristic value of ${\dot{M}}_t$ also holds for winds that might not be driven by the hot iron bump is currently unclear. We calculated the transformed mass-loss rates for stars at the spectral transition from Of to WNh, which likely happens at a cooler temperature regime than studied in this work³. Their corresponding transformed mass-loss rates ${\dot{M}}_{\text{t,trans}}$ appear to be below −4.5, e.g. at $\log \left({\dot{M}}_{\text{t,trans}}\left[M_{\odot }{\text{yr}}^{-1}\right]\right)\approx -5.0$ in the Arches cluster (Martins et al. 2008; Vink & Gräfener 2012) and even lower in R136 (Bestenlehner 2020; Bestenlehner et al. 2020).

• The choice of the maximum clumping factor D_∞ does not affect our derived $\dot{M}\left(T_{2/3}\right)$ trends, but leads to an additional shift in the obtained relations with higher clumping factors corresponding to higher $\dot{M}$ for the same T_2/3. In contrast to all shifts introduced by varying fundamental stellar parameters or abundances, the shift due to D_∞ does not vanish when considering ${\dot{M}}_{\text{t}}\left(T_{2/3}\right)$ instead of $\dot{M}\left(T_{2/3}\right)$.

• In the limit of optically thick winds, we obtain a linear relation between log T_2/3 and log ${\dot{M}}_{\text{t}}$, independent of chemical composition (but for a fixed clumping factor). Combined with the also Z-independent result from Sander & Vink (2020) that log ${\dot{M}}_{\text{t}}\propto \log \left(L/M\right)$, this could provide an easy-to-use prediction for WR effective temperatures in stellar structure and evolution models.

• Classical WR stars and non-WR helium stars with $\log \left({\dot{M}}_{\text{t}}\left[M_{\odot }{\text{yr}}^{-1}\right]\right)<-4.6$ are strong emitters of He II ionizing flux (with Q_{He II} > 10⁴⁸ s⁻¹). Helium stars with stronger winds are (mostly) opaque to radiation above 54 eV and thus should not be considered as sources of hard ionizing radiation.

• Albeit being limited in comparability to particular observed targets, our findings indicate that high clumping factors (D ≈ 50) might be necessary to reproduce the observed combinations of $\dot{M}$ and v_∞. This is in sharp contrast to the first results obtained from 3D wind modeling by Moens et al. (2022) arguing for much lower clumping factors of D ≈ 2. At present, the reason for this discrepancy is unclear. Various solutions are possible, e.g., missing opacities in our wind models – where the presently assumed high clumping would act as a “fudge factor” to make up for that. Alternatively, the sharp contrast in the clumping factor might simply be the result of a mismatch between the considered regimes. Currently, the 3D models from Moens et al. (2022) probe only the wind onset where even in our 1D models we assume D(r) ≪ D_∞ with D(τ_crit) typically ranging between 1.5 and 4.

• When comparing empirically obtained results in the T_2/3-${\dot{M}}_{\text{t}}$-plane to our derived curves, we find a significant fraction of stars to have lower values of ${\dot{M}}_{\text{t}}$ than predicted by our curves using hydrodynamic models. It is currently unclear whether this is due to a deviation from the theoretical setup in this work (e.g., different clumping stratification, other L-M combinations) or inherent simplifications in the empirical analyses (e.g., the use of a β-type velocity law). A dedicated analysis of individual objects with hydrodynam-ical model atmospheres will be necessary to uncover the origin of this discrepancy.

• The limits of driving optically thick winds crucially depend on our knowledge of opacities. In case of considerable changes – e.g., a higher iron opacity as reported by Bailey et al. (2015) for the so-called “deep iron bump” at T_e ≈ 2 × 10⁶ K – wind quantity predictions such as $\dot{M}$ and v_∞ could shift significantly. Moreover, our understanding of the limits of radiative driving would be affected as well, e.g., due to strengthening or weakening the bumpy radius dependency of the flux-weighted mean opacity ℌ_F. Beside the impact of multi-D effects, higher (Fe) opacities could play an important role to resolve current discrepancies, such as the lower luminosity end of the LMC WN population or the aforementioned need for higher D_∞ to reach the observed terminal velocities.

With these conclusions, our study underlines the complexity of radiation-driven mass loss, revealing both parameter regimes with a clear scalings and characteristic transitions as well as more obscure parameter regions where $\dot{M}$ appears to break down suddenly. We provide an adjustment of the recent $\dot{M}$-description from Sander & Vink (2020) to account for different radii (or effective temperatures respectively) and emphasize that the model efforts presented there as well as in this work were limited to the regime where winds are launched by the hot iron bump. We thus consider our work as an intermediate step on the way toward a more comprehensive understanding of WR-type mass loss and will expand to other regimes in future studies.

Acknowledgements

The authors would like to thank the anonymous referee for their careful and constructive comments and suggestions. AACS and VR acknowledge support by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) in the form of an Emmy Noether Research Group – Project-ID 445674056 (SA4064/1-1, PI Sander). RRL is funded by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) – Project-ID 138713538 – SFB 881 (“The Milky Way System”, subproject P04). LP acknowledges support by the Deutsche Forschungsgemeinschaft – Project-ID 496854903 (SA 4046/2-1, PI Sander). JSV is supported by STFC funding under grant number ST/V000233/1. This publication has benefited from discussions in a team meeting (PI: Oskinova) sponsored by the International Space Science Institute (ISSI) at Bern, Switzerland. A significant number of figures in this work were created with WRPLOT, developed by W.-R. Hamann.

Appendix A $\dot{M}$-temperature Fit

For all of our model sequences, the data points in the log $\dot{M}$-log T_eff(τ_crit)-plane suggest a linear relation between the two quantities with deviations occurring only close to the wind driving limit (corresponding to the maximum $\dot{M}$ in Fig. A.1). In the fits, we thus exclude the uppermost 0.15 dex in log $\dot{M}$. The resulting fit coefficients for the slope including their error margins are given in Table A.1. For each individual sequence, the mass loss can be well described for T_eff(τ_crit) > T_eff,crit,min and $\dot{M}<{\dot{M}}_{\max }$ by $$\log \left(\frac{\dot{M}}{M_{\odot }{\text{yr}}^{-1}}\right)=-6\log \left(\frac{T_{\text{eff,crit}}}{141\text{\hspace{0.17em}kK}}\right)+\log \left(\frac{{\dot{M}}_{\text{offset}}}{M_{\odot }{\text{yr}}^{-1}}\right)$$(A.1)

with the value for ${\dot{M}}_{\text{offset}}$ being different for each model sequence. Table A.1 also lists these coefficients together with the corresponding validity limits T_eff,crit,min and ${\dot{M}}_{\max }$.

Fig. A.1 Linear fits (solid lines) to the Mass-loss rate $\dot{M}$ as a function of Teff(τcrit) in a double-logarithmic-plane. The fit coefficients for the different datasets are given in Table A.1.

Appendix B T_2/з temperatures for the Sander & Vink (2020) sample

The effective temperatures T_2/з at τ_Ross = 2/3 resulting from the model sequences in Sander & Vink (2020) are depicted in Fig. B.1. Similar to what we obtain when varying T_*, cooler values of T_2/з require higher mass-loss rates. As T_* is fixed to ≈ 141 kK in Sander & Vink (2020), higher luminosities or l/m-ratios are required to reach higher mass-loss rates for the same Z. At lower metallicity, stars have to get closer to the Eddington Limit to reach sufficient mass loss, shifting the onset of the drop in T_2/з to higher l and steepening in particular this drop. The differences in T_2/з and the range of luminosities covered has quite some interesting implications on the spectral appearance of the stars and consequently also which WR subtypes one would expect in a certain environment. Assuming that at least all the winds of early-type WR stars are launched at the hot iron bump, the temperatures of the lowest luminosity WN stars should get hotter at low Z. This seems to be the case when comparing the WN populations in the Milky Way and the LMC (e.g., Hamann et al. 2019; Shenar et al. 2020), but further – ideally hydrogen-free – WN populations in other Galaxies need to be studied to draw any firm conclusions. In the SMC, apart from the small sample size, all WN stars contain hydrogen and might not align with the L-M relation we assume in Sander & Vink (2020) and this work.

Fig. B.1 HRD with the effective temperature T2/з defined at a Rosseland optical depth of τRoss = 2/3 and the model luminosity L for our sets of He ZAMS models at different Z. For comparison, the HeZAMS (gray, dashed) and ZAMS (gray, solid) are shown as well.

Appendix C The effect of enhanced clumping on the radiative acceleration

The effect of clumping on our hydrodynamic wind solutions is not trivial. Given that we only use the so-called microclumping approximation assuming optically thin clumps and the solution of radiative transfer in the comoving frame is performed with the average density and not the clumped density, one might expect that the choice of D_∞ could have no effect at all. However, this is not the case. Different values of D_∞ affect the population numbers which in turn affect the radiative transfer. In particular, higher choices of D_∞ favor recombination in the wind. For most elements, lower ionization stages provide more opacity and thus more radiative acceleration.

In Figs. C.1 and C.2 we present the resulting acceleration contributions for the hydrogen-free 20 m_⊙ WN models with D_∞ = 50 and 10, respectively. To obtain these curves, the individual opacities resulting from the different ions are stored in addition to the total opacity. Beside the total calculation of $$a_{\text{rad}}\left(r\right)=\frac{4\pi }{c}{\displaystyle \underset{0}{\overset{\infty }{\int }}{\varkappa }_{\nu }}\left(r\right)H_{\nu }\left(r\right)\text{d}\nu$$(C.1)

similar integrals to Eq. (C.1) are calculated using only the ion-specific opacities (e.g. $h_{\nu }^{\text{Fe\hspace{0.17em}V}}$) instead of the total ℌ_v, yielding the specific acceleration contribution for each ion. The corresponding wind parameters of the two displayed models and two other sets with varying D_∞ are listed in Table C.1. With out fixed characteristic velocity for the clumping onset of v_cl = 100 km s⁻¹, the depicted models increase from almost no clumping to D_∞ within the hot iron bump. Thus, the mass-loss rates are barely affected, but the terminal velocity increases significantly from 1186 km s⁻¹ to 1754 km s⁻¹.

The comparison of Fig. C.2 with Fig. C.1 confirms that the additional opacity to reach the higher terminal velocity is provided by lower ions, most notably Fe IV, which is the leading accelerator in the outer wind for the model with D_∞ = 50, while Fe V remains in the lead for the model with D_∞ = 10. In both cases Fe V is the most populated Fe ion in the outermost wind, while the “fresh” opacity provided by the lesser populated Fe IV is most efficient for the line acceleration in the case of D_∞ = 50. In the case of D_∞ = 10, the population of Fe IV instead is too low to contribute significantly. The change in v_∞ is further enlarged by the significantly smaller deceleration zone in the D_∞ = 50 model. In the deeper wind layers, the higher clumping boosts the contribution from the iron M-shell opacities and leads to an increased bound-free contribution (i.e. recombination) from He II.

The two other examples in Table C.1 illustrate that in some cases also the mass-loss rate can be notably affected by changes of D_∞. In our models, this is a consequence of the fixed value of v_cl. For regimes where generally lower values of v_∞ are reached, e.g. in lower metallicity model set, often the whole amount of acceleration is reduced, shifting also the region with v ≈ 100 km s⁻¹. This can then have two effects leading to a lower $\dot{M}$ for lower values of D_∞: first, a direct reduction of opacities in the region that determines $\dot{M}$. In addition, the reduced wind density could push the star out of the regime where the critical point is in a totally optically thick region (cf. Sander & Vink 2020), which would lead to a further reduction in $\dot{M}$.

In the last column of Table C.1, we provide the resulting transformed mass-loss rates ${\dot{M}}_{\text{t}}$. While there is already a scaling with $\dot{M}\sqrt{D_{\infty }}$ in these, it does not compensate the clumping changes as the square root of the D_∞-ratios is much larger than the changes in v_∞ (and $\dot{M}$). For example, the hydrogen free models differ by $\sqrt{50/10}\approx 2.24$, while v_∞ only increases by a factor of 1.48. Therefore, the models with higher D_∞ posses a higher ${\dot{M}}_{\text{t}}$, despite larger terminal velocities reducing its value.

Fig. C.1 Contributions of the different ions to the radiative acceleration of a hydrodynamically consistent, hydrogen-free WN model with T* = 130 kK, log L/L⊙ = 5.7, M = 20 M⊙, and D∞ = 50: Different ions are denoted by a combination of different color and symbol. The total radiative acceleration (arad), the Thomson acceleration from free electrons (aThom = Γe · g), and the contribution from gas (and turbulence) pressure (apress) are also shown for comparison. The loosely dashed horizontal line denotes the total Eddington limit that needs to be overcome to launch a wind.

Fig. C.2 Contributions of different ions to the radiative acceleration, plotted similar to Fig. C.1, but now for a model employing D∞ = 10. The wind reaches a lower terminal velocity and the supersonic region with Γrad = arad/g < 1 is more pronounced than in the model with D∞ = 50.

Appendix D Estimating effective temperatures in stellar structure models

For stellar structure models, the occurrence of optically thick winds usually spoils the straight-forward prediction of the observable effective temperature T_2/3 (see, e.g., Groh et al. 2014, for a more detailed discussion). In the previous Sect. 3.3, we obtained that for a given clumping factor D_∞, our model sequences collapse almost perfectly to a single line in the ${\dot{M}}_{\text{t}}$-T_2/3-plane, yielding $$T_{2/3}\propto {\dot{M}}_{\text{t}}^{-1/2}$$(D.1)

for log $\left({\dot{M}}_{\text{t}}\left[M_{\odot }{\text{yr}}^{-1}\right]\right)>-4.5$, thereby providing us with a potential path to predict the observable effective temperature. The relation (D.1) seems to be approximately unaffected by abundance (X_i) changes, but there is a clear offset for different choices of D_∞. In our calculations, there is a difference of ∆ log (T_2/3 [K]) ≈ 0.08 … 0.1 between D_∞ = 4 and D_∞ = 10 and ∆ log (T_2/3 [K]) ≈ 0.1 … 0.15 between D_∞ = 10 and D_∞ = 50, but the current amount of data along the D_∞ plane is insufficient to provide a robust mathematical formula that could enable a scaling with D_∞.

A different slope of ≈ −2/3 was obtained in Sect. 3.3, when considering the sample of Sander & Vink (2020) instead of our new model sequences. The origin of the difference in the slopes must be rooted in the different nature of the sequences: In the new sequences calculated for this work, L and M are constant. For higher mass-loss rates $\dot{M}$ we then obtain lower values of v_∞ along a sequence. In the sequences from Sander & Vink (2020), where we proceed to higher L/M-ratios along each dataset, such a trend between $\dot{M}$ and v_∞ is only reached in the optically thin part, while we obtained $\dot{M}\propto {\upsilon }_{\infty }$ in the dense wind regime. As a consequence, models from the two different sources with approximately the same value of $\dot{M}$ will differ in their v_∞. When comparing the v_∞-values, the models from the T_*-sequences in this work will have lower terminal velocities and thus their ${\dot{M}}_t$ will be higher. Arguing that $\dot{M}$ is the major factor in setting the T_2/3 value, we can thus conclude that this difference in the v_∞ trends leads to the steeper slope for the sequences along the L/M-domain.

Given the focus of this work on the temperature trends and the fact that the steep linear trends in Fig. 12 do not provide a good description around the transition region of $\left({\dot{M}}_{\text{t}}\left[M_{\odot }{\text{yr}}^{-1}\right]\right)\approx -4.5$, we therefore suggest the more shallow formula $$\log \left(T_{2/3}\left[\text{K}\right]\right)=2.9-0.5\log \left({\dot{M}}_{\text{t}}\left[M_{\odot }{\text{yr}}^{-1}\right]\right)$$(D.2)

as a first attempt to approximate the observable effective temperature T_2/3 of a WR star in stellar structure models, which is essentially a rounded version of Eq. (3). This formula implicitly assumes D_∞ = 50 and is recommended for ${\dot{M}}_t>{10}^{-4.5}{M}_{\odot }{\text{yr}}^{-1}$. For lower estimates of D_∞, the offset value of 2.9 would need to be reduced by about 0.1 … 0.2. We emphasize that Eq. (D.2) is a first approach that needs to be tested and likely refined in future studies.

With Eq. (D.2) given, only the transformed mass-loss rate M_t needs to be known to determine T_2/3. In Sander & Vink (2020), we could show that for T_* = 140 kK the quantity ${\dot{M}}_{\text{t}}$ is practically independent of metallicity in the limit of optically thick winds (“pure WR regime”). There, ${\dot{M}}_{\text{t}}$ can be expressed as a linear function of L/M with a possible deviation only occurring for He stars with current masses above 50 M_⊙. To check whether this conclusion is independent of T_*, we calculate a small set of models sequences with different L/M values for different T_eff,crit ≈ T_*. The resulting trends are shown in Fig. D.1. It is clear from Fig. D.1 that there is some uncertainty in the slopes as well as a potential dependence of the slopes on T_* itself, but in general an approximately linear behavior is obtained for each choice of T_eff,crit ≈ T_*. Hence, we obtain a viable prediction method for stellar evolution models. This method is particularly elegant as it does not require any further assumptions about the flux-weighted mean opacity or the shape of the velocity field as for example necessary in the current wind-corrected temperatures in the GENEC models (see, e.g., Groh et al. 2014).

To get a formula for ${\dot{M}}_{\text{t}}$ that only depends on quantities which can be obtained from stellar structure calculations, we can use the result derived in Sect. 4.4. Considering that in the new model sequences calculated for this work both L and D∞ are constant within one sequence, we can conclude that ${\dot{M}}_{\text{t}}\propto \dot{M}\text{/}{\upsilon }_{\infty }$ and obtain $$\log \left({\dot{M}}_{\text{t}}\left[M_{\odot }{\text{yr}}^{-1}\right]\right)=-7.2\cdot \log \left(T_{\text{eff}}\left({\tau }_{\text{crit}}\right)\left[\text{K}\right]\right)+\text{offset}$$(D.3)

or ${\dot{M}}_{\text{t}}\propto R_{\text{crit}}^{3.6}$ in the regime of optically thick winds, i.e. for log $\left({\dot{M}}_{\text{t}}\left[M_{\odot }{\text{yr}}^{-1}\right]\right)>-4.5)$.

In a second step, we then merge Eq. (D.3), which has been determined for sequences of constant L/M, with the L/M-dependence obtained in Sander & Vink (2020). Together, we synthesize the formula $$\log \frac{\dot{M}}{M_{\odot }{\text{yr}}^{-1}}=1.25\log \frac{L/M}{L_{\odot }/M_{\odot }}+3.6\log \frac{R_{\text{crit}}}{R_{\odot }}+{\dot{M}}_{\text{t,off}}\left(X_i,D_{\infty }\right).$$(D.4)

Again, it might be more convenient to replace R_crit with T_eff,crit and gauge this with the 20 M_⊙ model at 141 kK. This then yields $$\log \frac{\dot{M}}{M_{\odot }{\text{yr}}^{-1}}=1.25\log \frac{L/M}{L_{\odot }/M_{\odot }}-7.2\log \frac{T_{\text{eff,crit}}}{141\text{kK}}-\mathrm{9.39.}$$(D.5)

Fig. D.1 Transformed mass-loss rate ${\dot{M}}_{\text{t}}$ as a function of L/M for different sequences varying in T*. All models use XH = 0.2 except the dashed-dotted sequences having XH = 0 for comparison. Apart from the red, dashed sequence (using D∞ = 10), all sequences employ D∞ = 50. The gray, dotted line is an interpolation of the solid sequences along the theoretical temperatures for the He ZAMS from Grassitelli et al. (2018) and Langer (1989). The light dotted linear curves in the background indicate the slope of 1.25 which is used in Eq. (D.4) and onward.

In Eq. (D.5), we also dropped the offet ${\dot{M}}_{\text{t,off}}\left(X_i,D_{\infty }\right)$ which contains further, uncertain dependencies. These can alter the resulting values of ${\dot{M}}_{\text{t}}$, e.g. by ≈ −0.2 dex when changing from D_∞ = 50 to 10. Inserting Eq. (D.5) into Eq. (D.2), we obtain the final formula for estimating T_2/3: $$\log \left(\frac{T_{2/3}}{\text{K}}\right)=7.595-0.625\log \frac{L/M}{L_{\odot }/M_{\odot }}+3.6\log \frac{T_{\text{eff,crit}}}{141\text{kK}},$$(D.6)

This formula is only valid for hydrogen-free WN stars as we have considerable offsets for other chemical compositions in ${\dot{M}}_{\text{t}}\left(L/M\right)$ (cf. Fig. D.1). While the conversion between ${\dot{M}}_{\text{t}}$ and T_2/з is unaffected by chemical composition (cf. Sect. 3.3), the resulting radiative acceleration is not. For example, the additional acceleration from free electrons in partially stripped stars with remaining surface hydrogen leads to higher mass-loss rates than in H-free stars of the same L/M-ratio (cf. Fig. 5), thereby substantially shifting the balance between $\dot{M}$ and υ_∞ and the resulting ${\dot{M}}_{\text{t}}$-values. In a future study, we thus plan to extend Eq. (D.6) by a hydrogen-dependent term. While various uncertainties, e.g., of the precise slopes in ${\dot{M}}_{\text{t}}\left(L/M\right)$ and $T_{2/3}\left({\dot{M}}_{\text{t}}\right)$ limit the accuracy of Eq. (D.6), it is sufficient enough to tell whether observed effective temperatures of predicted objects are in the range of e.g. 100, 50, or only 20 kK. Depending on the scientific context, such differences in T_2/з can have a big impact. While we do not have enough data to draw larger conclusions, the thick gray dotted line Fig. D.1 illustrates that on the He ZAMS, compact radii of the stars at the critical point might even outweigh an expected increase in ${\dot{M}}_{\text{t}}$ due to a larger L/M. Nonetheless, we can present the first estimate of T_2/з for WN stars derived from fundamental principles which can be readily applied in stellar evolution models and population synthesis. The validity of the assumptions made here will have to be tested within dedicated test calculations in stellar evolution models and benchmarked with WR observations analyzed with traditional as well as dynamically consistent models.

Appendix E Additional Figures

In addition to Fig. 5 discussed in Sect. 3, we plot the mass-loss rate $\dot{M}$ as a function of T_* in Fig. E.1. For very high mass-loss rates, we see a bending of the curves toward lower values of T_*. This is a consequence of the deeper wind launching, which in this regime happens further in than the defining optical depth for T_* (i.e. at τ_R,cont > 20). Numerically, these models use a higher optical depth as their inner boundary and we then determine T_*(τ_R,cont = 20) for a better comparison with the rest of the model calculations. In this regime, T_* is no longer a good approximation for T_eff,crit.

Fig. E.1 Analogous plot to Fig. 5, but now showing the mass-loss rate $\dot{M}$ as a function of T* for our model sequences.

Fig. E.2 Analogous plot to Fig. E.1, but now showing the transformed mass-loss rate ${\dot{M}}_{\text{t}}$ instead of the normal $\dot{M}$.

Fig. E.3 Mass-loss rates as a function of different temperature scales for a series of dynamically consistent atmosphere models with log L/L⊙ = 5.35, M = 12.9 M⊙, and XH = 0.2: The thick red dashed line denoted the classical effective temperatures defined at a Rosseland optical depth of τR = 2/3, while the green solid line and the blue dashed-dotted lines denotes the effective temperatures referring to τcrit and τR,cont = 20 respectively. The green dotted line on the right denotes the (electron) temperature at the critical point. Curves in lighter colors reflect models using the simple integration treatment suppressing negative velocity gradients.

Fig. E.4 Analogous plot to Fig. E.3, but now for the WC model series with log L/L⊙ = 5.7, M = 20 M⊙, and 0.5 Z⊙.

To eliminate the effect of different clumping factors and any remaining distance uncertainties, it is helpful to consider the transformed mass-loss rate ${\dot{M}}_{\text{t}}$ instead of $\dot{M}$. In Fig. E.2, we show the analogous plot to Fig. E.1 with $\dot{M}$ being replaced by ${\dot{M}}_{\text{t}}$. While the vertical spread in the observations is slightly reduced, the general temperature mismatch between the empirical T_* and our model sequences remains, highlighting once more the “Wolf–Rayet radius problem” seen in traditional atmosphere analyses.

In an extend to Fig. 6 and Fig. 7, we show similar plots for the model sequences with 12.9 M_⊙, X_H = 0.2 and Z = Z_⊙ in Fig. E.3 and the WC model sequence with 20 M_⊙ and Z = 0.5 Z_⊙ in Fig. E.4. Similar to the result obtained for the 15 M_⊙-sequence discussed in Sect. 3.2, there is a lower minimum temperature for obtaining wind solutions driven by the hot iron bump.

References

All Tables

All Figures

Fig. 1 Illustrating example of the possible velocity treatments in case that the integration of the hydrodynamic equation of motions yields a nonmonotonic υ(r) (solid blue curve): in the simple treatment, negative velocity gradients are suppressed during the integration, yielding the blue, dash-dotted curve. In some cases, such as illustrated in this plot, this can spoil the terminal velocity υ∞. In the more sophisticated treatment, negative gradients are therefore taken into account and υ(r) is modified such that the interpolated solution keeps the obtained υ∞.

In the text

Fig. 2 Resulting acceleration stratification of the converged models with two different treatments of nonmonotonic velocity fields: The dashed red and black lines illustrate the two sides of the hydrodynamic equation of motion for a model where negative velocity gradients are ignored in the solution. The corresponding solid lines show the result for the alternative method where the nonmonotonic υ(r) is interpolated afterwards. The small inlet shows the resulting velocity fields for both models. In this example we show models with T* = 125 kK, log L/L⊙ = 5.475 and M = 15 M⊙.

In the text

Fig. 3 Major contributions to the radiative acceleration for two hydrodynamically consistent, hydrogen-free WN models with T* = 130 kK, log L/L⊙ = 5.7, M = 20 M⊙, and D∞ = 10 (upper panel) or D∞ = 50 (lower panel): for the line contributions, all elemental contributions except Fe are summed over all ions. The total radiative acceleration (arad), the Thomson acceleration from free electrons (aThom = Γe · g), and the contribution from gas (and turbulence) pressure (apress) are also shown for comparison. The loosely dashed horizontal line denotes the total Eddington limit that needs to be overcome to launch a wind.

In the text

	Fig. 4 Example spectra from our WN model sequence with log L/L⊙ = 5.7 and M = 20 M⊙, XH = 0.2, Z = Z⊙, and D∞ = 10 for the optical wavelength regime with different colors corresponding to the different models.
In the text

Fig. 5 Mass-loss rate Ṁ as a function of Teff(τcrit) for our model sequences. For comparison, a set of empirically inferred temperatures T* for different types of WR stars is shown as well (gray symbols), illustrating the well-known “WR radius problem.” Since the empirical values are inferred from models without dynamical consistency, the values of T* defined at a Rosseland optical depth of τR,cont = 20, are shown. For our model sequences, the values of T* and Teff(τcrit) align very closely except for the highest mass-loss rates. A direct comparison with T* from the dynamically-consistent models is provided in Fig. E.1.

In the text

Fig. 6 Mass-loss rates versus different temperature scales for a series of dynamically consistent atmosphere models with log L/L⊙ = 5.7, M = 20 M⊙, and ХH = 0: the thick red dashed line denotes the effective temperatures defined at a Rosseland optical depth of τR = 2/3, while the green solid line and the blue dashed-dotted line denote the effective temperatures referring to τcrit and τR,cort = 20, respectively. The green dotted line on the right denotes the (electron) temperature at the critical point. Curves in lighter colors reflect models using the simple integration treatment suppressing negative velocity gradients (cf. Sect. 2). The black dashed line shows the hydrodynamic structure solutions by Grassitelli et al. (2018).

In the text

	Fig. 7 Same as Fig. 6, but for a model sequence with log L/L⊙ = 5.475, M = 15 M⊙, and ХH = 0. Contrary to the situation in Fig. 6 for a 20 M⊙ star, there is an abrupt breakdown of solutions beyond a minimum temperature and a maximum Ṁ. For the T2/з- and Teff(τcrit)-scales these points are marked with a vertical line attached to a gray-hatched area.
In the text

	Fig. 8 Mass-loss rates as a function of different temperature scales for the L/M model sequences presented in Sander & Vink (2020). Higher mass-loss rates correspond to higher L/M-ratios in this plot.
In the text

Fig. 9 Illustration of the radiative acceleration for a model with 10 M⊙ and T* ≈ 115 kK which is not capable of launching a wind from the “hot iron bump” and thus is dynamically not converged: In the inner part, the total radiative acceleration (red dashed line) approaches Γrad = 1, but does not surpass it sufficiently to launch a wind that could be maintained in the following deceleration region.

In the text

	Fig. 10 Transformed mass-loss rate as a function of T2/з for our model sequences.
In the text

Fig. 11 Effective temperature at a Rosseland optical depth of 2/3 as a function of the transformed mass-loss rate Ṁt for our new calculated model sequences. The sequences connected by solid lines all employ D∞ = 50, while those with dashed and dotted curves indicate sequences using D∞ = 10 and 4, respectively, as indicated in the plot. For comparison, also various empirical results from the literature are depicted by discrete, gray symbols.

In the text

Fig. 12 Effective temperature at a Rosseland optical depth of 2/3 as a function of the transformed mass-loss rate Ṁt for the whole set of models from Sander & Vink (2020). For log (Ṁt [M⊙ yr−1]) < −5.5, T2/3 is effectively independent of Ṁt. In Sander & Vink (2020), the value of T* is fixed for all models, but the difference in L/M still yields a wide range of T2/з values. The dashed-dotted line represents a linear fit of the temperature trend for log (Ṁt [M⊙ yr−1]) > −4.5 while the dotted curve represents the fit for the new model sequence illustrated in Fig. 11.

In the text

	Fig. 13 Mass-loss rate Ṁ as a function of the gravitational acceleration at the critical radius $g_{\text{crit}}=GM{R}_{\text{crit}}^{-2}$. Thin, dotted, gray lines indicate curves with $\dot{M}\propto g_{\text{crit}}^{-3/2}$.
In the text

	Fig. 14 Terminal velocity as a function of T2/3 for the different model sets. For comparison, also the derived trend for OB-star winds in the Milky Way (dashed line) and the LMC (dashed-dotted line) from Hawcroft et al. (in prep.) are shown.
In the text

Fig. 15 Transformed mass-loss rate ${\dot{M}}_{\text{t}}$ as a function of the gravitational acceleration at the critical radius $g_{\text{crit}}=GM{R}_{\text{crit}}^{-2}$. To reflect the expected trends from the Rcrit-fits, thin gray lines are plotted in the background. The dotted, gray lines indicate $\dot{M}\propto g_{\text{crit}}^{-2.5}$ (optically thin regime), while the dashed, gray lines correspond to $\dot{M}\propto g_{\text{crit}}^{-1.8}$ (optically thick regime).

In the text

	Fig. 16 Number of helium ionizing photons per second QHe II as a function of the transformed mass-loss rate ${\dot{M}}_{\text{t}}$ for the new model sequences calculated in this work.
In the text

	Fig. 17 Number of helium ionizing photons per second QHe II as a function of the transformed mass-loss rate ${\dot{M}}_{\text{t}}$ for the model sequences computed in Sander & Vink (2020).
In the text

	Fig. A.1 Linear fits (solid lines) to the Mass-loss rate $\dot{M}$ as a function of Teff(τcrit) in a double-logarithmic-plane. The fit coefficients for the different datasets are given in Table A.1.
In the text

	Fig. B.1 HRD with the effective temperature T2/з defined at a Rosseland optical depth of τRoss = 2/3 and the model luminosity L for our sets of He ZAMS models at different Z. For comparison, the HeZAMS (gray, dashed) and ZAMS (gray, solid) are shown as well.
In the text

Fig. C.1 Contributions of the different ions to the radiative acceleration of a hydrodynamically consistent, hydrogen-free WN model with T* = 130 kK, log L/L⊙ = 5.7, M = 20 M⊙, and D∞ = 50: Different ions are denoted by a combination of different color and symbol. The total radiative acceleration (arad), the Thomson acceleration from free electrons (aThom = Γe · g), and the contribution from gas (and turbulence) pressure (apress) are also shown for comparison. The loosely dashed horizontal line denotes the total Eddington limit that needs to be overcome to launch a wind.

In the text

	Fig. C.2 Contributions of different ions to the radiative acceleration, plotted similar to Fig. C.1, but now for a model employing D∞ = 10. The wind reaches a lower terminal velocity and the supersonic region with Γrad = arad/g < 1 is more pronounced than in the model with D∞ = 50.
In the text

Fig. D.1 Transformed mass-loss rate ${\dot{M}}_{\text{t}}$ as a function of L/M for different sequences varying in T*. All models use XH = 0.2 except the dashed-dotted sequences having XH = 0 for comparison. Apart from the red, dashed sequence (using D∞ = 10), all sequences employ D∞ = 50. The gray, dotted line is an interpolation of the solid sequences along the theoretical temperatures for the He ZAMS from Grassitelli et al. (2018) and Langer (1989). The light dotted linear curves in the background indicate the slope of 1.25 which is used in Eq. (D.4) and onward.

In the text

	Fig. E.1 Analogous plot to Fig. 5, but now showing the mass-loss rate $\dot{M}$ as a function of T* for our model sequences.
In the text

	Fig. E.2 Analogous plot to Fig. E.1, but now showing the transformed mass-loss rate ${\dot{M}}_{\text{t}}$ instead of the normal $\dot{M}$.
In the text

Fig. E.3 Mass-loss rates as a function of different temperature scales for a series of dynamically consistent atmosphere models with log L/L⊙ = 5.35, M = 12.9 M⊙, and XH = 0.2: The thick red dashed line denoted the classical effective temperatures defined at a Rosseland optical depth of τR = 2/3, while the green solid line and the blue dashed-dotted lines denotes the effective temperatures referring to τcrit and τR,cont = 20 respectively. The green dotted line on the right denotes the (electron) temperature at the critical point. Curves in lighter colors reflect models using the simple integration treatment suppressing negative velocity gradients.

In the text

	Fig. E.4 Analogous plot to Fig. E.3, but now for the WC model series with log L/L⊙ = 5.7, M = 20 M⊙, and 0.5 Z⊙.
In the text