© The Authors 2026

1 Introduction

Winds from hot, luminous stars provide a natural laboratory for testing complex radiative-transfer and hydrodynamic models. They play a central role in the evolution of massive stars, angular momentum loss, and the final stages of evolution (towards a neutron star or black hole). Moreover, such winds shape the stellar environment, injecting momentum, energy, and chemically processed material into the interstellar medium (ISM), thus regulating star formation and contributing to galactic chemical enrichment.

These winds are primarily driven by the transfer of momentum from the stellar radiation field to the plasma in the outer layers of the atmosphere, resulting in what are known as radiation-driven winds. These outflows were first described by Castor et al. (1975), whose work led to what is now known as the CAK theory. Subsequent extensions led to the modified CAK framework (m-CAK, Pauldrach et al. 1986; Friend & Abbott 1986), which models particle acceleration due to radiation pressure from a finite-size star. This theory also incorporates the effects of stellar rotation, applicable for rotation speeds up to 75% of the critical rotation velocity. Curé & Araya (2023) provided an updated overview of this theoretical framework.

A commonly adopted approach for modelling the winds of massive stars is to prescribe the velocity profile through the empirical ‘β-law’ (Pauldrach et al. 1986). Theoretical expectations place β in the 0.5–1 range, yet in practice it is treated as a free parameter and pushed to much higher values (β ≈ 2−3) to fit the Hα line profile observed in B supergiant stars (e.g. Markova & Puls 2008; Haucke et al. 2018). However, such extrapolation falls outside the wind outflow regime where the β-law is physically justified, and it turns the wind parameters (mass-loss rate and terminal velocity) into adjustable quantities rather than predictions from a hydrodynamic model. A more physically motivated alternative is to derive the wind structure directly by solving the hydrodynamic equations, thereby ensuring that the resulting wind parameters are consistent with the underlying physics instead of being imposed through an assumed velocity law.

Although the m-CAK theory reproduces the most observed wind features in the early-type stars, several key questions remain unresolved: the non-uniqueness of wind-structure solutions when modelling line spectral features (Venero et al. 2024), the lack of exact stationary solutions for the entire parameter space of the line-force multipliers (Venero et al. 2016), the existence and character of the bi-stability jump (e.g. Markova & Puls 2008; de Burgos et al. 2024), the interaction of slow and fast winds (Curé et al. 2011), the coupling between wind variability and stellar pulsations (Cidale et al. 2023), or the detailed properties of clumpy winds (e.g. Rubio-Díez et al. 2022). Having exact hydrodynamic solutions is therefore essential to address outstanding problems, as hydrodynamics yields accurate velocity and density profiles for studying wind variability and instabilities, and supplying robust initial conditions for time-dependent and multi-D wind models.

Within the framework of the m-CAK theory, Abbott (1982) parametrised radiative acceleration due to spectral lines using three key quantities: k, α, and δ. These parameters characterise, respectively, the effective number of contributing lines, the distribution of line strengths, and the sensitivity of the line-driving force to variations in the ionisation state. Depending on the value of the δ parameter, the hydrodynamic solutions may correspond to either fast (m-CAK–classical) or so-called δ-slow regimes (Curé et al. 2011). For instance, slowly rotating stars may develop either fast (for δ ≲ 0.2) or δ-slow (for δ ≳ 0.28) solutions, which differ markedly in their terminal wind velocities. However, for intermediate values of δ, a gap exists for a subset of line-force parameters where no stationary solutions connecting the fast and slow regimes have previously been found (Venero et al. 2016, 2024).

One of the primary objectives of this work is to model radiation-driven winds in 1D within the gap region. We employed the time-dependent hydrodynamic code ZEUS-3D (Clarke 1996, 2010), which is capable of converging towards stationary solutions by iteratively advancing from an initial configuration. This approach allowed us to probe and characterise the solutions emerging inside the gap region, thereby providing insights into the dynamics of radiation-driven winds in massive stars across this otherwise unexplored sector of parameter space. Once the hydrodynamic solutions within the gap region were obtained, we solved the radiative transfer equation under non-local thermodynamic equilibrium (non-LTE) conditions in the comoving frame, assuming spherical symmetry to compute synthetic line profiles. These profiles were then compared to quantify and assess the differences arising from the various wind solutions.

The paper is organised as follows: Sect. 2 reviews the stationary m-CAK theory, while Sect. 3 describes the numerical codes employed and the methodology followed. The resulting hydrodynamic solutions and the corresponding synthetic line profiles are presented in Sect. 4. Section 5 provides a discussion, and Sect. 6 summarizes our main conclusions.

2 The time-dependent m-CAK theory

This section provides a summary of the time-dependent m-CAK theory in the equatorial plane (one radial dimension) for radiation-driven winds (e.g. Feldmeier 1995; Araya et al. 2018). Following Araya et al. (2018), gas outflow from the outer layers of stars into space is mainly described by

