© The Authors 2026

1 Introduction

Winds from massive O-type stars in the range of 20-60 M_⊙ are a key ingredient in stellar evolutionary models, shaping both their structure and their final fate. However, for the most massive stars in the Universe with initial masses above 100 M_⊙ - the so-called very massive stars (VMSs) - mass loss is not merely important: It is the dominant factor that determines their evolution and ultimate fate (Hirschi 2015; Köhler et al. 2015).

O-star winds have long been described by the optically thin line-driving theory of Lucy & Solomon (1970) and Castor et al. (1975, hereafter CAK). Evolutionary models for VMSs have traditionally extrapolated these prescriptions into higher masses. However, such extrapolations likely underestimate the true VMS mass-loss rates. Vink et al. (2011) using Monte Carlo (MC) simulations predicted an upturn or a ‘kink’ in the mass-loss-luminosity relation, associated with the Eddington parameter Γ ∝ L/M. This kink suggests that winds strengthen dramatically and become optically thick; that is, the wind optical depth, τ_F, crosses unity, as stars approach the Eddington limit.

Concurrently with these VMS wind predictions, the widely assumed stellar upper-mass limit of 100-150 M_⊙ (Weidner & Kroupa 2004; Figer 2005; Oey & Clarke 2005) was revised to ~200-300 M_⊙ (Crowther et al. 2010; Bestenlehner et al. 2011; Martins 2015; Brands et al. 2022; Kalari et al. 2022). These values, however, are based on present-day masses extrapolated to the zero-age main sequence (ZAMS) using O-star rates, and may thus be underestimated. If VMS mass loss is indeed stronger than CAK extrapolations suggest, the intrinsic upper-mass limit may extend significantly higher - perhaps up to 1000 M_⊙ (Vink 2018; Sabhahit et al. 2022; Higgins et al. 2022). A closely related issue concerns the formation of VMSs (Krumholz 2015) and the upper-mass limit, which may be regulated by the strength of stellar winds in environments of different metallicity.

In the local Universe, the hydrogen-rich class of Wolf-Rayet stars of the nitrogen (N) sequence (WNh stars) in young (~1-2 Myr) massive clusters, which exhibit strong He II emission, are our best candidates for VMSs (Massey & Hunter 1998; de Koter et al. 1998; Martins et al. 2008; Crowther et al. 2016; Wofford et al. 2014; Upadhyaya et al. 2024). Such VMSs are not only key to stellar evolution theory but also to astrophysical environments, since they may drive the chemical enrichment of globular clusters (Vink 2018; Higgins et al. 2023, 2025; Gieles et al. 2025) and also provide crucial feedback at high redshift (Senchyna et al. 2024; Vink 2023; Charbonnel et al. 2023; Cameron et al. 2023; Isobe et al. 2023; Topping et al. 2024; Curti et al. 2025; Ji et al. 2025). Predicting their wind yields and final masses requires reliable mass-loss prescriptions, especially across the transition from optically thin (τ_F < 1) O-star winds to optically thick (τ_F > 1) WNh winds (Vink & Gräfener 2012).

The Vink et al. (2011) kink approximately coincides with the wind efficiency parameter η = Ṁv_∞c/L approaching unity. Above this threshold, the predicted mass-loss scaling with Γ steepens from approximately two to about five, far stronger than a CAK-style extrapolation. While CMFGEN analyses from the VLT-FLAMES Tarantula Survey (Bestenlehner et al. 2014) provide empirical support, several studies still neglect the η ≈ τ_F ≈ 1 criterion emphasised by Vink & Gräfener (2012).

Given the upturn or kink predicted by the MC simulations, wind calculations using a different method and code are important to independently verify the MC results, which motivates new predictions of VMS mass-loss rates both below and above the potential upturn or kink. Here we present new computations with the hydrodynamically consistent branch of the Postdam Wolf-Rayet (PoWR) code (Sander et al. 2017), as recently applied to the study of late-type WNh stars (Lefever et al. 2025) and some of the most massive stars in the R136 cluster in the Large Magellanic Cloud (LMC) (Sabhahit et al. 2025). These kinds of hydro-models are also ideal tools to study the transition from absorption-line dominated O-star spectra to emission-line dominated WNh spectra as the objects approach the Eddington limit. Along this transition, in particular, the equivalent width of the Wolf-Rayet (WR) He II emission feature increases (see Vink et al. 2011; Crowther & Walborn 2011).

The paper is organised as follows: in Sect. 2, we detail the basics of our hydrodynamical wind-atmosphere models and the necessary inputs and outputs. In Sects. 3 and 4, we present the radiative acceleration and mass-loss properties of VMSs, discussing the mass-loss kink feature and bistable behaviour at cooler temperatures. We further test our mass-loss predictions against the transition mass-loss rate in the Arches Cluster. A brief comparison to previous VMS wind investigations is presented in Sect. 5, and we conclude in Sect. 6.

2 PoWR^HD atmosphere models

We utilised the hydrodynamical branch of the non-local thermodynamic equilibrium (non-LTE) stellar atmospheric code PoWR (Gräfener et al. 2002; Hamann & Gräfener 2003; Gräfener & Hamann 2005, 2008; Sander 2015; Sander et al. 2015, 2017) to build our model atmosphere grid. PoWR models a onedimensional, spherically symmetric, expanding atmosphere by solving the equations of radiative transfer, statistical and thermal equilibrium, and continuity. The atmosphere computational framework has long relied on the Accelerated Lambda Operator iteration (for details, see Hamann 1985, 1986) to solve these fully coupled non-linear equations across both spatial and frequency domains. Radiative transfer calculations are performed using the co-moving frame approach, in which the opacity and emissivity profile functions become isotropic (Mihalas et al. 1975). The code takes relevant stellar and wind parameters as inputs and computes the emergent synthetic spectra. The input stellar parameters include luminosity, temperature (or radius), logg (or mass), and chemical abundances (for a detailed input parameters, see Sect. 2.1). The wind parameters are specified by the mass-loss rate, Ṁ, and the velocity stratification, v(r), which is typically prescribed using a standard β-type velocity law of the form v(r) ~ (1 - R_*/r)^β.

The fundamental advance and novelty of the PoWR^HD models over the base code stems from solving the stationary hydrodynamic equation of motion (cf. Eq. (1)). The code balances the time-independent advective acceleration term (υdυ/dr), against the various forces in the wind-atmosphere, such as gravity and the gas and radiation pressure gradients at each atmospheric layer, thereby integrating hydrodynamics into the overall iterative procedure of the atmosphere code. By solving the wind hydrodynamics, the models gain predictive power, as the mass-loss rate and the wind velocity stratification are no longer prescribed inputs. These quantities instead emerge as predicted outputs, calculated self-consistently from the atmospheric force balance at each iteration step.

Such next-generation hydrodynamically consistent atmosphere models therefore achieve what earlier codes could not: they simultaneously derive consistent wind structure and mass loss, and generate synthetic spectra through unified non-LTE atmosphere-hydrodynamics calculations, enabling empirical determination of fundamental stellar and wind parameters (see e.g. Sander et al. 2017; Sabhahit et al. 2025). The numerical implementation details and general input specifications are documented in detail in Sander (2015); Sander et al. (2017, 2020). Below, we briefly summarise the required inputs to build our model grid and the relevant outputs.

2.1 Key input parameters

In this work, we isolated the effects of stellar mass and inner boundary temperature on VMS wind properties. While a comprehensive parameter space exploration and theoretical massloss predictions for the most-massive WNh stars based on PoWR^HD models will follow in a forthcoming work, the present study maintains all parameters except mass and temperature at the fixed values detailed below.

We fixed the luminosity to a constant value of log(L_*/L_⊙) = 6.0 throughout our grid. The inner boundary radius, R_*, of our models is defined at a specific Rosseland continuum optical depth, τ_R,cont, of 20. The effective temperature, T_*, at the inner boundary is varied between 12-50 kK, primarily encompassing the temperature range characteristic of observed WNh stars and luminous blue variables (LBVs) at this luminosity. The stellar mass range extended from approximately 40-135 M_⊙, with the precise range depending on T_*¹. The outer boundary was set to R_out = 1000 R_*.

