Open Access
Issue
A&A
Volume 670, February 2023
Article Number A83
Number of page(s) 21
Section Stellar atmospheres
DOI https://doi.org/10.1051/0004-6361/202245110
Published online 09 February 2023

© The Authors 2023

Licence Creative CommonsOpen Access article, published by EDP Sciences, under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

This article is published in open access under the Subscribe to Open model. Subscribe to A&A to support open access publication.

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 arad 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 (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 (T2/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 T2/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 T2/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 (PoWRHD) 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 = arad g−1 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 (Rmax) 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−1, 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 (XH = 0.2) and one sequence where we switch to a WC-type (XC = 0.4, XN = 4 × 10−5, XO = 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 XH = 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.

thumbnail 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 υ.

thumbnail 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.

Table 1

Input parameters for our hydrodynamically consistent He-ZAMS models.

thumbnail 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.

Table 2

Overview of the calculated model sequences.

thumbnail 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 Teff(τ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 Teff(τcrit) (corresponding to decreasing radii Rcrit). 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 XH = 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 (Rcrit) 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 XH > 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 -Teff(τ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.

thumbnail 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.

thumbnail 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 T2/з at a Rosseland optical depth of τR = 2/3 (thick red dashed line), often simply denoted as Teff 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 Teff(τcrit), Teff(Rcrìt), or simply Teff,crit; and (iv) the electron temperature Te at the critical point (thin green dotted line). In contrast to the first three temperatures, Te is not an effective temperature. In the deeper layers of the atmosphere where the deviations from LTE become negligible, Te 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−4 M yr−1 in Fig. 6), the values of T* and Teff,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 Teff,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 (Teff,crit) and T2/з 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−1. 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е(Rcrit) 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е(Rcrit) 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е(Rcrit) increases mildly with lower until Teff,crit surpasses Tе(Rcrit).

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 T2/з-temperatures look very similar to those obtained in our new models with varying T*. On the other hand, Tе(Rcrit) varies much more than in any of our new model sequences. As we also see a shift in the Tе(Rcrit)-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е(Rcrit)-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е(Rcrit). In a zeroth-order approximation, one could state that Tе(Rcrit) is constant for a given L/M and chemical composition.

thumbnail 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.

thumbnail 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 Teff,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 Teff,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).

thumbnail 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.

Table 3

Comparison of mass-loss rates obtained with different methods for a hydrogen-free WN star with log L/L = 5.7 and 20 M⊙.

3.3 Scaling with the transformed mass-loss rate

Since the different calculated model sequences show very similar slopes for (T2/з), 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” (2)

originally introduced by Gräfener & Vink (2013), for our new model sequences as a function of T2/з. 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 . For our 20 M models with XH = 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 (T2/з < 130kK and log t > −4.5) can be reasonably well approximated by a linear fit, yielding (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 .

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 T2/з 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 T2/з 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 T2/з values, but employing more compact radii that better align with structural predictions.

Despite the generally better match when comparing T2/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 T2/з 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 T2/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 T2/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 (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.

thumbnail Fig. 10

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

thumbnail 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.

thumbnail 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 Teff,crit for our sequences with fixed L/M. Performing a linear fit, we obtain a relation in the form of (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 and (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 Rcrit-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 (7)

we can rewrite Eq. (6) as (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 gcrit is not surprising as an increased gravitational force needs to be overcome. Given the constant mass M along the model sequences, the change in gcrit expected from Eq. (8) is purely geometrical, that is only from the change in Rcrit. 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 Te(Rcrit) surpassing Teff,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 gcrit has to be considered as well, which we do in the following section.

thumbnail Fig. 13

Mass-loss rate as a function of the gravitational acceleration at the critical radius . Thin, dotted, gray lines indicate curves with .

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 T2/з 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).

thumbnail 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 T2/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 (XH = 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 atmosphere1. 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 T2/3 decreases when comparing models with the same T*. Interestingly, as we saw in Fig. E.1, the value of remains the same when comparing against T2/3. In a follow-up study, we will test whether this behavior is universal when considering surface hydrogen or whether the chosen fraction of XH = 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 in Fig. 5.

To get some insights on the general scaling of WR-type winds with effective temperature (here: T2/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 (T2/3 > 130 kK), when the winds become more optically thin as the mass-loss rates decrease, the steepness of the v(T2/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 T2/3, we further investigate the behavior of v as a function of T* or the physically more meaningful Teff,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 (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 can also be obtained from directly fitting v(Rcrit) in the optically thick limit. Using the critical radius to define the escape velocity (10)

we obtain (11)

This relation also holds for the effective escape velocity (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 Rcrit in the considered temperature range. Thus, in the optically thick wind regime we have a slight difference with along the T* -dimension compared to the well-known vvesc,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 vesc,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 vesc,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 vesc.

For lower mass-loss rates (log ), corresponding usually to Teff,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(T2/3), with the slopes known from OB-type winds. With considerable scatter in the exponent of up to ±0.5, we find (13)

corresponding to or , 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 vesc,eff is likely more complicated.

4.4 Influence on Mt

The scaling of v with Rcrit introduces an additional dependency when considering the transformed mass-loss rate as a function of Rcrit or gcrit.

From Eqs. (5) and (6), we know that or . From the definition of (Eq. (2)) we get (14)

With the different trends derived for the optically thick and thin limit in Sect. 4.3, we obtain and respectively. Using gcrit as defined in Eqs. (7) and (8), this yields (15)

for the optically thick limit and (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 -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 should be considered less robust than the -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 . While it is hard to draw strong conclusions, there at least appears to be one continuous slope for in the thinner wind regime, regardless of whether the critical point is located at temperatures below or above Teff,crit. Since we find a change for 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 Teff,crit.

thumbnail Fig. 15

Transformed mass-loss rate as a function of the gravitational acceleration at the critical radius . To reflect the expected trends from the Rcrit-fits, thin gray lines are plotted in the background. The dotted, gray lines indicate (optically thin regime), while the dashed, gray lines correspond to (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 Teff,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.12, 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 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 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 . 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 Teff,critT* ≈ 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 -values from Gräfener et al. (2017) and commonly used (semi-)empirical recipes. For Z we find a difference of ~0.2 dex in , 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 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 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 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 XH = 0.2 and 0.5 Z extends to much cooler temperatures (Teff,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 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 . 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 toward lower He star masses. On the other hand, we reach lower limits of Teff,crit for driving winds by the hot iron bump (cf. Sect. 3.2). In our 12.9 M sequence with XH = 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 LM 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 . Since Eq. (5) seems to be valid across both the optically thick and thin regime, we can approximate for WN winds driven by the hot iron bump via (17)

with denoting the mass-loss rate from Sander & Vink (2020) and Rcrit,T141 = Rcrit (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 Rcrit,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 Teff,crit, yielding (18)

Apart from small deviations for the highest mass-loss rates , the fixed value of 141 kK accurately represents the value of Teff,crit in Sander & Vink(2020), as illustrated previously in Fig. 8.

This adjusted -recipe requires the knowledge of either Teff,crit or Rcrit. As we did include only a small microturbulent velocity in our modeling efforts (30 km s−1), 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 (T2/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 Rmax of our atmosphere models to extremely large values, often up to 100000 R*. (Typical atmosphere models for spectral fitting require only Rmax = 100 … 1000 R*.)

In Fig. 16, we plot the rate of ionizing photons per second QHe II as a function of the transformed mass-loss rate . The absolute numbers in the regime with low QHe 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 Rmax, the value for QHe II can be underestimated, but should never exceed 1041 s−1. Given that this is many orders of magnitude below the actually strong QHe II emitters with rates >1047 s−1, the values of the shaded regime in Fig. 16 should have no practical consequences.

Displaying the He II ionizing flux as a function of 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 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 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.

thumbnail Fig. 16

Number of helium ionizing photons per second QHe II as a function of the transformed mass-loss rate for the new model sequences calculated in this work.

thumbnail Fig. 17

Number of helium ionizing photons per second QHe II as a function of the transformed mass-loss rate 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 (, v) as well as the characteristics in the outer wind were maintained. We draw the following conclusions:

  • The mass-loss rates depend significantly on the critical radius Rcrit and thus also on the assumed model temperature setting Teff(Rcrit). For model sequences with constant luminosity L and stellar mass M, we obtain 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 , 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 and above – the effective temperature at the critical point Teff(τ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 for the transition between the optically thin and thick regime. While there is some scatter between different model sequences, this characteristic value of (plus some error margin) provides a very convenient tool to distinguish between the regimes as 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 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 work3. Their corresponding transformed mass-loss rates appear to be below −4.5, e.g. at 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 trends, but leads to an additional shift in the obtained relations with higher clumping factors corresponding to higher for the same T2/3. In contrast to all shifts introduced by varying fundamental stellar parameters or abundances, the shift due to D does not vanish when considering instead of .

  • In the limit of optically thick winds, we obtain a linear relation between log T2/3 and log , independent of chemical composition (but for a fixed clumping factor). Combined with the also Z-independent result from Sander & Vink (2020) that log , 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 are strong emitters of He II ionizing flux (with QHe II > 1048 s−1). 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 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 T2/3--plane to our derived curves, we find a significant fraction of stars to have lower values of 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 Te ≈ 2 × 106 K – wind quantity predictions such as 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 appears to break down suddenly. We provide an adjustment of the recent -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 -temperature Fit

For all of our model sequences, the data points in the log -log Teff(τcrit)-plane suggest a linear relation between the two quantities with deviations occurring only close to the wind driving limit (corresponding to the maximum in Fig. A.1). In the fits, we thus exclude the uppermost 0.15 dex in log . 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 Teff(τcrit) > Teff,crit,min and by (A.1)

with the value for being different for each model sequence. Table A.1 also lists these coefficients together with the corresponding validity limits Teff,crit,min and .

thumbnail Fig. A.1

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

Table A.1

– Linear fit results for log versus log Teff,crit plus offsets and limitations for the temperature-dependent mass loss of our model sequences described by Eq. (A.1).

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

The effective temperatures T2/з 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 T2/з 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 T2/з to higher l and steepening in particular this drop. The differences in T2/з 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.

thumbnail 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 (C.1)

similar integrals to Eq. (C.1) are calculated using only the ion-specific opacities (e.g. ) 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 vcl = 100 km s−1, 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−1 to 1754 km s−1.

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.

Table C.1

Derived wind parameters for WN models with different D

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 vcl. 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−1. This can then have two effects leading to a lower for lower values of D: first, a direct reduction of opacities in the region that determines . 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 .

In the last column of Table C.1, we provide the resulting transformed mass-loss rates . While there is already a scaling with 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 ). For example, the hydrogen free models differ by , while v only increases by a factor of 1.48. Therefore, the models with higher D posses a higher , despite larger terminal velocities reducing its value.

thumbnail 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.

thumbnail 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 T2/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 -T2/3-plane, yielding (D.1)

for log , thereby providing us with a potential path to predict the observable effective temperature. The relation (D.1) seems to be approximately unaffected by abundance (Xi) changes, but there is a clear offset for different choices of D. In our calculations, there is a difference of ∆ log (T2/3 [K]) ≈ 0.08 … 0.1 between D = 4 and D = 10 and ∆ log (T2/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 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 and v is only reached in the optically thin part, while we obtained in the dense wind regime. As a consequence, models from the two different sources with approximately the same value of 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 will be higher. Arguing that is the major factor in setting the T2/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 , we therefore suggest the more shallow formula (D.2)

as a first attempt to approximate the observable effective temperature T2/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 . 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 Mt needs to be known to determine T2/3. In Sander & Vink (2020), we could show that for T* = 140 kK the quantity is practically independent of metallicity in the limit of optically thick winds (“pure WR regime”). There, 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 Teff,critT*. 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 Teff,critT*. 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 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 and obtain (D.3)

or in the regime of optically thick winds, i.e. for log .

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 (D.4)

Again, it might be more convenient to replace Rcrit with Teff,crit and gauge this with the 20 M model at 141 kK. This then yields (D.5)

thumbnail Fig. D.1

Transformed mass-loss rate 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 which contains further, uncertain dependencies. These can alter the resulting values of , 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 T2/3: (D.6)

This formula is only valid for hydrogen-free WN stars as we have considerable offsets for other chemical compositions in (cf. Fig. D.1). While the conversion between and T2/з 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 and υ and the resulting -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 and 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 T2/з 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 due to a larger L/M. Nonetheless, we can present the first estimate of T2/з 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 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 Teff,crit.

thumbnail Fig. E.1

Analogous plot to Fig. 5, but now showing the mass-loss rate as a function of T* for our model sequences.

thumbnail Fig. E.2

Analogous plot to Fig. E.1, but now showing the transformed mass-loss rate instead of the normal .

thumbnail 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.

thumbnail 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 instead of . In Fig. E.2, we show the analogous plot to Fig. E.1 with being replaced by . 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, XH = 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

  1. Aadland, E., Massey, P., John Hillier, D., et al. 2022, ApJ, 931, 157 [NASA ADS] [CrossRef] [Google Scholar]
  2. Abbott, B. R., Abbott, R., Abbott, T. D., et al. 2016, ApJ, 818, L22 [NASA ADS] [CrossRef] [Google Scholar]
  3. Bailey, J. E., Nagayama, T., Loisel, G. P., et al. 2015, Nature, 517, 56 [NASA ADS] [CrossRef] [Google Scholar]
  4. Bestenlehner, J. M. 2020, MNRAS, 493, 3938 [Google Scholar]
  5. Bestenlehner, J. M., Crowther, P. A., Caballero-Nieves, S. M., et al. 2020, MNRAS, 499, 1918 [Google Scholar]
  6. Björklund, R., Sundqvist, J. O., Puls, J., & Najarro, F. 2021, A&A, 648, A36 [EDP Sciences] [Google Scholar]
  7. Castor, J. I., Abbott, D. C., & Klein, R. I. 1975, ApJ, 195, 157 [Google Scholar]
  8. Chieffi, A., & Limongi, M. 2013, ApJ, 764, 21 [Google Scholar]
  9. Conti, P. S. 1991, ApJ, 377, 115 [NASA ADS] [CrossRef] [Google Scholar]
  10. Crowther, P. A., & Hadfield, L. J. 2006, A&A, 449, 711 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  11. Crowther, P.A., Dessart, L., Hillier, D.J., Abbott, J.B., & Fullerton, A.W. 2002, A&A, 392, 653 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  12. de Koter, A., Heap, S. R., & Hubeny, I. 1997, ApJ, 477, 792 [Google Scholar]
  13. Dray, L. M., Tout, C. A., Karakas, A. I., & Lattanzio, J. C. 2003, MNRAS, 338, 973 [NASA ADS] [CrossRef] [Google Scholar]
  14. Farmer, R., Laplace, E., de Mink, S. E., & Justham, S. 2021, ApJ, 923, 214 [NASA ADS] [CrossRef] [Google Scholar]
  15. Garcia, M., Herrero, A., Najarro, F., Lennon, D. J., & Alejandro Urbaneja, M. 2014, ApJ, 788, 64 [NASA ADS] [CrossRef] [Google Scholar]
  16. Georgy, C., Ekström, S., Meynet, G., et al. 2012, A&A, 542, A29 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  17. Gräfener, G., & Hamann, W.-R. 2005, A&A, 432, 633 [CrossRef] [EDP Sciences] [Google Scholar]
  18. Gräfener, G., & Vink, J. S. 2013, A&A, 560, A6 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  19. Gräfener, G., Koesterke, L., & Hamann, W.-R. 2002, A&A, 387, 244 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  20. Gräfener, G., Vink, J. S., de Koter, A., & Langer, N. 2011, A&A, 535, A56 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  21. Gräfener, G., Vink, J. S., Harries, T. J., & Langer, N. 2012, A&A, 547, A83 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  22. Gräfener, G., Owocki, S. P., Grassitelli, L., & Langer, N. 2017, A&A, 608, A34 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  23. Grassitelli, L., Langer, N., Grin, N. J., et al. 2018, A&A, 614, A86 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  24. Groh, J. H., Meynet, G., Ekström, S., & Georgy, C. 2014, A&A, 564, A30 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  25. Hainich, R., Rühling, U., Todt, H., et al. 2014, A&A, 565, A27 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  26. Hainich, R., Pasemann, D., Todt, H., et al. 2015, A&A, 581, A21 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  27. Hamann, W. R. 1985, A&A, 145, 443 [NASA ADS] [Google Scholar]
  28. Hamann, W.-R., & Gräfener, G. 2003, A&A, 410, 993 [CrossRef] [EDP Sciences] [Google Scholar]
  29. Hamann, W.-R., & Gräfener, G. 2004, A&A, 427, 697 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  30. Hamann, W., & Koesterke, L. 1998, A&A, 335, 1003 [NASA ADS] [Google Scholar]
  31. Hamann, W. R., Koesterke, L., & Wessolowski, U. 1995, A&A, 299, 151 [Google Scholar]
  32. Hamann, W., Gräfener, G., & Liermann, A. 2006, A&A, 457, 1015 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  33. Hamann, W.-R., Gräfener, G., Liermann, A., et al. 2019, A&A, 625, A57 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  34. Higgins, E. R., Sander, A. A. C., Vink, J. S., & Hirschi, R. 2021, MNRAS, 505, 4874 [NASA ADS] [CrossRef] [Google Scholar]
  35. Hillier, D. J. 1991, in Wolf-Rayet Stars and Interrelations with Other Massive Stars in Galaxies, eds. K.A. van der Hucht, & B. Hidayat, 143, 59 [NASA ADS] [CrossRef] [Google Scholar]
  36. Hillier, D. J., & Miller, D. L. 1999, ApJ, 519, 354 [Google Scholar]
  37. Hillier, D. J., Davidson, K., Ishibashi, K., & Gull, T. 2001, ApJ, 553, 837 [NASA ADS] [CrossRef] [Google Scholar]
  38. Köhler, K., Langer, N., de Koter, A., et al. 2015, A&A, 573, A71 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  39. Langer, N. 1989, A&A, 210, 93 [NASA ADS] [Google Scholar]
  40. Langer, N., Hamann, W. R., Lennon, M., et al. 1994, A&A, 290, 819 [NASA ADS] [Google Scholar]
  41. Lefever, R. R., Shenar, T., Sander, A. A. C., et al. 2022, ArXiv e-prints [arXiv:2209.06043] [Google Scholar]
  42. Leitherer, C., Vacca, W. D., Conti, P.S., et al. 1996, ApJ, 465, 717 [NASA ADS] [CrossRef] [Google Scholar]
  43. Limongi, M., & Chieffi, A. 2018, ApJS, 237, 13 [NASA ADS] [CrossRef] [Google Scholar]
  44. Maeder, A. 1983, A&A, 120, 113 [NASA ADS] [Google Scholar]
  45. Martinet, S., Meynet, G., Nandal, D., et al. 2022, A&A, 664, A181 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  46. Martins, F., Hillier, D. J., Paumard, T., et al. 2008, A&A, 478, 219 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  47. Moens, N., Poniatowski, L. G., Hennicker, L., et al. 2022, A&A, 665, A42 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  48. Moriya, T.J., & Yoon, S.-C. 2022, MNRAS, 513, 5606 [Google Scholar]
  49. Nakauchi, D., & Saio, H. 2018, ApJ, 852, 126 [NASA ADS] [CrossRef] [Google Scholar]
  50. Nugis, T., & Lamers, H. J. G. L. M. 2000, A&A, 360, 227 [NASA ADS] [Google Scholar]
  51. Nugis, T., & Lamers, H. J. G. L. M. 2002, A&A, 389, 162 [EDP Sciences] [Google Scholar]
  52. Petrovic, J., Pols, O., & Langer, N. 2006, A&A, 450, 219 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  53. Plat, A., Charlot, S., Bruzual, G., et al. 2019, MNRAS, 490, 978 [Google Scholar]
  54. Poniatowski, L. G., Sundqvist, J. O., Kee, N. D., et al. 2021, A&A, 647, A151 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  55. Ro, S. 2019, ApJ, 873, 76 [NASA ADS] [CrossRef] [Google Scholar]
  56. Ro, S., & Matzner, C. D. 2016, ApJ, 821, 109 [Google Scholar]
  57. Sander, A. A. C., & Vink, J. S. 2020, MNRAS, 499, 873 [Google Scholar]
  58. Sander, A., Hamann, W.-R., & Todt, H. 2012, A&A, 540, A144 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  59. Sander, A., Shenar, T., Hainich, R., et al. 2015, A&A, 577, A13 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  60. Sander, A. A. C., Hamann, W.-R., Todt, H., Hainich, R., & Shenar, T. 2017, A&A, 603, A86 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  61. Sander, A. A. C., Fürst, F., Kretschmar, P., et al. 2018, A&A, 610, A60 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  62. Sander, A. A. C., Hamann, W.-R., Todt, H., et al. 2019, A&A, 621, A92 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  63. Sander, A. A. C., Vink, J. S., & Hamann, W. R. 2020, MNRAS, 491, 4406 [Google Scholar]
  64. Schaerer, D., Contini, T., & Kunth, D. 1999, A&A, 341, 399 [NASA ADS] [Google Scholar]
  65. Schmutz, W., Leitherer, C., & Gruenwald, R. 1992, PASP, 104, 1164 [NASA ADS] [CrossRef] [Google Scholar]
  66. Shaviv, N. J. 1998, ApJ, 494, L193 [NASA ADS] [CrossRef] [Google Scholar]
  67. Shaviv, N. J. 2000, ApJ, 532, L137 [NASA ADS] [CrossRef] [Google Scholar]
  68. Shenar, T., Hainich, R., Todt, H., et al. 2016, A&A, 591, A22 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  69. Shenar, T., Sablowski, D. P., Hainich, R., et al. 2019, A&A, 627, A151 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  70. Shenar, T., Gilkis, A., Vink, J. S., Sana, H., & Sander, A. A. C. 2020, A&A, 634, A79 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  71. Smith, L. J., Norris, R. P. F., & Crowther, P. A. 2002, MNRAS, 337, 1309 [CrossRef] [Google Scholar]
  72. The LIGO Scientific Collaboration, The Virgo Collaboration, The KAGRA Collaboration. 2021, ArXiv e-prints [arXiv:2111.03634] [Google Scholar]
  73. Vink, J. S., & de Koter, A. 2005, A&A, 442, 587 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  74. Vink, J. S., & Gräfener, G. 2012, ApJ, 751, L34 [Google Scholar]
  75. Vink, J. S., & Sander, A. A. C. 2021, MNRAS, 504, 2051 [NASA ADS] [CrossRef] [Google Scholar]
  76. Vink, J. S., Higgins, E. R., Sander, A. A. C., & Sabhahit, G. N. 2021, MNRAS, 504, 146 [NASA ADS] [CrossRef] [Google Scholar]
  77. Woosley, S. E., Sukhbold, T., & Janka, H. T. 2020, ApJ, 896, 56 [NASA ADS] [CrossRef] [Google Scholar]
  78. Yoon, S.-C. 2017, MNRAS, 470, 3970 [NASA ADS] [CrossRef] [Google Scholar]
  79. Yoon, S. C., Langer, N., & Norman, C. 2006, A&A, 460, 199 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  80. Yusof, N., Hirschi, R., Eggenberger, P., et al. 2022, MNRAS, 511, 2814 [NASA ADS] [CrossRef] [Google Scholar]

1

A small hydrogen-layer on the surface has also the structural consequence of an increased stellar radius, which would again affect the wind parameters. Here, we discuss only the immediate atmospheric consequences for a fixed set of stellar parameters.

2

In this discussion, we do not consider structure models that show hydrostatic envelope inflation for more massive He stars (see, e.g., Fig. 19 in Köhler et al. 2015) as Grassitelli et al. (2018) demonstrated that such an inflation likely does not occur if a strong wind can be launched.

3

The transformed mass-loss rate as such should not be confused with the transition mass-loss rate from Vink & Gräfener (2012). Nonetheless, if the other necessary parameters are known, one can estimate the corresponding transformed mass-loss rates for stars defining the transition mass-loss rate, denoted as .

All Tables

Table 1

Input parameters for our hydrodynamically consistent He-ZAMS models.

Table 2

Overview of the calculated model sequences.

Table 3

Comparison of mass-loss rates obtained with different methods for a hydrogen-free WN star with log L/L = 5.7 and 20 M⊙.

Table A.1

– Linear fit results for log versus log Teff,crit plus offsets and limitations for the temperature-dependent mass loss of our model sequences described by Eq. (A.1).

Table C.1

Derived wind parameters for WN models with different D

All Figures

thumbnail 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
thumbnail 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
thumbnail 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
thumbnail 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
thumbnail 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
thumbnail 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
thumbnail 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
thumbnail 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
thumbnail 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
thumbnail Fig. 10

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

In the text
thumbnail 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
thumbnail 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
thumbnail Fig. 13

Mass-loss rate as a function of the gravitational acceleration at the critical radius . Thin, dotted, gray lines indicate curves with .

In the text
thumbnail 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
thumbnail Fig. 15

Transformed mass-loss rate as a function of the gravitational acceleration at the critical radius . To reflect the expected trends from the Rcrit-fits, thin gray lines are plotted in the background. The dotted, gray lines indicate (optically thin regime), while the dashed, gray lines correspond to (optically thick regime).

In the text
thumbnail Fig. 16

Number of helium ionizing photons per second QHe II as a function of the transformed mass-loss rate for the new model sequences calculated in this work.

In the text
thumbnail Fig. 17

Number of helium ionizing photons per second QHe II as a function of the transformed mass-loss rate for the model sequences computed in Sander & Vink (2020).

In the text
thumbnail Fig. A.1

Linear fits (solid lines) to the Mass-loss rate 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
thumbnail 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
thumbnail 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
thumbnail 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
thumbnail Fig. D.1

Transformed mass-loss rate 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
thumbnail Fig. E.1

Analogous plot to Fig. 5, but now showing the mass-loss rate as a function of T* for our model sequences.

In the text
thumbnail Fig. E.2

Analogous plot to Fig. E.1, but now showing the transformed mass-loss rate instead of the normal .

In the text
thumbnail 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
thumbnail 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

Current usage metrics show cumulative count of Article Views (full-text article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.

Data correspond to usage on the plateform after 2015. The current usage metrics is available 48-96 hours after online publication and is updated daily on week days.

Initial download of the metrics may take a while.