• the mass-continuity equation,

  $$\frac{\partial \rho }{\partial t}+\frac{1}{r^2}\frac{\partial }{\partial r}(\rho r^{2}v)=0,$$(1)

• and the momentum equation,

  $$\frac{\partial v}{\partial t}+v\frac{\partial v}{\partial r}=-\frac{a^2}{\rho }\frac{\partial \rho }{\partial r}+g_{eff}(r)+g_{rad}^L(r,v,dv/dr),$$(2)

  where r, ρ(r, t), v(r, t), ∂v(r, t)/∂r, and a represent the radial coordinate, mass density, flow velocity, velocity gradient, and isothermal sound speed, respectively.

The effective gravity at the equatorial plane, g_eff, is given by

$$g_{eff}(r)=-\frac{G{M}_{*}(1-\Gamma )}{r^2}\left(1-{\Omega }^2\frac{R_{*}}{r}\right),$$(3)

where G is the gravitational constant, M_∗ is the mass of the star, R_∗ is its radius, and Ω = v_rot/v_crit is the ratio of the equatorial rotational velocity of the star to the critical velocity at the equator (0 < Ω < 1). The parameter Γ represents the Eddington factor that accounts for continuum radiative acceleration dominated by Thomson scattering.

Acceleration due to the spectral lines, $g_{\text{rad}}^L$, is given by

$$g_{rad}^L(r,v,dv/dr)=CF\frac{\Gamma G{M}_{*}k}{r^2}{\left(\frac{1}{{\sigma }_e{v}_{th}}\right)}^{\alpha }{\left(\frac{1}{\rho }\frac{\partial v}{\partial r}\right)}^{\alpha }{\left(\frac{n_e}{W(r)}\right)}^{\delta },$$(4)

where σ_e is the absorption coefficient due to Thomson scattering and v_th is the thermal speed of the protons. In addition, n_e is the electron density, and CF(r, v, ∂v/∂r) is the correction factor due to the finite stellar disk. In this work, we adopted the standard geometrical dilution factor given by

$$W(r)=\frac{1}{\text{2}}\left(1-\sqrt{1-\frac{R_{*}^2}{r^2}}\right),$$(5)

which accounts for the finite angular extent of the stellar disk as seen from radius r. The finite-disk correction factor CF was computed following the m-CAK formalism (e.g. Pauldrach et al. 1986), using

$$CF=\frac{{(1+\sigma )}^{1+\alpha }-{(1+\sigma {\mu }_{*}^2)}^{1+\alpha }}{(1+\alpha )\sigma {(1+\sigma )}^{\alpha }(1-{\mu }_{*}^2)},$$(6)

where ${\mu }_{*}=\sqrt{1-R_{*}^2/r^2}$ and σ = d ln v/d ln r − 1.

An alternative description of the radiative line force was introduced by Gayley (1995), who reformulated the CAK parametrisation in terms of the dimensionless line-strength parameter $\overline{Q}$ and the cut-off parameter Q₀, instead of the line-force parameter k. In this formalism, line acceleration $g_{\text{rad}}^L$ takes the form

$$g_{rad}^L(r,v,\text{d}v/\text{d}r)=CF\frac{\Gamma G{M}_{*}\overline{Q}}{(1-\alpha )r^2}{\left(\frac{1}{Q_0{\sigma }_{e}c}\right)}^{\alpha }{\left(\frac{1}{\rho }\frac{\partial v}{\partial r}\right)}^{\alpha }{\left(\frac{n_e}{W(r)}\right)}^{\delta }.$$(7)

The relation between the Gayley parameters and the standard CAK line-force parameters can be written as

$$k=\frac{1}{1-\alpha }{\left(\frac{v_{th}}{c}\right)}^{\alpha }\overline{Q}{Q}_0^{-\alpha }.$$(8)

We adopted $Q_0\approx \overline{Q}$, treating both parameters as representative of the same underlying line-strength distribution. This allowed us to derive $\overline{Q}$ from the adopted values of k and α.

The line-force parameters α, k, and δ define the analytical parametrisation of the line-driving acceleration in the m-CAK theory. Especially relevant for our work is the parameter δ, introduced by Abbott (1982), which accounts for changes in ionisation throughout the wind.

The m-CAK theory delivers three physical stationary solutions (Curé & Araya 2023) depending on δ and the rotation rate, Ω:

• The fast solution has a slow rotation rate (Ω ≲ 0.75) and almost no changes in the ionisation of the wind with distance (i.e. δ ≲ 0.2). This is known as the classical m-CAK solution (Pauldrach et al. 1986; Friend & Abbott 1986) and predicts a high terminal velocity.