For the chemical abundances, we adopted ZAMS values following solar-scaled patterns from Grevesse & Sauval (1998), with a total metal mass fraction of Z = 0.02 (Appendix B details the mass fraction breakdown for individual elements, along with levels and lines per ion). The hydrogen (H) mass fraction was fixed at X = 0.7, while helium (He) followed from the unity constraint X + Y + Z = 1.

Density inhomogeneities - ubiquitous in massive star winds - were incorporated through the micro-clumping formalism of Hillier et al. (2003), with the clumping onset occurring at a specified Rosseland optical depth (see Sander et al. 2017, for the implementation). The clumping factor, D_cl, was increased from unity (smooth wind) to an asymptotic value, D_cl,∞ = 10, with the onset at Rosseland optical depth τ_R = 0.1. This choice of onset is motivated by multi-wavelength spectral fits to archetypal VMSs in the Tarantula Nebula, specifically R136a1 (WN5h) and the binary system R144 (WN5/6h+WN6/7h) using PoWR^HD models (Sabhahit et al. 2025). Turbulent pressure contributions were included through adoption of a radially constant turbulent velocity υ_turb = 70km s⁻¹. Such high values are motivated from time-dependent, 2D simulations of massive O-star atmospheres, where substantial photospheric turbulence of the order 30-100 km s⁻¹ is expected to originate from the hot iron (Fe) bump (Debnath et al. 2024).

2.2 Output quantities

The PoWR^HD models enforce the stationary hydrodynamic equation of motion throughout the atmosphere: $$\begin{array}{r}\upsilon (r){\displaystyle \frac{\mathrm{d}\upsilon (r)}{\mathrm{d}r}}=-{\displaystyle \frac{1}{\rho (r)}}{\displaystyle \frac{\mathrm{d}P(r)}{\mathrm{d}r}}-{\displaystyle \frac{G{M}_{\star }}{r^2}}+a_{\mathrm{r}\mathrm{a}\mathrm{d}}(r),\end{array}$$(1)

where P(r) is the combined gas and turbulent pressure, while ρ(r) and υ(r) denote the mass density and velocity stratifications which are linked to the mass-loss rate Ṁ through the continuity equation: Ṁ = 4πr²ρ(r)υ(r). The total radiative acceleration a_rad(r), arising from momentum transfer on transition lines, scattering, and other continuum processes, fundamentally couples the hydrodynamic equation to the radiation field and population numbers.

Equation (1) fundamentally represents the force balance equation, with the advective acceleration term on the left-hand side counterbalanced by the cumulative atmospheric forces on the right - gravity and the gradients of gas, radiation, and turbulent pressure. The apparent simplicity of Eq. (1) is however misleading. The complexity lies masked within the radiative acceleration term a_rad(r), given by $$\begin{array}{r}{a}_{\mathrm{r}\mathrm{a}\mathrm{d}}(r)={\displaystyle \frac{1}{c}}{\int }_0^{\mathrm{\infty }}{\aleph }_{\nu }(r)F_{\nu }(r)\mathrm{d}\nu ={\displaystyle \frac{{\aleph }_{\mathrm{F}}(r)F(r)}{c}},\end{array}$$(2)

where ϗ_F(r) and F(r) represent the flux-weighted mean opacity and total flux at radius r, respectively. The quantity ϗ_F(r) cannot be determined locally; it must instead incorporate the detailed frequency-dependent opacity, which depends on radiation and velocity fields, along with ionisation and excitation states of absorbing ions throughout the atmosphere. Such complexity demands a brute-force, iterative global integration to achieve self-consistent hydrodynamic atmosphere models.

The flux-weighted mean opacity calculation is however essential because it accounts for Doppler effects in the expanding wind and ensures accurate radiative acceleration through conservation of the flux integral in Eq. (2). Unlike the classical Thomson acceleration - well constrained at model initialisation, particularly when stellar temperatures are high enough for complete H and He ionisation - the total radiative acceleration is calculated iteratively. The radial stratification of this total radiative acceleration from converged models is examined in Sect. 3.

Using the continuity equation in Eq. (1) and rearranging the hydrodynamic equation such that all terms containing the velocity gradient are grouped together on one side, we can show that the hydrodynamic equation features a critical point at radius R_crit (see derivation in Sander et al. 2017). The critical point in our model occurs where the flow velocity equals the isothermal sound speed corrected for turbulent velocity², so that the critical radius R_crit roughly coincides with the sonic point radius R_sonic. Each hydrodynamic update recalculates the velocity stratification throughout the atmosphere via inward and outward integration from R_crit. We therefore predict the complete velocity stratification υ(r), including the terminal velocity at the outer boundary: υ_∞ = υ(r = R_out). Mass-loss rate updates preserve the inner boundary condition to maintain the specified T_R,cont.

From the predicted Ṁ and υ_∞, the wind efficiency parameter is given by: $$\eta ={\displaystyle \frac{\dot{M}{\upsilon }_{\mathrm{\infty }}}{L_{\star }/c}}.$$(3)

The parameter η quantifies the ratio of wind momentum to available photon momentum. For powerful winds characteristic of WNh and classical WR stars, where multiple photon scattering becomes significant, η can easily exceed the single-scattering limit of unity.

Another important wind-related quantity outputted by our models is the wind optical depth, specifically, the flux-weighted mean wind optical depth at the sonic point. It is calculated by integrating the product of flux-weighted mean opacity and density inward from the outer boundary: $${\tau }_F(r_{\mathrm{s}\mathrm{o}\mathrm{n}\mathrm{i}\mathrm{c}})={\displaystyle {\int }_{r_{\mathrm{s}\mathrm{o}\mathrm{n}\mathrm{i}\mathrm{c}}}^{\mathrm{\infty }}{\aleph }_F(r)\rho (r)\mathrm{d}r.}$$(4)

Since we obtain a detailed solution to the radiative transfer problem, we predict the flux-weighted mean opacity and therefore output the wind optical depth.

Finally, converged models generate emergent synthetic spectra through formal integration in the observer’s frame, incorporating relevant broadening mechanisms. We additionally computed ionising photon fluxes (photons per second) by integrating the emergent spectrum shortward of relevant continuum edges, for instance, the Lyman continuum edge yields the H-ionising flux. These predicted quantities from our models are detailed below, while Appendix B provides a comprehensive table of wind-related outputs from our PoWR^HD models.

3 Radiative acceleration

A major advance in our understanding of hot-star winds over recent decades has been recognising the importance of proximity to the Eddington limit (Gräfener & Hamann 2008; Vink et al. 2011; Sander & Vink 2020). Yet predicting a priori how metal-licity and other stellar parameters affect VMS winds remains challenging, which is precisely why detailed radiative transfer calculations, such as those performed here are essential.

Nevertheless, we can infer qualitative trends to first order from the force balance in stellar winds (see Sabhahit et al. 2025, for a didactic study). The balance of relevant forces in massive star winds is characterised by the total radiative Eddington parameter, defined as the ratio of total radiative acceleration to gravity: $${\mathrm{\Gamma }}_{\mathrm{r}\mathrm{a}\mathrm{d}}(r)={\displaystyle \frac{a_{\mathrm{r}\mathrm{a}\mathrm{d}}(r)}{G{M}_{\star }/r^2}}={\displaystyle \frac{{\aleph }_F(r)L(r)}{4\pi Gc{M}_{\star }}}\approx {\displaystyle \frac{{\aleph }_F(r)L_{\star }}{4\pi Gc{M}_{\star }}},$$(5)

where the radial variation in total Γ_rad primarily reflects changes in the flux-weighted mean opacity. The luminosity L(r) technically varies with radius r, since a fraction of the photon energy is used to lift the material out of the stellar gravitational potential and to accelerate it to a certain kinetic energy. However, even for powerful WR-type winds, the approximation, L(r) ≈ L_*, is justified, as the deviations typically remain below ~5% (e.g. Sander et al. 2020).

