© The Authors 2025

1 Introduction

Massive stars play an important role in understanding the formation and structure of the universe (Davidson & Humphreys 2012; Langer 2012; Vink et al. 2015; Kashi et al. 2016; Sabhahit et al. 2022; Farrell et al. 2022; Mukhija & Kashi 2024, 2025). Due to their enormous masses and powerful stellar winds, these stars significantly contribute to the chemical enrichment of the universe by synthesising heavy elements (Langer 2012; Hurley et al. 2013; Nomoto et al. 2013; Johnson 2019; Volpato et al. 2024). The mechanism driving these powerful winds was first described by Castor et al. (1975) and referred to as the CAK-theory of radiation-driven winds. According to this theory, the radiative acceleration in the atmosphere arises from a multitude of spectral absorption lines, leading to the outward expulsion of material (Kudritzki & Puls 2000; Puls et al. 2008). Later on Friend & Abbott (1986); Pauldrach et al. (1986) introduced the modified-CAK (m-CAK) model, which incorporates two key enhancements: it accounts for the finite size of the star through the finite disk correction factor, and it includes the effects of stellar rotation, specifically the centrifugal force resulting from rapid rotation. This m-CAK model has successfully described both terminal wind velocities (V_∞) and mass-loss rates (Ṁ) observed in massive stars (Puls et al. 2008). They found that stellar rotation has a significant impact on the observable wind parameters of massive stars. However, they could not find wind solutions for stars rotating at high rotational speeds (i.e. when υ_rot¹ is more than 75% of its critical velocity, υ_crit). Stars with such rotation rates (e.g. Oe, Be, Bn, WR, among others) are categorised as rapid rotators, as their high rotational speed significantly influences their wind dynamics and evolution (Lamers & Pauldrach 1991; Bjorkman & Cassinelli 1993; Owocki et al. 1996; Maeder & Meynet 2000; Townsend et al. 2004; Curé & Araya 2023).

Curé (2004) revisited the effect of stellar rotation on radiation-driven winds using 1D stationary models and identified a novel wind solution that emerges when the rotational speed exceeds ~75% of the critical limit. This solution, referred to as the Ω-slow solution (here Ω is a dimensionless rotation parameter), is defined as the ratio of the stellar angular velocity to the critical angular velocity, $$\Omega =\frac{{\omega }_{\text{rot}}}{{\omega }_{\text{crit}}}=\frac{{\upsilon }_{\text{rot}}\left(\Omega ,\text{eq}\right)}{{\upsilon }_{\text{crit}}}\frac{R_{\text{eq}}^{\max }}{R\left(\Omega ,\text{eq}\right)},$$(1)

where υ_rot(Ω, eq) is the equatorial rotational velocity, R(Ω, eq) is the equatorial radius for a given Ω, and $R_{\text{eq}}^{\max }=3R_{\text{pole}}/2$ is the maximum equatorial radius at critical rotation. A slower and denser wind characterises the Ω-slow solution compared to the standard m-CAK solution (henceforth referred to as the fast solution). It predicts significantly enhanced densities near the equator, particularly when combined with bistability effects Lamers & Pauldrach (1991), see also Björklund et al. (2021, 2023); de Burgos et al. (2024), resulting in a density contrast of up to 10⁴ between the equatorial and polar regions. This suggests that the Ω-slow solution could play a critical role in the formation of equatorial disks observed in B[e] stars (Curé et al. 2005). Madura et al. (2007) conducted detailed time-dependent hydro-dynamic simulations to investigate the behaviour of radiation-driven winds in rapidly rotating stars. Their findings suggested the inviability of the Ω-slow solution as a mechanism to reproduce the high equatorial densities inferred for B[e] star disks. They also revealed the presence of abrupt kink transitions in the wind velocity profiles at rotational speeds ranging from 75% to 85% of the critical rotational velocity.