• The Ω-slow solution corresponds to models with high rotation rates (Ω ≳ 0.75) and yield a low terminal velocity (Curé 2004).

• The δ-slow solution results from high values of δ (usually δ ≳ 0.28) and leads to low values for the terminal velocity (Curé et al. 2011).

The fast and δ-slow solutions are separated by a region in the δ-parameter space, or ’gap’ (Venero et al. 2016), where no steady hydrodynamic solution could be found using the Hydwind code (Curé 2004). As the Hydwind code solves hydrodynamic equations under spherical symmetry and in a stationary state, both fast and δ-slow solutions can be achieved without significant convergence issues for most radiation force parameters. However, when δ values approach the critical region (or gap), the solutions fail to converge. This gap interrupts the smooth transition of hydrodynamic regimes along the domain of solutions. To find the solution within the gap, we employed the time-dependent code ZEUS-3D (Clarke 1996, 2010), following the scheme developed by Araya et al. (2018).

3 Methodology

3.1 Time-dependent hydrodynamic models

The code ZEUS and its subsequent extensions, ZEUS-2D and ZEUS-3D, have been among the most influential numerical tools in computational astrophysics since their introduction by Stone & Norman (1992) and Stone et al. (1992). Its robust algorithms for solving the equations of hydrodynamics and magnetohydrodynamics, together with its flexibility for multidimensional simulations, enabled pioneering studies of numerous astrophysical flows. The code has been widely applied to problems such as accretion disks and magnetorotational turbulence, protostellar jets and outflows, turbulence in the ISM, star formation, and the hydrodynamics of stellar winds (e.g. Hawley et al. 1995; Ouyed & Pudritz 1997; Fujita et al. 2003; Barai et al. 2012; Venero et al. 2016; Araya et al. 2018).

The ZEUS-3D code employs an Eulerian finite-difference, operator-split scheme with consistent advection for transport (Clarke 1996, 2010). Line acceleration was implemented following the parametrisation of Gayley (1995), as adapted by Araya et al. (2018). The radial velocity gradient was computed using a first-order forward finite difference, and the resulting line acceleration was then added explicitly as a source term in the momentum equation.

At the inner boundary, we imposed the base density, while the inflow velocity was allowed to float. It is typically obtained by linearly extrapolating the velocity from the two innermost grid zones of the computational domain. In practice, the base density was adjusted in order to obtain a stable outflow and allow the hydrodynamic solution to converge towards a stationary state. At the outer boundary, we imposed a standard outflow (zero-gradient) condition. No fixed terminal velocity was enforced. The simulations were performed in one radial dimension using a grid of 500 points extending from the stellar radius up to 100 R_∗. The mesh is non-uniform, employing logarithmic radial coordinates to enhance the grid resolution near the stellar surface, where strong gradients are expected, following the prescription of Araya et al. (2018).

To initialise ZEUS-3D, we supplied converged stationary solutions, obtained from the Hydwind code, adopting the required stellar and line-force parameters. Our procedure is as follows. Utilising the Hydwind code, we first sought the position of the gap in the δ-parameter space. Starting from a low δ value in the fast wind regime, we increased δ until reaching the highest value that achieves a converged solution, marking the lower boundary of the gap, hereafter referred to as δ_min. Then, a δ value corresponding to the slow regime was chosen (for example, δ = 0.4), which was gradually decreased. The final δ value at which a converged solution was achieved defines the upper boundary of the gap, denoted as δ_max. Thus, for a given set of stellar and line-force parameters (α, k), together with the rotation rate Ω, the gap is defined as the interval δ_min < δ < δ_max.

ZEUS-3D solves the time-dependent hydrodynamic equations and, unlike Hydwind, does not impose the critical-point condition to construct a stationary solution. Instead, the equations are integrated numerically from the stellar surface outwards, and the system evolves in time until a steady state is reached. In this approach, the final solution is determined solely by the adopted stellar and line-force parameters.

To illustrate the relaxation towards the stationary regime, Fig. 1 shows the temporal evolution of the velocity field v(r) computed with ZEUS-3D. The curves correspond to different times expressed in units of the flow crossing time, defined as t_flow = R_max/v_∞, where R_max is the outer radius of the computational domain and v_∞ is the terminal velocity of the wind. For the model shown here, characterised by T_eff = 18 000 K, log g = 2.5, R_∗ = 23 R_⊙, Ω = 0.27, α = 0.515, k = 0.104, δ = 0.1, and $\overline{Q}=67.25$, we adopted a base wind density of ρ₀ = 10⁻¹² g cm⁻³. For consistency in the comparison, the stationary solution computed with Hydwind was obtained using the same base density, instead of the usual condition in which the stellar radius is defined at the optical depth τ = 2/3. For this model, the flow crossing time is approximately t_flow ≈ 4 × 10⁶ s. The figure displays the velocity profile at 0.1, 0.4, 1, and 10 t_flow. At early times, the inner wind regions already approach their stationary structure, while the outer layers are continue to evolve towards a steady state. In particular, within approximately r ≲ 2R_∗, the solution is already close to the final profile at 0.1 t_flow. As the simulation proceeds, the outer wind gradually relaxes, and by t ≈ t_flow the entire domain essentially reaches the stationary configuration. The simulation nevertheless continues up to 10 t_flow, which largely exceeds the time required for convergence. For comparison, Fig. 1 also includes the stationary solution computed with Hydwind for the same stellar and line-force parameters. The time-dependent solution obtained with ZEUS-3D clearly converges towards the same velocity profile, confirming the consistency between both approaches.