Figure 1 illustrates the radial variation of total Γ_rad for a representative very massive star model with M_* = 105 M_⊙ and T_* = 35 kK. We decompose this into contributions from boundbound line transitions, Γ_lines, and continuum processes, Γ_cont, which includes free-free (ff) and bound-free (bf) continuum terms plus classical electron (Thomson) scattering Γ_e: $$\begin{array}{rl}{\mathrm{\Gamma }}_{\mathrm{r}\mathrm{a}\mathrm{d}}& ={\mathrm{\Gamma }}_{\mathrm{l}\mathrm{i}\mathrm{n}\mathrm{e}\mathrm{s}}+{\mathrm{\Gamma }}_{\mathrm{c}\mathrm{o}\mathrm{n}\mathrm{t}}\\ & ={\mathrm{\Gamma }}_{\mathrm{l}\mathrm{i}\mathrm{n}\mathrm{e}\mathrm{s}}+{\mathrm{\Gamma }}_{\mathrm{f}\mathrm{f}}+{\mathrm{\Gamma }}_{\mathrm{b}\mathrm{f}}+{\mathrm{\Gamma }}_{\mathrm{e}}.\end{array}$$(6)

The classical Γ_e provides a useful measure of the star’s proximity to the Eddington limit, maintaining near-constancy throughout the atmosphere (solid magenta in Fig. 1). The classical Γ_e is given by: $${\mathrm{\Gamma }}_{\mathrm{e}}={\displaystyle \frac{a_{\mathrm{t}\mathrm{h}\mathrm{o}\mathrm{m}}}{G{M}_{\star }/r^2}}\approx {\displaystyle \frac{{\sigma }_{\mathrm{e}}{L}_{\star }}{4\pi Gc{M}_{\star }}},$$(7)

where σ_e denotes the electron scattering opacity. For sufficiently high temperatures capable of fully ionising H and He (electron temperatures above approximately 10 and 50kK, respectively), the electron scattering opacity remains essentially constant. While this generally holds throughout the inner, sub-critical region of the atmosphere where such conditions are realised, the low-temperature sequences in our grid can show a slight decrease in Γ_e near the critical point. Therefore, the Γ_e values reported in the subsequent sections are averages taken within the sub-critical region.

Two important conclusions emerge from examining the continuum component. First, the true continuum contribution from ff+bf transitions remains vanishingly small relative to electron scattering, evident from the minimal difference between Γ_cont and Γ_e terms. Second, although the continuum component dominates total Γ_rad in the sub-critical inner atmosphere - accounting for over 70% of the total - it cannot launch the wind. Wind launching requires Γ_rad to exceed unity. As the figure demonstrates, continuum processes alone provide insufficient acceleration to overcome gravity; line transitions are essential for wind initiation.

The dominant contributors to line acceleration in the sub-critical region where mass-loss rates are established are line opacities of FeII-FeVII (depending on T_*). Immediately below the critical point, Fe provides the necessary acceleration to increase total Γ_rad to unity. Despite its relatively low abundance compared to other surface metals, the efficiency of Fe wind driving is due to the spectral coincidence between millions of Fe transitions and the ultraviolet peak characteristic of hot stellar photospheres.

Beyond the critical point, line opacities from several additional elements become important, extending beyond the CNO group. Notable contributors to outer-wind driving include Si, P, S , Cl, and Ar (see also Vink et al. 1999; Sander et al. 2017). Contributions from Si, P, and S are significant in our WNh models, as well as in the B hypergiant models applied to ζ¹ Sco by Bernini-Peron et al. (2025) and in late-WNh models from Lefever et al. (2025). This contrasts with pure-He classical WR winds, where these contributions are negligible (Sander et al. 2020), meaning the contribution of Si, P, and S is generally important for stars on or rightward (cooler) of the ZAMS.

Iron’s pivotal role in inner-wind driving is further demonstrated by its absolute and relative contributions to Γ_lines. Figure 2 presents a stacked plot comparing contributions from various components to total Γ_rad and Γ_lines at the critical point as a function of stellar mass for our T_* = 35 kK temperature sequence. In the top sub-panel, as stellar mass decreases, the Γ_e component at the critical point increases monotonically while line transition contributions decrease. The Fe contribution remains roughly constant at approximately 40% across most of the parameter space before dropping to roughly 20% at the lowest masses. In the lowest-mass model, despite Fe contributing only 20% of total Γ_rad, wind launching remains almost entirely Fe-driven. In the bottom sub-panel, Fe contributes more than half of the total line acceleration and remains approximately constant with mass. As we demonstrate in Sect. 4.3, insufficient Fe driving can lead to steep drops in the mass-loss rate in certain temperature and mass regimes.

Fig. 1 Total radiative acceleration and its contributions from continuum processes and line transitions. Decomposition of the total radiative acceleration (solid black) into continuum (dashed green) and line (dashed red) components along with individual contributions from constituent elements in the atmosphere (coloured lines with symbols, see legend). Each of these element contributions consists of the sum of their line and continuum accelerations. All accelerations are normalised to gravity. The dotted and dashed vertical grey lines mark the sonic and critical points, Rsonic and Rcrit, respectively. The model stellar parameters are log(L*/L⊙) = 6.0, M* = 105 M⊙, T* = 35kK, X = 0.7, and Z = 0.02.

Fig. 2 Contributions from electron scattering, free-free and bound-free continuum transitions, and line transitions of individual elements at the critical point, shown as a function of mass. The contributions are normalised to gravity and expressed as Eddington parameters. Top : absolute contributions from different opacity sources. Bottom : relative line contributions from different elements. The model sequence corresponds to a temperature of T* = 35 kK.

4 Winds of VMSs

4.1 Testing against the transition mass-loss rate

Determining mass-loss rates for hot, massive stars presents significant challenges, with several complementary approaches available (see Vink 2022, for a recent review). Empirical determinations typically involve fitting observed wind-sensitive spectral lines such as Hα or UV diagnostics using non-LTE expanding atmosphere models. In this iterative process, stellar and wind parameters are varied until optimal spectral fits are achieved, from which one can infer the best-fit empirical mass-loss rates. Alternatively, we can predict theoretical massloss rates using analytical models based on CAK theory and its subsequent improvements accounting for complex line-lists and finite-disc effects (Castor et al. 1975; Friend & Abbott 1986; Pauldrach et al. 1986), semi-analytical approaches such as MC simulations (e.g. Lucy & Abbott 1993; Vink et al. 2001; Müller & Vink 2008), or hydrodynamically consistent stellar windatmosphere models (Gräfener & Hamann 2005; Sander & Vink 2020), such as those investigated here.

Each approach comes with their own inherent limitations. Optical wind-line diagnostics - mainly recombination lines - are inherently model-dependent and strongly sensitive to assumptions about wind clumping in the underlying atmosphere models. As the adopted clumping factor increases, empirically derived rates systematically decrease (Hamann & Koesterke 1998). Furthermore, degeneracies between stellar and wind parameters can make it impossible to uniquely determine best-fit values. Theoretical predictions also suffer from clumping uncertainties, with predicted mass-loss rates depending on clumping assumptions near the critical point (Gräfener & Hamann 2005; Muijres et al. 2011; Sabhahit et al. 2025; Bernini-Peron et al. 2025).

To address these challenges, Vink & Gräfener (2012) developed the concept of the transition mass-loss rate. The underlying idea is straightforward: at the observed spectral morphological transition from absorption-dominated O star spectra to emissionline-dominated WNh star winds (as seen in the Hβ transition described by Crowther & Walborn 2011), the wind optical depth crosses approximately unity.

Building on the hydrodynamic equation of motion, Vink & Gräfener (2012) derived an analytical formula linking wind efficiency and optical depth to describe this spectral transition condition: $$\begin{array}{c}\eta \approx f\cdot {\tau }_F(r_{\mathrm{s}\mathrm{o}\mathrm{n}\mathrm{i}\mathrm{c}}),{\tau }_F(r_{\mathrm{s}\mathrm{o}\mathrm{n}\mathrm{i}\mathrm{c}})\approx 1.\end{array}$$(8)