Araya et al. (2018) investigated the transition region in a 1D equatorial model of a rapidly rotating star, focusing on where the fast radiation-driven wind solution switches to the Ω-slow solution at high rotational speeds. They found a co-existing region inside a small interval of Ω, where both solutions are simultaneously present at the equatorial plane, and satisfy the same boundary conditions at the stellar surface. Their calculations were based on the stationary code HYDWIND (Curé 2004) and time-dependent hydrodynamical code ZEUS−3D (Clarke 1996). They suggested that the topology of the non-linear m-CAK differential equation predicts two wind solutions as a function of Ω for a given set of line force parameters. The fast solution ceases to exist when Ω is above a certain threshold value that depends on the stellar and line force parameters. In contrast, the Ω-slow solution begins to exist from this interval up to Ω = 1. They also showed that time-dependent solutions are highly sensitive to the topology branch to which the initial velocity profile belongs. If the initial velocity profile deviates significantly from the steady-state asymptotic solution, the calculations can result in non-monotonic ‘kink’ solutions. To achieve globally monotonic acceleration and obtain either the fast or Ω-slow solutions, an initial velocity profile must accurately represent the respective fast or slow regime (topology branch). They also revealed that higher velocities and lower densities characterise the fast solution regime, whereas the Ω-slow solution regime showed lower velocities and higher densities. They also found that stellar and line-force parameters influenced the position and width of the co-existence region. An increase in the value of the line-force parameter α shifts the co-existence region to higher Ω values, while the width of the region remains unchanged. Additionally, they did not account for gravity-darkening effects and the distortion of the star’s shape due to high rotational speeds; both have significant implications for the structure and dynamics of the wind.

When accounting for high rotational velocities, a star becomes oblate, adopting a shape resembling a rotating ellipsoid. According to von Zeipel (1924), the radiative flux F at a given co-latitude on a rotating star is directly proportional to the local effective gravity, i.e. T_eff ∝ 𝑔_eff^1/4. The oblateness caused by rapid rotation modifies the local effective gravity, which is expressed as 𝑔_eff = 𝑔_grav + 𝑔_rot. Consequently, this variation in 𝑔_eff alters the surface temperature distribution of the star, with polar regions becoming hotter and equatorial regions cooler, a phenomenon referred to as Gravity Darkening (hereafter GD, see von Zeipel 1924). Recent interferometric observations of nearby rapidly rotating stars by Berger et al. (2005); Aufdenberg et al. (2006); van Belle et al. (2006); Zhao et al. (2010) indicate that GD does not conform well to von Zeipel’s model, which appears to overestimate the temperature difference between the poles and the equator. This discrepancy is often expressed with a power law, $T_{\text{eff}}\propto g_{\text{eff}}^{\beta }$, where the exponent β is less than 1/4. The reduced variation in temperature with latitude, compared to von Zeipel’s model, is traditionally attributed to a thin convective layer at the star’s surface, as described by Lucy (1967), who found β = 0.8 for stars with convective envelopes. Similar observations were made by Lovekin et al. (2006); Espinosa Lara & Rieutord (2007, 2011); Gagnier et al. (2019) in their 2D numerical models of rotating stars.

In this study, we extend the work of Araya et al. (2018) by examining the co-existence region where the fast solution transits to the Ω-slow solution, in the presence of the GD and oblateness. Our analysis focuses on a 1D equatorial model with high rotational speeds. Initially, we recalculated the co-existing interval, consistent with the findings of Araya et al. (2018). We then explored wind solutions while accounting for the influence of GD and oblateness effects. Furthermore, we also explored the dependence of the co-existing region on various stellar and line-force parameters.

This paper is organised as follows: Section 2 introduces our hydrodynamic wind model, starting with our assumptions. Then, we outline the time-dependent and steady-state equations used to describe the rotating m-CAK theory. In Section 3, we solve the stationary equations in the star’s equatorial plane for a wide range of rotational speeds where fast and Ω-slow solutions are achieved in the presence of GD and oblateness. In addition, a co-existence region is found where both solutions are simultaneously present. The influence of the line-force parameters on the co-existing regions is also studied. Finally, in Section 4, we present the discussion and summary of our findings.