In some cases, transient features appear during the relaxation phase. For instance, when the initial condition corresponds to a wind solution belonging to a different regime, the intermediate velocity profiles display temporary distortions while the flow readjusts to the correct stationary solution (see discussion in Araya et al. 2018). These features disappear once the entire domain relaxes and does not affect the final converged state, provided that the simulation is evolved for a sufficiently long time.

To verify that the final solution is independent of the initial velocity profile, we used different initial conditions, including solutions from different wind regimes. In all cases, the simulations converge towards the same stationary solution determined by the input stellar and line-force parameters. This property is particularly relevant for exploring the gap region in parameter space, where no Hydwind solution is available.

The time integration in ZEUS-3D uses explicit time stepping with a Courant–Friedrichs–Lewy coefficient of 0.5. In practice, the simulations were evolved well beyond the time required for convergence. Specifically, the equations were integrated from t = 0 up to t_max = 9 × 10⁷ s, using a fixed time step of ∆t = 9 × 10⁴ s, resulting in 1000 temporal steps in total. Rather than identifying the exact moment at which the steady state is reached, we simply evolved the system for a sufficiently long time and adopted the final snapshot as the stationary solution. Convergence was further verified by checking that the mass-loss rate remained constant throughout the computational domain.

Fig. 1 Temporal evolution obtained with ZEUS-3D (solid lines) from an initial solution (dashed line). The times corresponding to the different velocity laws are given in units of the flow crossing time.

3.2 Radiative transfer calculations for spectral lines

The atmospheric structure of the star consists of a photospheric temperature and density stratification smoothly connected to an isothermal stellar wind. The transition between the two regions was implemented through a continuous interpolation. Subsequently, following the hydrodynamic calculations described above, we computed synthetic spectra to assess which stationary solutions most accurately reproduced the observed line profiles. To this end, we solved the radiative transfer equation in an expanding atmosphere under non-LTE conditions, assuming spherical symmetry and stationarity in the comoving frame (following Mihalas & Kunasz 1978). The radiative transfer and statistical equilibrium equations are iterated simultaneously until convergence was achieved.

These calculations were performed with the MULITAS code (MUlti LIne Transfer for Active Stars), which includes a multilevel radiative transfer treatment for He II and Mg II ions in the comoving frame (Venero et al. 2000; Cidale 1998), relevant for massive-star winds. The MULITAS code was also adapted by Cidale (1993) to model the Si IV atom (six bound levels plus continuum), allowing us to compute its ultraviolet (UV) resonance doublet.

In this work, we extended the MULITAS code to include a multilevel He I atom, incorporating both singlet and triplet systems. A total of 23 bound levels plus a continuum were implemented. Energy levels and oscillator strengths for radiatively permitted transitions were obtained from the NIST database (Kramida et al. 2024). We also incorporated the corresponding collisional rates (bound–bound and bound–free) and radiative cross-sections for the relevant levels (cf., Rohrmann 2001). This extension enables the computation of He I lines in both the optical and infrared (IR), complementing the UV diagnostics provided by Si IV.

To model the hydrogen lines, we employed the Appel code (Catala & Kunasz 1987; Cidale & Ringuelet 1993), which solves the radiative transfer equation in the comoving frame for continua and selected spectral lines under the same physical assumptions of sphericity and stationarity. Appel provides the H line profiles corresponding to the hydrodynamic velocity fields obtained from ZEUS-3D and is used in the same manner as MULITAS for the calculation of Si IV and He I lines.

Together, these radiative transfer calculations enabled us to test the hydrodynamic models across the UV, optical, and IR spectral ranges, using a consistent set of stationary wind structures. This approach also allowed us to assess which wind solution produces the most distinct spectral signatures.

4 Results

To model a B supergiant star, we adopted the T19 model from Venero et al. (2016) as a representative case to explore the new set of hydrodynamic solutions. The parameters of this model are: T_eff = 19 000 K, log g = 2.5, R_∗ = 40 R_⊙, α = 0.5, k = 0.32, and $\overline{Q}=433.1$. In addition, several values of the rotation rate Ω and the line-force parameter δ were considered to sample the full family of solutions, including the fast regime, the δ-slow regime, and the transition from fast to δ-slow solutions (i.e. in the gap region).