Through extensive testing across a range of stellar parameters, Vink & Gräfener (2012) found that the factor f depends primarily on the terminal-to-escape velocity ratio υ_∞/υ_esc and on the velocity-law exponent β. For Galactic metallicity and υ_∞/υ_esc = 2.5, they obtained f ≈ 0.6, assuming a β-type velocity law with relatively large values of β ≈ 1.5, typical of very massive-star winds (Vink et al. 2011; Sabhahit et al. 2025). The dependence of f on the υ_∞/υ_esc ratio was parametrised in Sabhahit et al. (2023) as $f\sim 0.75\cdot (1+{\upsilon }_{\mathrm{e}\mathrm{s}\mathrm{c}}^2/{\upsilon }_{\mathrm{\infty }}^2{)}^{-1}$ (for an analytical derivation, see Gräfener et al. 2017). Although clumping assumptions can in principle influence the υ_∞/υ_esc ratio, and therefore the value of f, the impact is small because f approaches a plateau once the υ_∞/υ_esc ratio is sufficiently high.

This means the transition from optically thin O-star winds to optically thick WNh-star winds approximately coincides with both the wind efficiency parameter η and τ_F(r_sonic) crossing order unity. Beyond this transition, Vink et al. (2011) discovered a steep upturn in predicted mass loss with the classical Eddington parameter, creating the characteristic kink in the Ṁ(Γ_e) relation discussed in the introduction.

Vink & Gräfener (2012) rearranged Eq. (8) to formulate the concept of the transition mass-loss rate: $${\dot{M}}_{\mathrm{t}\mathrm{r}\mathrm{a}\mathrm{n}\mathrm{s}}=f\cdot {\displaystyle \frac{L_{\star }}{c{\upsilon }_{\mathrm{\infty }}}}.$$(9)

As Eq. (9) reveals, the transition mass-loss rate is essentially a luminosity estimate expressed in mass-loss units. The luminosity and terminal velocity of transition objects can be determined with greater accuracy than empirical measures of the mass-loss rate as the latter depend on model specifics such as clumping assumptions. Uncertainties in f typically fall within the error bars of luminosity measurements. This gives the transition massloss rate a crucial advantage: it is both accurate and significantly less model-dependent than traditional estimates. However, this accuracy applies specifically to the transition region between O and WNh stars, not beyond it.

Consider the Arches Cluster which hosts several O stars and WNh stars. The transition objects - whose spectral morphology bridges absorption-dominated O stars and emissionline-dominated WNh stars - have luminosities and terminal velocities of approximately log(L_*/L_⊙) ≈ 6.06(±0.1) and υ_∞ ≈ 2000(±400) km s⁻¹ (Martins et al. 2008). Applying Eq. (9) yields a transition mass-loss rate of log(Ṁ_trans) ≈ −5.16 (Vink & Gräfener 2012; Sabhahit et al. 2022).

Given the accurate and model-independent nature of this approach, any credible mass-loss prescription - whether empirical or theoretical - should reproduce this rate at the luminosity and temperature of the transition stars. While agreement at the transition does not guarantee accuracy elsewhere, it provides an essential first-order sanity check: if a prescription fails to match the transition mass-loss rate, we can confidently consider it inaccurate within this regime.

Figure 3 shows how our PoWR^HD models perform against the Arches transition mass-loss rate. We compare our predicted rates for the T_* = 35 kK model sequence with the transition massloss rate of log(Ṁ_trans) ≈ −5.16. The T_* = 35 kK sequence is selected such that its surface temperatures, defined at Rosseland optical depth τ_R = 2/3 as T_eff(τ_R = 2/3), is close to the observed T_eff ≈ 33.4kK of the transition objects. Our PoWR^HD models successfully reproduce the transition mass-loss rate at stellar parameters corresponding to the Arches transition objects.

We also compare the Arches transition mass-loss rate with the widely used Vink et al. (2000) mass-loss predictions. As shown in Figure 3 - and as already noted by Vink & Gräfener (2012) - the Vink et al. (2000) rates agree well with the transition mass-loss rates for the relevant stellar parameters. However, the Vink et al. (2000) rates diverge from our predictions both above and below the transition point. This deviation becomes particularly pronounced at lower masses or higher Γ_e, where our predicted mass loss exhibits a steeper scaling with mass and consequently with Γ_e above the transition. The presence of a kink in the mass-loss rate agrees with the results of Vink et al. (2011).

Our models predict a terminal velocity of 1500 km s⁻¹ at the transition, approximately 500 km s⁻¹ lower than observed in Arches transition objects. However, terminal velocity predictions are notoriously sensitive to outer wind conditions, including temperature stratification and clumping, both of which remain uncertain. We can achieve better agreement by increasing outerwind clumping, which enhances outer-wind driving and boosts terminal velocity. Indeed, such an increase in outer-wind clumping was required to reproduce the UV+optical spectra of R136a1 and R144 using hydrodynamical models (Sabhahit et al. 2025).

Figure 3 also compares our predicted wind efficiency parameter and optical depth with transition values from Vink & Gräfener (2012). Additionally, Figs. 4 and 5 display synthetic optical and UV spectra for the T_* = 35 kK model sequence, highlighting crucial wind-sensitive recombination and P-Cygni lines. At the point where our models match the Arches transition mass-loss rate, we find η ≈ 0.45. This value falls below the Vink & Gräfener (2012) transition value of η ≈ 0.6, likely reflecting our underestimated terminal velocities. The wind optical depth crosses unity precisely at the transition, coinciding with the spectral morphological change where Hβ shifts from absorption to emission. We observe similar behaviour in Hγ, while Hα transitions from absorption to emission at higher masses (i.e, lower Γ_e) compared to other Balmer lines. In contrast, He II λ4686 remains in emission throughout the sequence, strengthening as mass decreases. In the UV, one can see the expected sensitivity of the the CIV λ1550 and SiIV λ1400 resonance doublets to the terminal velocity in our model. Additional features also appear in the UV, for example around λ1470, which is an iron feature that strengthens as the winds become optically thick.

Finally, we can estimate the stellar mass required in our hydrodynamic models to match the transition mass-loss rate. We intend to closely match the (averaged) observed stellar luminosity and surface temperature of the transition object in the Arches. To this end, we interpolate between our model sequences to precisely match the observed T_eff(τ_R = 2/3) ≈ 33.4 kK (Martins et al. 2008). We find that M_* ≈ 60 M_⊙ (corresponding to Γ_e ≈ 0.43) and T_* ≈ 38.2 kK is required to match the transition. This M_* ≈ 60 M_⊙ coincides with the kink location in our models and agrees well with predictions from Sabhahit et al. (2022) for the location of the kink in Γ_e-space.

To summarise our findings: our PoWR^HD models successfully reproduce the transition mass-loss rate in the Arches Cluster, confirm the prediction of a mass-loss kink above the transition, and capture the morphological change in Hβ from absorption to emission as the wind optical depth crosses unity. While our model-predicted wind efficiency at the transition (η ≈ 0.45) is lower than Vink & Gräfener (2012) predicted η ≈ 0.6, this discrepancy likely stems from our underestimated terminal velocities. Overall, these results provide strong validation of our hydro-predicted mass-loss rate and its ability to capture the essential physics of the O to WNh star morphology transition.

Fig. 3 Testing the PoWRHD model predictions against the transition mass-loss rate in the Arches Cluster. Left : predicted mass-loss rates versus stellar mass from PoWRHD models for the T* = 35 kK sequence (solid black) and from Vink et al. (2000) (dotted blue), compared to the transition mass-loss rate of log(M) ≈ −5.16 in the Arches (horizontal red line) from Vink & Gräfener (2012). The PoWRHD models are computed with fixed log(L*/L⊙) = 6.0, X = 0.7, and Z = 0.02, roughly matching the observed properties of the Arches transition objects (Martins et al. 2008). The T* = 35 kK sequence shown has surface temperatures close to the observed Teff(τR = 2/3) ≈ 33.4 kK of the transition objects. The stellar masses of the transition objects span ~55—75 M⊙ corresponding to their averaged luminosity and temperature. Middle : predicted wind efficiency parameter versus mass compared to the transition value η ≈ 0.6 (dashed red) from Vink & Gräfener (2012). Right : predicted wind optical depth versus mass compared to the transition value τF(rsonic) ≈ 1 (dashed red) from Vink & Gräfener (2012).

Fig. 4 Synthetic spectra showing relevant wind lines in the optical for the T* = 35 kK model sequence presented in Fig. 3. The spectra are colour-coded based on the stellar mass.