2 1D steady state wind solutions for fast rotating stars

The work of Araya et al. (2017) thoroughly developed the complete set of equations governing the m-CAK radiation-driven wind theory, including the oblate finite disk correction factor and the gravity darkening effect models from von Zeipel (1924) and Espinosa Lara & Rieutord (2011), and incorporating these effects into the hydrodynamic wind model implemented in the HYDWIND code. It is worth noting that, for the m-CAK theory, we have also implemented a time-dependent version of the Equation of Motion (EoM) using the ZEUS-3D code (Clarke 1996, 2010, and references therein). Although these two numerical methods differ significantly (HYDWIND solves the stationary m-CAK EoM, and ZEUS-3D handles its time-dependent version), the steady-state solutions they produce are in close agreement. Given this consistency, and considering that HYDWIND is explicitly developed for the m-CAK theory, we decided to implement gravity darkening and stellar oblateness effects within this code. While ZEUS-3D is a well-established, versatile astrophysical (Magneto-)hydrodynamic code widely used for various applications, HYDWIND’s specialisation makes it our preferred tool for m-CAK stationary wind modelling.

The reader is referred to Araya et al. (2017) for the comprehensive set of equations described therein, as this work lays the essential groundwork and provides the mathematical formalism necessary to understand the advancements presented here. The primary goal of this work is to investigate and delineate the regions of co-existence of wind solutions subject to the combined effects of gravity darkening and stellar oblateness, thereby advancing our understanding of rotationally modified radiatively driven stellar winds. All solutions presented in this work correspond to the equatorial plane (co-latitude θ = 90°).

2.1 Radiation-driven wind equations

The 1D m-CAK hydrodynamical equations for rotating line-driven winds are formulated from the continuity and momentum equations, with the radiative acceleration parametrised following the standard CAK formalism (Castor et al. 1975; Friend & Abbott 1986). The basic equations are: $$4\pi r^2\rho \upsilon =\dot{M},$$(2)

and $$\upsilon \frac{d\upsilon }{dr}=-\frac{1}{\rho }\frac{dP}{dr}-\frac{GM\left(1-{\Gamma }_{\text{E}}\right)}{r^2}+\frac{{\upsilon }_{\varphi }^2\left(r\right)}{r}+g_{\text{rad}}^{\text{line}}\left(\rho ,d\upsilon /dr,n_{\text{E}}\right),$$(3)

where r is the radial coordinate, M the stellar mass, ρ the density, v the radial velocity, P is the gas pressure for an isothermal wind, and Γ_E the Eddington parameter due to electron scattering. For an ideal gas, the equation of state is given by P = a² ρ, being a the isothermal sound speed.

Following Curé (2004) and Araya et al. (2017), we adopt the usual variable transformations u = −R_*/r, w = υ/a, and w′ = dw/du, with a_rot = υ_rot/a, to write the compact form of the EoM as $$\begin{array}{c}F\left(u,w,w^{\prime }\right)=\left(1-\frac{1}{w^2}\right)w\frac{dw}{du}+A+\frac{2}{u}+a_{\text{rot}}^{2}u\\ -C^{\prime }{f}_{\text{OFD}}\left(\Omega ,\theta ,u\right)g\left(u\right)w^{-\delta }\left(w\frac{dw}{du}\right)=0,\end{array}$$(4)

where f_OFD is the oblate finite-disk correction factor including gravity darkening. Details of the derivation and definitions of the parameters are given in Curé (2004) and Araya et al. (2017). Here, α, k, and δ are the standard m-CAK line-force parameters (Abbott 1982; Puls et al. 2000). The parameter α describes the ratio between the line force from optically thick lines and the total line force. The parameter k (absorbed into the eigenvalue C′ in Eq. (4)) is related to the number of lines contributing effectively to the driving of the wind. Finally, δ accounts for changes in ionisation throughout the wind.