4.1 Hydrodynamic solutions

Figure 2 displays a sequence of hydrodynamic solutions for model T19, computed with the ZEUS-3D code, considering rotation rates of Ω = 0.0, 0.2, 0.4, and 0.6. In the plot, each curve corresponds to a different value of δ, ranging from 0.0 to 0.4 with a step of 0.01. The solutions with the highest terminal velocities correspond to the fast regime, while those with the lowest terminal velocities correspond to the δ-slow regime. In between, a new family of solutions emerges; these are the new solutions reported in this work (dashed lines in the plot).

Despite the uniform sampling in δ, the transition-region solutions exhibit noticeably wider spacing compared to both the fast and δ-slow regimes. This effect becomes more evident as the rotation rate increases. The widening of the spacing reflects the rapid structural changes occurring in the wind as the flow transitions between the two classical regimes, highlighting the sensitive dependence of the solution on the ionisation parameter in this interval.

These new solutions form a continuous transition between the fast and δ-slow regimes. Near the photosphere, they behave similarly to the fast solutions, but as the distance from the star increases, they gradually resemble the δ-slow wind regime. This transition introduces a change in the slope dv/du (where u = − R_∗/r), which produces a distinct kink in the velocity profile. When comparing different rotation rates, this kink becomes steeper as the stellar rotation increases. However, at higher values of Ω, the situation becomes more critical: the flow approaches conditions where the Ω-slow regime is expected to appear (a regime beyond the scope of this work). In addition, B supergiants are generally slow rotators, so such high rotation rates are not expected to be representative of these stars.

A key result is that these kinks are stable structures of the wind, preserved even after long temporal integrations. As the simulations evolve for longer times, the solutions retain their shape, including the kink, confirming that they correspond to stationary states. Moreover, as δ increases, the kink appears progressively farther from the stellar surface. Therefore, if the ionisation conditions were to vary (see Sect. 5.3 for a discussion), effectively modifying the value of δ, the kink may propagate outwards through the stellar wind. This behaviour could be related to wind variability and to the formation and evolution of discrete absorption components, as reported in OB stars (e.g. Kaper & Henrichs 1994; Kaper et al. 1999).

Figure 3 presents the terminal velocity and mass-loss rate as functions of the δ parameter for all rotation rates considered (from 0.0 to 0.6, in steps of 0.1). The filled symbols mark the fast and δ-slow solutions, while the open circles denote the intermediate solutions found within the transition region. The left panel shows the variation in terminal velocity with δ, while the right panel displays the corresponding mass-loss rates. The result predicts that a decrease in terminal velocity is related to an increase in mass-loss rate. This figure extends the results of Venero et al. (2016), now fully covering the region where no stationary solutions were previously found. The ZEUS-3D simulations fill this gap, providing a continuous sequence of hydrodynamic solutions that smoothly bridge the fast and δ-slow regimes. The previously missing intermediate region is now fully resolved, revealing that the transition between both regimes is continuous, both in the structure of the flow (through the velocity profiles) and in the resulting global wind parameters (v_∞ and Ṁ).

Furthermore, Table A.1 lists the terminal velocities and mass-loss rates predicted by both Hydwind and ZEUS-3D for the different δ values, taking Ω = 0.2 as a representative case. While Hydwind shows no solutions within the interval 0.21 < δ < 0.29, ZEUS-3D successfully finds stationary solutions for all δ values. Although we have explored only a single model, it is important to emphasize that both codes yield similar terminal velocities and mass-loss rates in the fast regime. In contrast, in the δ-slow regime, the predicted mass-loss rates differ by up to a factor of two, with ZEUS-3D systematically yielding lower values.

Fig. 2 Hydrodynamic solutions for the T19 model (Teff = 19 000 K), adopting values from 0 < δ < 0.4 and ∆δ = 0.01, using different rotational rates, obtained with the ZEUS-3D code. The dashed lines represents the solutions in the gap region.

4.2 Synthetic profiles

Figure 4 presents the synthetic line profiles computed for three representative hydrodynamic solutions from model T19: a fast solution (δ = 0.00), a transition-region solution (δ = 0.26), and a δ-slow solution (δ = 0.35), all computed for Ω = 0. The first column shows the corresponding velocity profiles, while the remaining columns display the synthetic profiles for Si IV doublet λλ1394, 1403 Å, He I 5876 Å, Hα, and the He I line + Brα blend near 4.05 µm.

Across all wavelengths, the profiles form a coherent morphological sequence that reflects the continuity of the underlying hydrodynamic structure. The transition-region solution consistently shows intermediate behaviour, matching its position between the fast and δ-slow regimes.