Fig. 5 Synthetic spectra showing relevant wind lines in the UV for the T* = 35 kK model sequence presented in Fig. 3. The spectra are colour-coded based on the stellar mass.

4.2 Wind properties: Ṁ, υ_∞,η, and τ_F(r_sonic)

Figure 6 presents relevant wind parameters, including massloss rates and terminal wind velocities, as functions of stellar mass, with temperature sequences colour-coded by T_*. Across all temperature sequences, we observe a consistent fundamental trend: mass-loss rates increase systematically towards lower stellar masses when other parameters remain fixed. This behaviour reflects the underlying physics: For fixed luminosity and composition, decreasing the stellar mass decreases the gravitational force in the wind and drives up the total Γ_rad in the sub-critical region (r < R_crit), directly boosting mass-loss rates (cf. Castor et al. 1975; Lamers & Cassinelli 1999; Vink et al. 2001; Sander & Vink 2020; Sabhahit et al. 2025).

Terminal velocities follow a complementary pattern, generally decreasing with decreasing mass (cf. Castor et al. 1975; Friend & Abbott 1986), although the behaviour can be nonmonotonic sometimes (see e.g. the sequence in Lefever et al. 2025; Bernini-Peron et al. 2025). To first order, an increase in the total Γ_rad in the sub-critical regime (r < R_crit) drives a higher mass loss but produces denser, slower winds. Conversely, increased Γ_rad in the supercritical region (r > R_crit) accelerates winds to higher velocities without affecting mass-loss rates (see also Vink et al. 2001; Sander & Vink 2020).

More intriguing are the temperature-dependent variations in these scaling relations. Hotter models with T_* ≳ 30 kK reveal a distinctive upturn - a characteristic kink - in the log Ṁ-log M (and therefore also the log Ṁ - log Γ_e) relation. The scaling behaviour transitions from the shallow slope of ~ 1.5-3 at higher masses (lower Γ_e) to a dramatically steeper slope of ~10 at lower masses (higher Γ_e). This behaviour qualitatively matches the MC predictions of Vink et al. (2011) and aligns with empirical evidence of such a kink in the 30 Doradus cluster reported by Bestenlehner et al. (2014).

In contrast, for the cooler models with T_* ≲ 25 kK, their log Ṁ - log M relations maintain consistently steep slopes throughout our parameter range, showing no pronounced kink behaviour. Particularly notable are the T_* = 25 and 15 kK sequences, which maintain distinctly steep slopes across the entire mass range.

Terminal velocity behaviour further distinguishes hot and cool regimes. Hotter models show terminal velocities that generally decline with decreasing mass, while cooler models exhibit weak mass dependence, remaining nearly constant at ~500 km s⁻¹. Similarly low velocities were also obtained in the late WNh models of Lefever et al. (2025).

Both wind efficiency parameters and optical depths increase monotonically with decreasing mass across all temperatures. At kink locations in hotter models, wind efficiency reaches approximately 0.4-0.45 - notably lower than the Vink & Gräfener (2012) transition value of η ≈ 0.6, likely reflecting our underestimated terminal velocities. On the other hand, wind optical depths cross unity precisely at these kink locations, supporting our physical interpretation.

4.3 Role of iron and the bistability jump

The steep mass-loss scaling at T_* = 25 and 15 kK demands a deeper investigation into the physical driving mechanisms. Figure 7 demonstrates the fundamental importance of Fe-line acceleration at the critical point. The upper sub-panel reveals a clear trend: absolute Fe acceleration decreases with mass, consistent with the increasing dominance of Γ_e over Γ_lines towards lower masses (see Sect. 3).

Temperature sequences at T_* = 25 and 15 kK, however, stand out with Fe-line contributions declining even at higher masses. While subtle indications of high-mass decline appear at various other temperatures, the effect becomes particularly pronounced at T_* = 25 and 15kK. This behaviour is even more striking in the relative Fe-line contributions to the total Γ_lines (lower sub-panel). At these specific temperatures, the relative Fe contribution drops toward higher masses, contrasting with the approximately constant or weakly declining trends observed at other temperatures. Although Fe continues to dominate sub-critical line acceleration across all models presented, an abrupt reduction occurs precisely at the critical point of outflow initiation for the T_* = 25 and 15 kK sequences. This sudden Fe-driving deficit produces the steep mass-loss rate decline at these characteristic temperatures. This effect was also seen in some of the late-type WNh sequences in Lefever et al. (2025).

The correspondence between the lack of Fe-line driving and reduced mass-loss rates at temperatures consistent with Vink et al. (1999) - where OB-star wind bistability occurs - is not coincidental. The top sub-panel in Figure 8 presents massloss rates as functions of T_* for various stellar masses. Stars close to their Eddington limit (corresponding to low mass for a fixed luminosity) exhibit straightforward monotonic behaviour: mass-loss rates increase with decreasing T_*, reflecting pure geometrical effect from the increase in radius, that is, for a fixed luminosity, lower temperatures mean larger radii and shallower gravitational potential wells (cf. Sander et al. 2023; Lefever et al. 2025). The low-temperature, high Γ_e regime proves particularly relevant for massive LBVs, which exhibit strong winds and, at times, eruptive mass-loss episodes. These models approaching their Eddington limit naturally predict strong winds in the temperature range characteristic of LBVs, consistent with empirical mass-loss rates determined during their quiescent phases.

Models with higher masses, that is, models with lower Γ_e, conversely, display complex non-monotonic temperature scaling with sudden jumps in the mass loss below T_* ≈ 25 and 15kK. The bottom sub-panel in Figure 8 illustrate the contributions of FeV, FeIV, FeIII, and FeII to the critical-point line driving. The ionisation sequence evolves predictably with temperature: Fe IV dominates wind driving at T_* ≈ 30 kK. As temperatures decreases, Fe IV recombines into Fe III which becomes the primary contributor near T_* ≈ 20 kK. At intermediate temperatures around T_* ≈ 25kK, neither FeIV nor FeIII dominates. Below T_* ≈ 25 kK, Fe III begins to dominate the wind driving and we see a rapid increase in the mass-loss rate.

These results show excellent qualitative agreement with the first bistability jump predicted for B-star winds using MC wind models, both with global and local consistency (Vink et al. 1999; Vink 2018; Vink & Sander 2021), as well as other wind codes such as CMFGEN and METUJE (Petrov et al. 2016; Krticka et al. 2021). Empirical evidence appears to exist for the presence of the bistability jump in individual LBV objects (e.g. Groh et al. 2009); but such evidence does not appear to be present in empirical modelling results of O and B supergiant groups (Crowther et al. 2006).

The bistability pattern repeats at lower temperatures: Fe III recombines below T_* ≈ 20 kK, creating another driving minimum near T_* ≈ 15 kK where neither Fe III nor Fe II dominates. This produces our second mass-loss rate minimum. Finally, Fe II takes over below 15 kK, giving rise to yet another bistability jump. Crucially, the Fe-driving minima at T_* ≈ 25 and 15 kK coincide with our models maintaining a steep log Ṁ-log M scaling across the entire investigated mass range (see top sub-panel in Fig. 6).

These results demonstrate how ionisation physics governs wind behaviour: the transitions between dominant Fe ionisation states create temperature windows where Fe driving becomes ineffective, which ultimately gives rise to a steep decline in the mass-loss rate in this regime.

Fig. 6 Mass-loss rates, terminal wind velocities, wind efficiencies, and wind optical depths predicted from hydrodynamically consistent PoWRHD models as a function of stellar mass M*. The plots are colour-coded according to the different input inner boundary temperatures, ranging from T* = 50 to 12kK.

Fig. 7 Contribution of Fe-line acceleration at the critical point as a function of mass for different temperature sequences. Top : absolute contribution of Fe. Bottom : relative contribution of Fe to the total line acceleration. The contributions are normalised to gravity and expressed as Eddington parameters. The figure is colour-coded based on the temperature T*.

Fig. 8 Mass-loss rate as a function of stellar temperature, T*, for different masses. Top : predicted mass-loss rates. Bottom : relative contributions from relevant Fe wind-driving ions to the total line acceleration at the critical point. The radiative acceleration contributions are normalised to gravity and expressed as Eddington parameters. The curves are colour-coded by stellar mass.