To solve the m-CAK equation of motion (Eq. (4)) while considering GD and the oblate distortion of a rapidly rotating star, we used the simplified approximation from Araya et al. (2017) for the f_OFD that makes it possible to solve the EoM in HYDWIND. This approximation preserves the known topology, and Araya et al. (2017) solved it in both equatorial and polar directions.

2.2 Model setup

Following Araya et al. (2017), we used the hydrodynamic code HYDWIND (Curé 2004) to solve the EoM. This code uses a β-law as initial trial solution, with β values in the range ≈0.8–1.2 for fast solutions and ≈3.0–3.5 for Ω-slow solutions (Araya et al. 2018), and terminal velocity speeds of 1000 km/s and 400 km/s, respectively. We begin our analysis by calculating the wind solution in the presence of the GD for the interval 0.68 < Ω < 0.99. Our calculations are typically carried out for main-sequence B-type stars, specifically B0 IV, B2.5 V, and B3 I. For the B2.5V star, we follow Madura et al. (2007), with stellar parameters set as: T_eff=20 000 K, R_*, = 4 R_⊙, and log 𝑔 = 4.11 dex. See Table 1, last row, for the base density value. The line force is parametrised following Gayley (1995), with the values α = 0.5, δ = 0.0, and k = 0.6098 or line strength parameter $\overline{Q}=1533$. For the B0 IV star, we use the parameters from Sigut & Jones (2007), while the line force parameters are taken from Pauldrach et al. (1986). The stellar and line-force parameters for the B3 I model are adopted from the equatorial bi-stable wind model developed by Pelupessy et al. (2000). The stellar and line force parameters of these three models are summarised in Table 1.

The assumption of invariant m-CAK parameters (α, k, and δ) with varying stellar rotation (Ω) is an inherent approximation. These parameters are fundamentally derived from atomic line transitions, which are sensitive to local physical conditions that can be altered by stellar rotation through oblateness and gravity darkening. There are a few works that address the calculation of the line-force parameters for the fast solutions (Pauldrach et al. 1986; Puls et al. 2000; Noebauer & Sim 2015; Gormaz-Matamala et al. 2019), but no studies have performed these calculations for the Ω-slow solution. In this work, following Stee & de Araujo (1994) and de Araujo (1995), we vary the value of the line force parameter α, which represents the ratio between the line-force originating from optically thick lines and the total line-force that controls the contribution of optically thick and thin lines.

3 Co-existence of fast and Ω-slow regimes

3.1 Numerical method for mapping solution regimes

The existence of the two distinct solution branches (fast and Ω-slow solutions) requires a robust numerical strategy to map each regime and accurately identify the interval of Ω where they may co-exist. To achieve this, we employed the following systematic procedure using the HYDWIND code.

Tracing the Ω-slow Solution Branch:

We begin at a near-critical rotational speed (e.g. Ω = 0.995), where only the Ω-slow solution is expected. At this starting point, a β-velocity law with a large exponent (typically β ≈ 3.0–3.5) is used as the initial trial solution. Once this solution converges, we gradually reduce Ω and, if necessary, fine-tune β. This process is continued until a value of Ω is reached, below which a converged Ω-slow solution can no longer be found.

Tracing the Fast Solution Branch:

Conversely, to map the fast solution branch, we start at a moderate rotational speed where the fast wind is dominant (e.g. Ω = 0. 68). A β-velocity law provides the initial guess with a small exponent (typically β ≈ 0.8). After convergence, we incrementally increase the value of Ω. This process is repeated until the fast solution no longer exists.

The co-existence interval is then defined as the range of Ω for which this procedure yields valid solutions for both the fast and Ω-slow branches. This systematic approach ensures that each solution branch is followed continuously, providing a robust determination of the boundaries for each regime.

3.2 GD effect

In this section, we focus on the co-existence region when only the GD effect is present. Our analysis is applied to the three stellar models (see Table 1). For each model, the base density, ρ₀, is kept constant, as indicated in the last row of this table. The calculations begin with Ω = 0.68, incrementally increasing to explore the fast wind solution. For the Ω-slow solution, the calculations utilise an initial trial velocity profile with β values ranging from 3 to 3.5. These are known to represent the characteristics of slower winds.