The Si IV doublet displays the expected P Cygni morphology in all three cases. As the terminal velocity decreases from the fast to the δ-slow solution, the width of the blue absorption trough becomes progressively narrower, while the emission component also weakens along the same sequence, even when the mass loss rate increases from 1.06 × 10⁻⁶ to 2.35 × 10⁻⁶ M_⊙yr⁻¹.

The He I λ5876 Å line begins as a pure absorption feature in the fast solution and gradually develops a weak red emission wing for the δ-slow solution. Although the variation is modest, the intermediate profile is clearly bracketed by the fast and δ-slow solutions, providing a clean spectroscopic counterpart to the underlying velocity structure and mass-loss rate.

The Hα line shows a much stronger response to changes in the wind structure. The fast solution remains dominated by absorption, while the transition model already forms a clear P Cygni profile. In the δ-slow case, the emission lobe becomes very strong, consistent with the substantial increase in mass-loss rate (from 1.06 × 10⁻⁶ to 2.35 × 10⁻⁶ M_⊙yr⁻¹).

In the IR, both the He I λ4.04 µm line and Brα respond sensitively to the structural changes induced by increasing δ. The He I line evolves from pure emission to a well-defined P Cygni shape in the δ-slow regime. The Brα line evolves from absorption to a strong P Cygni profile, with the emission peak dominating in the slow-wind case. When combined, both lines appear largely in emission for the δ-slow solution, as the He I emission masks the blue absorption component of the Brα P Cygni profile. For instance, both lines were reported in emission in the B supergiant 55 Cygni (Cidale et al. 2023).

Overall, all diagnostics consistently place the transition-region solution between the fast and δ-slow regimes, reinforcing the continuity of the physical sequence revealed by the ZEUS-3D models. This behaviour indicates that the transition solution naturally connects both regimes within a continuous family of wind solutions.

Fig. 3 Terminal velocity (left panel) and mass-loss rate (right panel) as functions of the δ parameter for the hydrodynamic solutions obtained with ZEUS-3D for the T19 model, for several rotation rates. The filled symbols correspond to the fast and δ-slow solutions, while the open circles denote the solutions found within the transition region.

Fig. 4 Synthetic line profiles computed for the three hydrodynamic wind regimes using the T19 model: the fast solution (δ = 0, top row), transition solution (δ = 0.26, middle row), and δ-slow solution (δ = 0.35, bottom row). The first column indicates the hydrodynamic solution obtained with ZEUS-3D, followed by the synthetic profiles for Si IV (second column), He I 5876 Å (third column), Hα (fourth column), and the He I + Brα blend around 4.05 µm (fifth column). The fast solution corresponds to Ṁ = 1.06 × 10−6 M⊙ yr−1 and v∞ = 720 km s−1; the solution in the transition region to Ṁ = 1.67 × 10−6 M⊙ yr−1 and v∞ = 322 km s−1; and the δ-slow solution to Ṁ = 2.35 × 10−6 M⊙ yr−1 and v∞ = 253 km s−1.

5 Discussion

5.1 Why the Hydwind code fails to find solutions within the gap

Under particular values of the radiation-force parameters, especially for δ values within the interval (δ_min, δ_max), the Hydwind code is unable to obtain a physically consistent wind solution. This limitation arises mainly from two factors. First, the inner boundary condition in the Hydwind code is imposed at an optical depth τ = 2/3, or equivalently at the corresponding mass density at that layer. As a consequence, once the inner boundary condition is fixed, the velocity gradient dv/dr cannot smoothly pass through the critical point. Since the inner boundary conditions cannot be varied, the solution may therefore fail to satisfy the regularity and continuity conditions at the critical point. Although one could in principle adjust the boundary conditions to enforce these constraints, doing so within a stationary framework becomes impractical.

A time-dependent approach, however, offers an important advantage: the flow can dynamically adjust and relax towards a physically consistent solution, allowing the system to naturally approach the critical configuration. Nevertheless, the ZEUS-3D code does not always converge to a hydrodynamic solution. This issue was resolved by lowering the mass density at the inner boundary condition by one order of magnitude with respect to the photospheric density obtained from the Hydwind solution, where the photosphere was defined by imposing τ = 2/3. Once a solution was obtained, we identified the radius at which this density is reached (typically around 1.01 R_∗) using an appropriate Tlusty or Kurucz atmospheric model (Hubeny & Lanz 1995; Kurucz 1979) matching the stellar parameters. Then, we reran the ZEUS-3D simulations imposing the boundary condition at that corresponding photospheric layer. In this way, we ensured that the photospheric density stratification and the wind were smoothly connected.

5.2 Comparison of line profiles computed with the traditional β-law and the hydrodynamic models