4.4 Mass-loss scaling with mass and temperature

In this section, we present a fitting relation to predict the massloss rate as a function of M_* and T_*. Since our investigation varies only these parameters, the following formula does not represent a complete mass-loss prescription. Rather, it is designed to capture the scaling relations and functional form, which will serve as a foundation for future work, as we plan to extend this relation by incorporating other relevant parameters in our grid such as luminosity and surface abundances.

Our fit combines a double-exponential mass scaling - capturing both shallow and steep scalings - with quadratic temperature dependence and Gaussian dips reproducing the bistability jumps. These Gaussian terms capture the mass-loss rate minima near T_* ≈ 25kK and l5kK³: $$\begin{array}{rl}\log \dot{M}(M_{\star },T_{\star })& =a+\log [{10}^{f_{\mathrm{l}\mathrm{o}\mathrm{w}}(M_{\star })}+{10}^{f_{\mathrm{h}\mathrm{i}\mathrm{g}\mathrm{h}}(M_{\star })}]\\ & +b(M_{\star })\log \left(T_{\star }/T_{\mathrm{r}\mathrm{e}\mathrm{f}}\right)+c[\log \left(T_{\star }/T_{\mathrm{r}\mathrm{e}\mathrm{f}}\right){]}^2\\ & -G_1(M_{\star })\exp [-{\left({\displaystyle \frac{T_{\star }-T_1}{{\sigma }_T}}\right)}^2]\\ & -G_2(M_{\star })\exp [-{\left({\displaystyle \frac{T_{\star }-T_2(M_{\star })}{{\sigma }_T}}\right)}^2],\end{array}$$(10)

where M_* is in solar masses and T_* is in kK. The bistability temperatures T₁ and T₂(M_*) mark the centres of the Gaussian dips, while σ_T controls their width. We fix T₁ = 25 kK and σ_T = 2 kK during fitting, as our models demonstrate negligible variation of these parameters with mass or temperature. Stellar mass and temperature are normalised to references values, M_ref and T_ref, with M_ref = 60 M_⊙ at T_ref = 38 kK. These normalisation values are not arbitrary - as demonstrated in Sect. 4.1, they roughly correspond to the M_* and T_* values required by our hydrodynamic models to reproduce the Arches Cluster transition mass-loss rate. The mass-dependent coefficients take the following form: $$\begin{array}{rl}a& =-5.527(\pm 0.022)\\ f_{\mathrm{l}\mathrm{o}\mathrm{w}}(M_{\star })& =-1.864(\pm 0.198)\cdot \log (M_{\star }/M_{\mathrm{r}\mathrm{e}\mathrm{f}}),\\ f_{\mathrm{h}\mathrm{i}\mathrm{g}\mathrm{h}}(M_{\star })& =-9.865(\pm 0.273)\cdot \log (M_{\star }/M_{\mathrm{r}\mathrm{e}\mathrm{f}})\\ b(M_{\star })& =-2.062(\pm 0.521)+5.671(\pm 1.536)\cdot \log (M_{\star }/M_{\mathrm{r}\mathrm{e}\mathrm{f}})\\ c& =3.974(\pm 0.459)\\ G_1(M_{\star })& =0.178(\pm 0.042)\exp [9.611(\pm 1.293)\cdot \log (M_{\star }/M_{\mathrm{r}\mathrm{e}\mathrm{f}})]\\ G_2(M_{\star })& =0.291(\pm 0.051)\exp [8.658(\pm 1.048)\cdot \log (M_{\star }/M_{\mathrm{r}\mathrm{e}\mathrm{f}})]\\ T_2(M_{\star })& =16.941(\pm 0.268)-7.274(\pm 1.506)\cdot \log (M_{\star }/M_{\mathrm{r}\mathrm{e}\mathrm{f}}),\\ M_{\mathrm{r}\mathrm{e}\mathrm{f}}& =60-0.521(\pm 0.133)\cdot (T_{\star }-T_{\mathrm{r}\mathrm{e}\mathrm{f}})\\ T_1& =25\mathrm{k}\mathrm{K},{\sigma }_T=2\mathrm{k}\mathrm{K},T_{\mathrm{r}\mathrm{e}\mathrm{f}}=38\mathrm{k}\mathrm{K}.\end{array}$$

Several simplifications help elucidate the different scaling relations in Eq. (10). Our fitting function ensures that at high temperatures (T_* ≳ 30 kK) - beyond the influence of bistability features - the Gaussian corrections vanish, leaving only the baseline mass and temperature scaling. The double-exponential term 10^f_high(M_*) + 10^f_low(M_*) captures a critical feature: the transition from shallow to steep power-law scaling of log Ṁ versus log M_* at higher temperatures. For the simplified case of T_* = 38kK, Eq. (10) reduces to $$\log \dot{M}(M_{\star },38\mathrm{k}\mathrm{K})\approx a+\log [{10}^{f_{\mathrm{l}\mathrm{o}\mathrm{w}}(M_{\star })}+{10}^{f_{\mathrm{h}\mathrm{i}\mathrm{g}\mathrm{h}}(M_{\star })}],$$(11)