To compute the Ω-slow solutions, the process starts at a near-critical rotational speed (Ω = 0.995) and is gradually reduced. This ensures the solution follows the Ω-slow wind regime, as the velocity and density structures adapt to the lower rotation rate. This iterative approach allows for identifying the transition region where the fast solution ceases to exist and the Ω-slow solution emerges. Fig. 1 illustrates the mass loss and terminal velocity as a function of Ω, for all three stellar models. The different behaviour of fast and Ω-slow solutions is clearly seen. Moreover, co-existing regions where both types of solutions overlap are identified. As illustrated in Fig. 1 and detailed in Table 2, the co-existence interval between the fast and Ω-slow solutions varies for the three stellar models, highlighting the effects of stellar parameters and rotational speed. For the B2.5 V star, the co-existence interval is observed within the range 0.936 < Ω < 0.942, whereas for the B0 IV star, it occurs between 0.932 < Ω < 0.937.

In contrast, for the B3 I star, the interval significantly occurs on lower Ω, ranging from 0.899 < Ω < 0.91. The presence of GD causes the fast wind solutions to shift towards higher rotational velocities compared to the results of Araya et al. (2018), which did not account for GD effects. This shift in the fast solution also causes the co-existence interval to move towards higher rotational speeds. The inclusion of GD modifies the stellar surface flux distribution, which directly influences the effective gravity and subsequently alters the dynamics of the wind solutions. Our analysis also reveals a significant difference in the wind properties between the two solutions. Specifically, the terminal velocities of the fast solutions exceed those of the Ω-slow solutions by approximately a factor of ~2. In Araya et al. (2018), it was observed that for all three models, the mass loss rate remained nearly constant for the Ω-slow solutions. However, in the presence of GD, a different trend emerges in our calculations: the mass loss rate gradually increases as the rotational speed of the Ω-slow solutions increases.

This finding suggests that GD, which arises due to the influence of centrifugal forces causing the star’s equator to rotate faster than its poles, significantly affects its wind structure. As a result, the temperature distribution becomes dependent on co-latitude, leading to variations in the mass-loss rate of stellar winds. As rotation speed increases, the enhanced GD effect causes the star’s surface temperature to vary more significantly between the equator and poles. This temperature variation affects the atmospheric structure and consequently impacts the mass loss rate of the star. Therefore, the original models (Araya et al. 2018) without GD showed a plateau in mass loss rate for Ω-slow solutions, and the inclusion of GD effects resulted in a progressive increase in mass loss rate with increasing rotational speed.

Furthermore, we investigate the effect of the line-force parameter α on the co-existence region between the fast and Ω-slow solutions. Figures 2, 3, and 4 present the mass-loss rates and terminal velocities as a function of the rotational speed Ω for all stellar models. These solutions are computed with a value of k = 0.3 and by decreasing the value of a from 0.5 to 0.35. Our results show that as the value of α decreases, the co-existence region shifts towards the lower values of Ω. This shift indicates that for the higher values of a, higher rotational speeds are required for transition into the Ω-slow solution. This behaviour is consistent with the findings of Curé (2004); Araya et al. (2018), which state that for a weaker line force (lower a), lower rotational speeds are sufficient to achieve the Ω-slow solution. Since α is a line force parameter that influences the strength of the radiation-driven line force, a smaller value of α corresponds to a weaker line force, which requires less rotational speed to transition into the Ω-slow solution. Conversely, when α is larger, the line force becomes stronger, requiring higher rotational speeds to reach the Ω-slow regime. This shift in the co-existence region demonstrates the sensitivity of the wind structure to the value of a and how this parameter controls the rotational speeds at which the transition from the fast solution to the Ω-slow solution occurs.

Fig. 1 Mass-loss rate and terminal velocity as a function of the rotational speed. The upper panel shows the mass-loss rate, normalised to the non-rotating case. The lower panel shows the terminal velocities. These values are obtained using a stationary hydrodynamic solution with the HYDWIND code, accounting for the gravity darkening (GD) effect.