Figure 5 shows the line profiles computed with the hydrodynamic solutions (obtained with ZEUS-3D) together with those calculated from β-laws with β = 1 and β = 2. In all cases, the synthetic spectra were computed by adopting the same mass-loss rate, terminal velocity, and initial velocity as those of the corresponding hydrodynamic model.

Reproducing the Hα line profiles typically observed in the winds of B supergiants generally requires β values in the ∼1.5−2.5 range (e.g. Markova & Puls 2008; Haucke et al. 2018). However, such velocity laws deviate significantly from any hydrodynamic wind solution. Nevertheless, a degeneracy exists when reproducing the overall shape of line profiles when different velocity laws are adopted (Venero et al. 2024). In other words, both the traditional β-law and hydrodynamic solutions produce similar spectral signatures under certain conditions. Therefore, comparing the synthetic line profiles obtained from both velocity prescriptions provides a useful way to evaluate their impact on the emergent spectra and to assess which hydrodynamic regimes can account for the high β values inferred from observations.

It is worth noting that the Hα profiles computed with a high β value resemble the profile obtained from the δ-slow hydrodynamic solution. This may have implications for previous studies based on β-laws, as similar Hα profiles could be associated with different underlying wind solutions, but with the same mass-loss rate and terminal velocity, although a more systematic analysis would be required to assess the generality of this result.

On the other hand, the IR He I emission at 4.049 µm (typically seen in BSG stars; see e.g. Cidale et al. 2023) can, in some cases, only be reproduced by adopting β-laws with a low β value, while the hydrodynamic solutions are able to reproduce this feature more naturally, indicating a higher sensitivity of this diagnostic to the adopted velocity structure. At the same time, the resulting emission predicted by each velocity law shows a different qualitative behaviour depending on the spectral range considered. These results reflect the interplay between velocity structure and mass-loss rate in shaping the emergent line profiles as well as the need to combine diagnostics from different spectral ranges to better describe the wind dynamics.

5.3 The ionisation structure of the wind

There is growing evidence for anomalous ionisation conditions in the winds of B supergiants. One example is the bi-stability jump, where both the mass-loss rate and terminal velocity change for stars with effective temperatures around 21 000 K (Lamers et al. 1995; Markova & Puls 2008). Additional evidence comes from the superionisation effect observed in the UV, which has led to the suggestion of hot regions in the wind and of X-ray emission (Cassinelli et al. 1981; Krtička & Kubát 2016). At the same time, optical spectra show lines from relatively low ionisation stages, and wind-clumping models are often required to explain the resulting opacity structure to fit the line features observed in different spectral ranges.

At present, no self-consistent model demonstrates how the ionisation parameter δ varies with spectral type, density distribution, and radial distance in the wind. Nevertheless, the fact that the empirical β-law required to reproduce the observed line profiles resembles a slow wind solution (i.e. associated with relatively high δ values) suggests that the ionisation structure of the wind may play an important role.

Although most applications of the standard m-CAK formalism adopt relatively low δ values, higher values are not excluded on physical grounds. Analytical considerations indicate that δ ≳ 1/3 can arise in media where neutral hydrogen behaves as a trace element, implying strong ionisation variations throughout the wind (Puls et al. 2000). In addition, approximate non-LTE calculations suggest that δ may approach unity in low-density winds or in low-metallicity environments (Kudritzki 2002). In this context, the range of δ values explored in this work can be regarded as physically plausible, particularly in view of the complex and not yet fully constrained ionisation structure of B supergiant winds. This indicates that further work in this direction is needed, especially using models that relax the assumption of isothermal winds.

Fig. 5 Comparison between synthetic line profiles computed using the hydrodynamic solutions obtained with ZEUS-3D (using the T19 model) and those derived from standard β-type velocity laws. The colour coding of the hydrodynamic solutions is the same as in Fig. 4: fast solution (δ = 0), transition solution (δ = 0.26), and δ-slow solution (δ = 0.35). The dashed lines correspond to the profiles computed adopting β laws with β = 1 (green) and β = 2 (red). From left to right: Hydrodynamic velocity structure, followed by the synthetic profiles for Si IV, Hα, and the He I + Brα blend around 4.05 µm.

6 Conclusions

In this work, we modelled radiation-driven winds of massive stars within the m-CAK framework, focusing on the transition region between the fast and δ-slow regimes. Using time-dependent ZEUS-3D simulations, we demonstrated that B super-giant winds admit a continuous sequence of stationary hydrodynamic solutions that extends beyond the traditionally recognised regimes. The region previously thought to lack stationary solutions is now fully resolved: we uncover a well-defined set of intermediate solutions that smoothly connects the fast and δ-slow regimes. Several of these solutions develop a distinct kink in the velocity law, whose location shifts outwards as δ decreases and becomes steeper with increasing rotation. These results show that the apparent discontinuities obtained with stationary approaches do not represent a physical gap in the parameter space, but rather the limitations of earlier methods.