which simplifies further depending on which term dominates based on the stellar mass relative to 60 M_⊙. For high-mass stars with M_* > 60 M_⊙ or low Γ_e, the first term dominates, while for lower-mass stars, M_* < 60 M_⊙, with a high Γ_e, the second term takes over: $$\begin{array}{r}\log \dot{M}(M_{\star },38\mathrm{k}\mathrm{K})\approx \{\begin{array}{ll}a-1.864\cdot \log (M_{\star }/60)& \text{if\hspace{0.17em}}{M}_{\star }>60\\ & \text{or low\hspace{0.17em}}{\mathrm{\Gamma }}_{\mathrm{e}}\\ a-9.865\cdot \log (M_{\star }/60)& \text{if\hspace{0.17em}}{M}_{\star }<60\\ & \text{or high\hspace{0.17em}}{\mathrm{\Gamma }}_{\mathrm{e}}.\end{array}\end{array}$$(12)

The double-exponential term therefore accurately captures the shallow-to-steep transition in the log Ṁ versus log M_* scaling predicted by our models. Figure 9 shows our fit relations, where for T_* ≳ 30 kK, the log Ṁ versus log M_* scaling steepens below M_* ≈ 60 M_⊙.

To understand how the baseline mass-loss rate scales with temperature, we fix M_* = 60 M_⊙. Neglecting the minor variation of M_ref with T_* and the double-Gaussian bistability dips, the fit relation reduces to $$\begin{array}{rl}\log \dot{M}(60M_{\odot },T_{\star })& \approx a+\mathrm{l}\mathrm{o}\mathrm{g}2-2.062\cdot \mathrm{l}\mathrm{o}\mathrm{g}(T_{\star }/38)\\ & +3.974\cdot [\log \left(T_{\star }/38\right){]}^2.\end{array}$$(13)

In Fig. 10, we show our fitted relation as a function of T_* from Eq. (10), compared to our model predictions. The baseline scaling of log Ṁ(T_*) with log T_* exhibits a steep negative dependence. For example, at T_* = 38 kK, differentiating Eq. (13) gives an instantaneous slope of $\dot{M}\propto T_{\star }^{-2.062}$. As T_* decreases, this slope further steepens. We also compare our predictions with reference lines showing a $\dot{M}\propto T_{\star }^{-6}$ scaling, which arise from the geometrical effect of increasing radii and decreasing surface gravity. We find that in certain parameter regimes, the mass-loss scaling with temperature is dominated by this purely geometrical effect, particularly at higher Γ_e and lower T_*.

We acknowledge a practical limitation of Eq. (10): the temperature T_* is defined at an arbitrary Rosseland continuum optical depth of τ_R,cont = 20, which may not be readily accessible in stellar structure and evolution models (see e.g. the recent discussion in Josiek et al. 2025). Most structure codes compute the surface temperature T_eff at the photospheric surface where τ_R = 2/3, assuming a static, grey atmosphere. Furthermore, due to the arbitrary choice of τ_R,cont, defining the effective temperature at τ_R = 2/3 provides a more physically meaningful measure of stellar temperature.

We therefore provide a supplementary relation connecting mass-loss rate with both T_* and T_eff(τ_R = 2/3) for our atmosphere models, enabling temperature parameter mapping based on the predicted Ṁ. Figure 11 illustrates the relationship between T_* and T_eff(τ_R = 2/3). The difference between the two temperatures increases with increasing mass-loss rate (cf. Smith et al. 2004; Sander et al. 2023), as the wind becomes more optically thick. We adopt a linear relation between T_* and T_eff (τ_R = 2/3), with both slope and y-intercept scaling with mass-loss rate: $$T_{\star }=a(\dot{M})+b(\dot{M})\cdot T_{\mathrm{e}\mathrm{f}\mathrm{f}}({\tau }_{\mathrm{R}}=2/3),$$(14)

where $$\begin{array}{rl}a(\dot{M})& =1-0.85(\pm 0.277)({\displaystyle \frac{\dot{M}}{{\dot{M}}_{\mathrm{r}\mathrm{e}\mathrm{f}}}}{)}^{0.901(\pm 0.097)}\\ b(\dot{M})& =1+0.152(\pm 0.14)({\displaystyle \frac{\dot{M}}{{\dot{M}}_{\mathrm{r}\mathrm{e}\mathrm{f}}}}{)}^{0.925(\pm 0.319)}\\ {\dot{M}}_{\mathrm{r}\mathrm{e}\mathrm{f}}& =6.92\cdot {10}^{-6}{M}_{\odot }\mathrm{y}{\mathrm{r}}^{-1}.\end{array}$$(15)

Here, Ṁ is in M_⊙ yr⁻¹, while temperatures are in kK. The coefficients a(Ṁ) and b(Ṁ) are constructed such that at very low mass loss, the difference saturates to ∆T = T_* - T_eff (τ_R = 2/3) ~ 1 kK in agreement with our model predictions in Fig. 11. While the typical error from our T_* fits is of the order of 1 kK, especially for the highest mass-loss rate models, it successfully captures the overall trend of increasing difference between the temperatures towards higher mass-loss rate. Therefore, given M_* and T_eff(τ_R = 2/3), an iterative solution of Eqs. (10) and (14) simultaneously yields Ṁ and T_*.

While we present fits capturing the temperature mapping between T_* and T_eff(τ_R = 2/3) for the current set of atmosphere models, we emphasise that these relations are preliminary and should not be over-interpreted. For example, the T_eff(τ_R = 2/3) values calculated in evolution models assume a grey static atmosphere connected to the outer boundary. Consequently, the T_eff(τ_R = 2/3) from evolution models could be inaccurate, due to effects such as inflation and turbulence associated with the hot iron bump (see e.g. Josiek et al. 2025, for potential issues with connecting structure and atmosphere models). Additionally, since we varied only M_* in this work, further testing with variations in other stellar parameters is necessary to establish whether this correction is universally applicable.

Finally, we perform a sanity check of our fitting relations using the transition mass-loss rate obtained from applying Eq. (9) to the transition objects in the Arches Cluster. Section 4.1 demonstrated that the mass-loss kink occurs at M_* ≈ 60 M_⊙. With mass dependencies normalised to 60 M_⊙ and T_eff(τ_R = 2/3) = 33.4kK, solving Eqs. (10) and (14) gives log Ṁ ≈ −5.22 (and T_* = 38.05 kK). This value closely matches the Arches transition mass-loss rate, validating our fitting relations.

Fig. 9 PoWRHD predicted mass-loss rates (symbols with dashed lines) compared to our best-fit relation (solid lines) as functions of stellar mass, M*. The fit derives from Eq. (10), with lines colour-coded by stellar temperature, T*.

Fig. 10 Same as Fig. 9, but presenting mass-loss rates as functions of stellar temperature, T*, colour-coded by stellar mass, M*. The dashed black lines indicate reference scalings of $\dot{M}\propto T_{\star }^{-6}$ scaling, which arise from geometric effects.

Fig. 11 Change in the inner boundary temperature T* as a function of the effective temperature defined at Rosseland optical depth τR = 2/3 as winds become optically thick. Symbols show our PoWRHD model predictions, colour-coded by mass-loss rate. Contours of constant massloss rate are shown with dashed lines, while solid lines denote the fits from Eq. (14). The dotted black line marks where the two temperatures are equal.

5 Discussion: Comparison to previous works

We now place our model predictions in context with previous VMS wind investigations. Gräfener & Hamann (2008) examined VMS winds across a comparable parameter space using a precursor to our PoWR^HD version. Our models improve upon theirs through expanded atomic data and broader coverage of the Γ_e - T_* parameter space (see Appendix B). While Gräfener & Hamann (2008) identified steep scaling of mass loss with Γ_e, their predicted mass-loss rates lacked kink behaviour - possibly because their models did not cover the Γ_e < 0.4 regime where slope shallowing occurs in our calculations. We further check the absolute mass-loss rate predictions for M_* = 60 M_⊙ and T_* = 38 kK where we identified our kink to match the transition mass-loss rate. There is roughly 0.2 dex difference between our predictions, with our rates higher compared to Gräfener & Hamann (2008) (see Fig. A.1). This difference could be due to our expanded atomic data, as well as our inclusion of turbulent pressure in the atmosphere.

For high-Γ_e VMSs, our models predict a steep mass-loss scaling with mass and consequently Γ_e of Ṁ ~ Γ_e¹⁰ (see Fig. 9). Although further work must establish whether such scaling persists in the Ṁ - Γ_e relation when varying luminosity and H content simultaneously, this exceeds the −5 scaling reported by Gräfener & Hamann (2008) and Vink et al. (2011).

The overlapping temperature range (30-50 kK) between our work and Gräfener & Hamann (2008) also reveals qualitative agreement between both studies: the baseline mass-loss rate increases with decreasing T_*.

VMS winds were further studied using MC techniques by Vink et al. (2011), who predicted the kink behaviour. While qualitative agreement exists between our predictions, the location of the kink in our models occurs at Γ_e close to 0.4, lower compared to Vink et al. (2011). However, this is consistent with Vink et al. (2011) expectations that their kink would shift toward lower values.

6 Conclusions

In this paper, we have investigated VMS wind properties using hydrodynamically consistent PoWR^HD stellar atmosphere models. We examined how stellar mass and temperature govern predicted mass-loss rates, terminal velocities, and wind characteristics across the parameter space relevant to WNh stars. Our principal findings are s follows:

1. Our models predict a characteristic upturn in the log Ṁ-log M relation for sequences with T_* ≳ 30 kK, producing a kink feature consistent with Monte Carlo simulations of Vink et al. (2011) and empirical evidence from Bestenlehner et al. (2014). The location of the kink in the MC models roughly coincides with the crossing of the singlescattering limit and the wind optical depth reaching order unity. With multiple scattering in the wind, the star can drive a stronger mass-loss rate, resulting in the appearance of a kink. Our CMF-based models corroborate the MC results; for instance, our models also predict the characteristic kink feature where the wind optical depth crosses unity, and the synthetic spectra transition from absorption-dominated O-type to emission-dominated H-rich WR morphologies. The wind efficiency parameter reaches η ~ 0.45 at the kink - slightly below the η - 0.6 predicted by Vink & Gräfener (2012). This discrepancy likely reflects our underprediction of terminal velocities compared to those inferred for transition-type stars in the Arches Cluster. For the chosen luminosity, the kink occurs at M_* ≈ 60 M_⊙ corresponding to Γ_e ≈ 0.43;

2. We tested our predictions against the model-independent transition mass-loss rate in the Arches Cluster near the Galactic center, offering a crucial first-order sanity check for any prescribed mass loss at the O to WNh boundary. Our mass-loss rates match the Arches transition value of log(Ṁ_trans) ≈ −5.16 (Vink & Gräfener 2012) at our predicted kink location. Importantly, we also reproduce the associated spectroscopic transition from absorption- to emission-dominated Hβ profiles as the wind optical depth passes unity, thereby confirming the physical interpretation of the transition mass-loss concept and its connection to optically thick wind onset in VMSs;

3. In addition to the kink, our models reveal dramatic temperature-dependent jumps in the predicted mass-loss rates, most prominently at T_* ≈ 25 and 15kK. These features arise from shifts in the dominant Fe ionisation stage - FeIV → FeIII near 25kK and FeIII → FeII near 15kK -creating radiative driving deficits where both upper and lower ionisation stages dominate. The resulting bistability jumps agree qualitatively with Monte Carlo line-driven wind predictions from both globally and locally consistent models (Vink et al. 1999; Vink 2018; Vink & Sander 2021), as well as CMFGEN and METUJE computations (Petrov et al. 2016; Krticka et al. 2021). At these specific temperatures, our predicted log Ṁ-log M relation maintains steep scaling throughout the investigated mass range;

4. A radiative acceleration analysis reveals the decisive role of Fe-group elements, which - despite low abundances -contribute over 50% of total line driving through myriad UV transitions. In the sub-critical wind regime where Ṁ is set, continuum processes alone are insufficient to overcome gravity and Fe line transitions initiate the wind;

5. To facilitate further applications, we provide fitting prescriptions (Equations (10) and (14)) that reproduce both the kink and the bistability behaviour. These relations combine a double-exponential mass scaling with Gaussian terms capturing the ionisation-driven minima, thereby encoding the complex temperature and mass dependencies. Although derived for fixed luminosity and composition, these formulae establish a quantitative foundation for future incorporation of additional parameters including luminosity, surface H fraction, and metallicity.

Together, these results represent significant progress toward predictive VMS mass-loss prescriptions. By capturing both the Eddington-parameter associated kink behaviour and Feionisation shift driven bistable winds, our hydrodynamically consistent models address fundamental gaps in VMS wind theory, with implications for stellar evolution, chemical feedback, and stellar fate.

Acknowledgements

We thank the anonymous referee for constructive comments that helped improve the paper. GNS and JSV are supported by STFC funding under grant number ST/Y001338/1. AACS is supported by the German Deutsche Forschungsgemeinschaft, DFG in the form of an Emmy Noether Research Group - Project-ID 445674056 (SA4064/1-1, PI: Sander). AACS further acknowledges financial support by the Federal Ministry for Economic Affairs and Climate Action (BMWK) via the German Aerospace Center (Deutsches Zentrum für Luft- und Raumfahrt, DLR) grant 50 OR 2503 (PI: Sander). This project was co-funded by the European Union (Project 101183150 - OCEANS).

References

Appendix A Additional plots

Fig. A.1 Mass-loss rate versus the classical Eddington parameter Γe, comparing our predictions with those of Gräfener & Hamann (2008). Our model predictions are indicated by coloured symbols connected with dashed lines. The best-fit relation from Eq (10) is shown as solid lines, while the fit from Gräfener & Hamann (2008) is shown as dotted lines. Absolute mass-loss rates predicted by the two fits for M* = 60 M⊙ and T* = 38 kK are indicated by red and black star symbols, respectively.

Appendix B Parameter space and predicted wind properties

In our PoWR^HD models, we use an extended set of 17 elements with their respective mass fractions provided in Table B.1. The total metal mass fraction adds up to Z = 0.02 with individual metal mass fractions distributed following solar-scaled abundances from Grevesse & Sauval (1998). The detailed breakdown of the relevant ion list and their respective level and line numbers are provided in Table B.2.

All Tables

All Figures

Fig. 1 Total radiative acceleration and its contributions from continuum processes and line transitions. Decomposition of the total radiative acceleration (solid black) into continuum (dashed green) and line (dashed red) components along with individual contributions from constituent elements in the atmosphere (coloured lines with symbols, see legend). Each of these element contributions consists of the sum of their line and continuum accelerations. All accelerations are normalised to gravity. The dotted and dashed vertical grey lines mark the sonic and critical points, Rsonic and Rcrit, respectively. The model stellar parameters are log(L*/L⊙) = 6.0, M* = 105 M⊙, T* = 35kK, X = 0.7, and Z = 0.02.

In the text

Fig. 2 Contributions from electron scattering, free-free and bound-free continuum transitions, and line transitions of individual elements at the critical point, shown as a function of mass. The contributions are normalised to gravity and expressed as Eddington parameters. Top : absolute contributions from different opacity sources. Bottom : relative line contributions from different elements. The model sequence corresponds to a temperature of T* = 35 kK.

In the text

Fig. 3 Testing the PoWRHD model predictions against the transition mass-loss rate in the Arches Cluster. Left : predicted mass-loss rates versus stellar mass from PoWRHD models for the T* = 35 kK sequence (solid black) and from Vink et al. (2000) (dotted blue), compared to the transition mass-loss rate of log(M) ≈ −5.16 in the Arches (horizontal red line) from Vink & Gräfener (2012). The PoWRHD models are computed with fixed log(L*/L⊙) = 6.0, X = 0.7, and Z = 0.02, roughly matching the observed properties of the Arches transition objects (Martins et al. 2008). The T* = 35 kK sequence shown has surface temperatures close to the observed Teff(τR = 2/3) ≈ 33.4 kK of the transition objects. The stellar masses of the transition objects span ~55—75 M⊙ corresponding to their averaged luminosity and temperature. Middle : predicted wind efficiency parameter versus mass compared to the transition value η ≈ 0.6 (dashed red) from Vink & Gräfener (2012). Right : predicted wind optical depth versus mass compared to the transition value τF(rsonic) ≈ 1 (dashed red) from Vink & Gräfener (2012).

In the text

	Fig. 4 Synthetic spectra showing relevant wind lines in the optical for the T* = 35 kK model sequence presented in Fig. 3. The spectra are colour-coded based on the stellar mass.
In the text

	Fig. 5 Synthetic spectra showing relevant wind lines in the UV for the T* = 35 kK model sequence presented in Fig. 3. The spectra are colour-coded based on the stellar mass.
In the text

	Fig. 6 Mass-loss rates, terminal wind velocities, wind efficiencies, and wind optical depths predicted from hydrodynamically consistent PoWRHD models as a function of stellar mass M*. The plots are colour-coded according to the different input inner boundary temperatures, ranging from T* = 50 to 12kK.
In the text

	Fig. 7 Contribution of Fe-line acceleration at the critical point as a function of mass for different temperature sequences. Top : absolute contribution of Fe. Bottom : relative contribution of Fe to the total line acceleration. The contributions are normalised to gravity and expressed as Eddington parameters. The figure is colour-coded based on the temperature T*.
In the text

	Fig. 8 Mass-loss rate as a function of stellar temperature, T*, for different masses. Top : predicted mass-loss rates. Bottom : relative contributions from relevant Fe wind-driving ions to the total line acceleration at the critical point. The radiative acceleration contributions are normalised to gravity and expressed as Eddington parameters. The curves are colour-coded by stellar mass.
In the text

	Fig. 9 PoWRHD predicted mass-loss rates (symbols with dashed lines) compared to our best-fit relation (solid lines) as functions of stellar mass, M*. The fit derives from Eq. (10), with lines colour-coded by stellar temperature, T*.
In the text

	Fig. 10 Same as Fig. 9, but presenting mass-loss rates as functions of stellar temperature, T*, colour-coded by stellar mass, M*. The dashed black lines indicate reference scalings of $\dot{M}\propto T_{\star }^{-6}$ scaling, which arise from geometric effects.
In the text

Fig. 11 Change in the inner boundary temperature T* as a function of the effective temperature defined at Rosseland optical depth τR = 2/3 as winds become optically thick. Symbols show our PoWRHD model predictions, colour-coded by mass-loss rate. Contours of constant massloss rate are shown with dashed lines, while solid lines denote the fits from Eq. (14). The dotted black line marks where the two temperatures are equal.

In the text

Fig. A.1 Mass-loss rate versus the classical Eddington parameter Γe, comparing our predictions with those of Gräfener & Hamann (2008). Our model predictions are indicated by coloured symbols connected with dashed lines. The best-fit relation from Eq (10) is shown as solid lines, while the fit from Gräfener & Hamann (2008) is shown as dotted lines. Absolute mass-loss rates predicted by the two fits for M* = 60 M⊙ and T* = 38 kK are indicated by red and black star symbols, respectively.

In the text