Fig. 2 Line-force parameter effects on wind solutions. The upper panel displays mass-loss rates normalised to the non-rotating case, and the lower panel shows terminal velocities. Solutions were derived from the B0 IV model in the presence of the GD, with systematically varied line-force parameter α at constant k = 0.3.

Fig. 3 Line-force parameter effects on wind solutions. Same as Fig. 2, but for the B2.5 V model.

Fig. 4 Line-force parameter effects on wind solutions. Same as Fig. 2, but for the B3 I model.

3.3 Oblateness effect

In this section, we focus on a single stellar model: a B0 IV star with the stellar and line-force parameters listed in Table 1. The calculations presented here consider only the effect of stellar oblateness.

As shown in Fig. 5, the fast wind solution exists in the interval 0.68 < Ω < 0.875, while the Ω-slow solution emerges in the range 0.876 < Ω < 0.95. Thus, when considering oblateness alone, no co-existence region is found between the fast and Ω-slow regimes. Compared to the case with gravity darkening, the transition between the two regimes occurs at lower rotational velocities.

We further explored the influence of varying line-force parameters. Lower values of α shift the transition to even lower rotational speeds, and Ω-slow solution emerges at even lower velocities. For instance, with α = 0.5 and k = 0.3, the fast solution is found within Ω < 0.81, and the Ω-slow solution within 0.82 < Ω, as shown in Fig. 5. Similarly, for α = 0.45, the fast solution spans Ω < 0.79, while the Ω-slow solution appears in the range 0.80 < Ω. In both cases, no co-existence interval is observed between the two solution regimes.

Fig. 5 Stellar oblateness effect on wind parameters. The upper panel shows the mass-loss rates, and the lower panel shows the terminal velocities as a function of the rotation parameter Ω. Calculations use the original B0 IV model parameters.

3.4 Both GD and oblate effects

In this section, we examine the combined effects of both GD and the oblate for the B0IV stellar model, using the stellar and line-force parameters provided in Table 1. The fast solution is found in the interval 0 < Ω < 0.923, while the Ω-slow solution is present in the range of 0.924 < Ω < 0.95 as shown in Fig. 6. This suggests that when both GD and oblateness are considered simultaneously, the solution shifts towards higher rotational velocities compared to the oblateness effect alone, but towards lower rotational velocities than when only gravity darkening is taken into account. There is no co-existence interval in this scenario.

Notably, the behaviour of wind parameters also differs from earlier results. When only GD or oblateness was applied, the mass loss rate increased steadily during the fast solution phase. When both effects are present, the mass loss rate initially increases during the fast solution phase, reaching a maximum, and then declines. This different behaviour is a direct result of the combined influence of GD and oblateness. When the Ω-slow solution is present, the mass loss rate slightly increases compared to the fast solutions, and follows the tendency to decrease as Ω increases.We were unable to find Ω-slow solutions for Ω values higher than 0.93, possibly because we could not find a suitable initial trial solution that our numerical code required (Curé 2004).

Fig. 6 Combined GD and oblateness effects in the B0 IV model. The upper panel shows the mass-loss rates, and the lower panel shows the terminal velocities as a function of the rotation parameter Ω. The inset panels provide a magnified view highlighting the region corresponding to the Ω-slow wind solutions.

4 Discussion and conclusion

This work aims to identify the regimes in which a star’s wind transitions from the fast solution to the Ω-slow solution and to determine the co-existence interval where both solutions can co-exist. These wind solutions are obtained by solving 1D stationary hydrodynamic equations at the stellar equatorial plane. By examining the impact of various stellar parameters, including rotational speed and line-force parameters, we aim to gain a deeper understanding of the conditions under which these distinct 1D wind solutions are present.