We also computed synthetic line profiles for H I, He I, and Si IV by solving the radiative transfer in the comoving frame for different hydrodynamic solutions to illustrate the effect of the wind structure. A key contribution of this work is the inclusion of a detailed atomic model for these elements within a moving-medium radiative transfer framework, enabling self-consistent predictions of UV and optical lines together with a comprehensive set of IR transitions.

The resulting profiles reflect the underlying hydrodynamic structure: the solutions in the transition region produce intermediate line morphologies between the fast and δ-slow winds, providing direct spectroscopic evidence for the continuity of the hydrodynamic sequence. To further assess the impact of the adopted velocity structure, we compared these profiles with those computed using standard β-type velocity laws. In the model explored, Hα profiles computed using standard β-type velocity laws, with β > 1, closely resemble the profile obtained from the δ-slow hydrodynamic solution. This suggests that similar Hα signatures may arise from different underlying wind structures, when the mass-loss rate and terminal velocity are the same, with potential implications for previous studies based on β-law prescriptions. However, a more systematic exploration of the parameter space is required to assess its general validity.

These results establish a unified picture in which the wind properties of B supergiants evolve smoothly with the ionisation parameter δ, without invoking additional physical mechanisms to reproduce the observed line morphologies. The existence of this continuous sequence opens the door to a more refined interpretation of spectroscopic diagnostics across the UV–optical–IR range, suggesting that the hydrodynamic state of B supergiant winds can be robustly constrained through simultaneous multi-wavelength observations. Furthermore, perturbations in the wind ionisation alter the line force and may trigger transitions between different hydrodynamic regimes. Such ionisation-driven regime changes offer a plausible explanation for structured and variable winds, including the emergence of large-scale features such as discrete absorption components. A systematic investigation of these effects will be the subject of future work.

Acknowledgements

We thank the referee, Nicolas Moens, for the careful reading and constructive comments, which helped to improve the clarity and quality of this work. M.C.F., R.O.J.V., and L.S.C. acknowledge financial support from CONICET (PIP 11220200101337CO) and from the Universidad Nacional de La Plata through the Programa de Incentivos (grant 11/G192). R.O.J.V. also acknowledges support from grant 11/G193. I.A. and M.C. thank the FONDECYT projects 1230131 and 1261498 for their support. This work has been partially co-funded by the European Union, Project 101183150 - OCEANS.

References

Appendix A Terminal velocity and mass-loss rates for the T19 model (Ω = 0.2)

All Tables

All Figures

	Fig. 1 Temporal evolution obtained with ZEUS-3D (solid lines) from an initial solution (dashed line). The times corresponding to the different velocity laws are given in units of the flow crossing time.
In the text

	Fig. 2 Hydrodynamic solutions for the T19 model (Teff = 19 000 K), adopting values from 0 < δ < 0.4 and ∆δ = 0.01, using different rotational rates, obtained with the ZEUS-3D code. The dashed lines represents the solutions in the gap region.
In the text

	Fig. 3 Terminal velocity (left panel) and mass-loss rate (right panel) as functions of the δ parameter for the hydrodynamic solutions obtained with ZEUS-3D for the T19 model, for several rotation rates. The filled symbols correspond to the fast and δ-slow solutions, while the open circles denote the solutions found within the transition region.
In the text

Fig. 4 Synthetic line profiles computed for the three hydrodynamic wind regimes using the T19 model: the fast solution (δ = 0, top row), transition solution (δ = 0.26, middle row), and δ-slow solution (δ = 0.35, bottom row). The first column indicates the hydrodynamic solution obtained with ZEUS-3D, followed by the synthetic profiles for Si IV (second column), He I 5876 Å (third column), Hα (fourth column), and the He I + Brα blend around 4.05 µm (fifth column). The fast solution corresponds to Ṁ = 1.06 × 10−6 M⊙ yr−1 and v∞ = 720 km s−1; the solution in the transition region to Ṁ = 1.67 × 10−6 M⊙ yr−1 and v∞ = 322 km s−1; and the δ-slow solution to Ṁ = 2.35 × 10−6 M⊙ yr−1 and v∞ = 253 km s−1.

In the text

Fig. 5 Comparison between synthetic line profiles computed using the hydrodynamic solutions obtained with ZEUS-3D (using the T19 model) and those derived from standard β-type velocity laws. The colour coding of the hydrodynamic solutions is the same as in Fig. 4: fast solution (δ = 0), transition solution (δ = 0.26), and δ-slow solution (δ = 0.35). The dashed lines correspond to the profiles computed adopting β laws with β = 1 (green) and β = 2 (red). From left to right: Hydrodynamic velocity structure, followed by the synthetic profiles for Si IV, Hα, and the He I + Brα blend around 4.05 µm.

In the text