As a necessary clarification, the core goal of this work is to delineate the regimes and co-existence interval of the Ω-slow solution, which appears when the rotational speed exceeds ≈ 75% of the critical speed (Curé 2004; Araya et al. 2018; Curé & Araya 2023). To achieve this, we utilised a 1D equatorial model, a methodological approximation that omits crucial multidimensional phenomena. Previous multi-dimensional studies, which focused primarily on the standard (fast) solutions and the Wind Compressed Disk (WCD) model (Bjorkman & Cassinelli 1993; Cranmer & Owocki 1995; Owocki et al. 1996), consistently demonstrated that non-radial line-force components (NRFs) inhibit WCD formation and induce a prolate (pole-enhanced) wind morphology, often contrasting with the predictions of equatorial compression inherent to 1D models. Crucially, these findings regarding the dominance of NRFs in shaping the prolate wind structure were derived without consistently considering the specific physics (or topology) of the Ω-slow solution, as the existence of this solution was established later (Curé 2004; Madura et al. 2007; Araya et al. 2018). Therefore, while our results map the stability and existence of the Ω-slow regime in 1D, the eventual global morphology (oblate or prolate) of this denser, slower outflow remains an open question that requires future comprehensive 2D/3D modelling, including NRFs, GD and the oblate shape of the star, within the Ω-slow regime. Therefore, while the 1D approximation constrains our results to the equatorial plane, it provides a physically consistent framework to isolate and characterise the onset and stability of the Ω-slow regime before extending to multidimensional treatments.

According to Araya et al. (2018), the variability observed in the line profiles of massive stars can be attributed to disturbances in the photosphere, such as non-radial pulsations, stellar spots, and other perturbations. These factors, individually or in combination, can significantly alter the photospheric density when the star’s rotational speed is within the co-existence range. Such density changes can potentially trigger a transition in the wind regime, particularly in the equatorial region. When the rotational velocity is within this critical range, the star may transit from a fast wind solution to an Ω-slow wind solution. This transition is thought to play a key role in creating a circumstellar disc, which can form around the star’s equatorial plane under these conditions. In this Ω-slow wind regime, radiation pressure is the primary mechanism for mass injection into the disc. However, if the disturbances within the star’s structure are not sufficiently high, the system may remain stable within the co-existence region, with the wind solution maintaining its current state. However, if the amplitude of these disturbances exceeds a certain threshold, the system might undergo another transition, reverting to the fast wind regime. This would result in the dissipation of the equatorial disc, implying that the existence and stability of these discs are closely linked to the nature of the stellar wind solution and the magnitude of the perturbations.

We analysed three different stellar models under three distinct scenarios: the presence of only GD, only the oblateness effect, and the combined effect of both GD and oblateness. Our results are compared to those obtained by Araya et al. (2018), who studied wind solutions without including these effects. We find that when GD and oblateness are included, the wind solutions exhibit a shift towards higher rotational velocities compared to the results of Araya et al. (2018). This suggests that both GD and oblateness play a substantial role in altering the wind dynamics. Furthermore, the co-existence interval, where both fast and Ω-slow solutions exist simultaneously, also shifts to higher rotational velocities under these conditions. Notably, the shift in wind solutions is most pronounced where only the GD effect is considered. The combined impact of GD and oblateness also leads to significant changes, but the GD effect alone induces the most noticeable shift towards higher rotational speeds.

We also investigated the effect of varying the line-force parameter α on the wind solution. Our analysis reveals that the location and width of the co-existence region depend on stellar and line-force parameters as the value of α increases, and the wind solution shifts towards higher rotational velocities. We also find that during the Ω-slow solution, the mass loss rate gradually increases in the presence of these effects. However, when accounting for oblateness and gravity darkening, the wind’s mass-loss rate is lower compared to cases without gravity darkening. As mentioned by Aufdenberg et al. (2006), the mass-loss rate is primarily influenced by the equatorial mass density, which significantly affects the strength of emission profiles. For a given set of line-force parameters, increasing the rotational velocity Ω generally leads to a higher emission line intensity, whether the star is assumed to be spherical or oblate with a constant temperature. However, in the case of an oblate star with gravity darkening, the maximum line intensity is reached at an Ω value less than 1. As Ω approaches 1, the emission decreases due to the reduction in effective temperature caused by gravity darkening. Thus, the mass-loss rate decreases.

In the context of the m-CAK theory, the presence of fast and Ω solutions, along with their co-existing interval, provides insight into stellar wind dynamics. Including GD and oblateness significantly impacts these solutions, particularly at lower line force parameter α values. GD reduces the radiative flux in the equatorial plane due to weaker surface gravity, while oblateness alters the stellar geometry, further weakening the driving force for winds. At lower α, the radiative acceleration becomes insufficient to sustain outflows, especially near the equator, where the velocity drops below the threshold required for the material to escape. The absence of wind solutions for lower α may result from the shift in critical points within the wind equations due to physical limitations or numerical challenges in resolving the topology of the wind equations under these conditions. The terminal velocity of the Ω-slow solutions is ≈200 km/s; thus, decreasing α leads to reduced equatorial velocities, ultimately inhibiting outflows in these regions.

In our work, we do not include the effect of viscosity in the wind equation of motion. Thus, angular momentum is conserved locally without viscosity, and the material follows purely radial outflow. In contrast, viscous forces enable the outward transfer of angular momentum, thereby reducing the terminal velocity of the wind. This is a key mechanism in the viscous decretion disk model (Lee et al. 1991; Okazaki 2000; Curé et al. 2022), which successfully explains the long-lived disks observed around Be stars. Thus, the inclusion of viscosity accounts for the internal angular momentum transport or the redistribution of material due to viscous forces, which can significantly influence the wind structure and disk formation. Considering this effect is beyond the scope of the present work. However, in future studies, we intend to investigate the impact of viscosity on wind solutions to gain a deeper understanding of how material flows in the equatorial region.

Acknowledgements

We thank the anonymous referee for the valuable comments and suggestions, which helped to improve the clarity and quality of this manuscript. MC, IA, and CA also thank the ANID FONDECYT project 1230131. MC and CA acknowledge support from Centro de Astrofísica de Valparaiso – CAV, CIDI N. 21 (Universidad de Valparaiso, Chile). AC acknowledges partial support from Centro de Estudios Atmosféricos y Cambio Climático, CIDI CR 212.115.013. AC thanks the support from Proyecto PUENTE, financiado por el proyecto PFE UVA20993, Fortalecimiento del Sistema de Investigación e Innovación de la Universidad de Valparaiso, Universidad de Valparaiso. This project was funded by the European Union (Project 101183150 - OCEANS). Powered@NLHPC: This research was partially supported by the supercomputing infrastructure of the NLHPC (ECM-02).

References

All Tables

All Figures

	Fig. 1 Mass-loss rate and terminal velocity as a function of the rotational speed. The upper panel shows the mass-loss rate, normalised to the non-rotating case. The lower panel shows the terminal velocities. These values are obtained using a stationary hydrodynamic solution with the HYDWIND code, accounting for the gravity darkening (GD) effect.
In the text

	Fig. 2 Line-force parameter effects on wind solutions. The upper panel displays mass-loss rates normalised to the non-rotating case, and the lower panel shows terminal velocities. Solutions were derived from the B0 IV model in the presence of the GD, with systematically varied line-force parameter α at constant k = 0.3.
In the text

	Fig. 3 Line-force parameter effects on wind solutions. Same as Fig. 2, but for the B2.5 V model.
In the text

	Fig. 4 Line-force parameter effects on wind solutions. Same as Fig. 2, but for the B3 I model.
In the text

	Fig. 5 Stellar oblateness effect on wind parameters. The upper panel shows the mass-loss rates, and the lower panel shows the terminal velocities as a function of the rotation parameter Ω. Calculations use the original B0 IV model parameters.
In the text

	Fig. 6 Combined GD and oblateness effects in the B0 IV model. The upper panel shows the mass-loss rates, and the lower panel shows the terminal velocities as a function of the rotation parameter Ω. The inset panels provide a magnified view highlighting the region corresponding to the Ω-slow wind solutions.
In the text