Co-existence and switching between fast and Ω-slow wind solutions in rapidly rotating massive stars

Abstract

Most radiation-driven winds of massive stars can be modelled with m-CAK theory, resulting in the so-called fast solution. However, the most rapidly rotating stars among them, especially when the rotational speed is higher than |${\sim } 75\hbox{ per cent}$| of the critical rotational speed, can adopt a different solution, the so-called Ω-slow solution, characterized by a dense and slow wind. Here, we study the transition region of the solutions where the fast solution changes to the Ω-slow solution. Using both time-steady and time-dependent numerical codes, we study this transition region for various equatorial models of B-type stars. In all cases, in a certain range of rotational speeds we find a region where the fast and the Ω-slow solution can co-exist. We find that the type of solution obtained in this co-existence region depends stongly on the initial conditions of our models. We also test the stability of the solutions within the co-existence region by performing base-density perturbations in the wind. We find that under certain conditions, the fast solution can switch to the Ω-slow solution, or vice versa. Such solution-switching may be a possible contributor of material injected into the circumstellar environment of Be stars, without requiring rotational speeds near critical values.

hydrodynamics, stars: early-type, stars: rotation, stars: winds, outflows

1 INTRODUCTION

Massive stars are characterized by strong outflowing stellar winds, driven by the transfer of momentum from the radiation field to the plasma via scattering processes in spectral lines. These winds are called radiation-driven winds and are described by the standard radiation-driven wind theory (or CAK theory), based on the pioneering work of Castor, Abbott & Klein (1975). Modifications to the CAK model (m-CAK) were carried out by Friend & Abbott (1986) and Pauldrach, Puls & Kudritzki (1986), who implemented the finite disc correction factor, which takes account of the cone angle subtended by the star instead of using a point star approximation. This m-CAK model has been successful in describing both terminal wind velocities (⁠|$v$|_∞) and mass-loss rates (⁠|$\dot{M}$|⁠) from hot stars (Puls, Vink & Najarro 2008).

The effect of stellar rotation on wind dynamics was also investigated by Friend & Abbott (1986) and Pauldrach et al. (1986). Both studies found that in general stellar rotation increases the mass-loss rate and decreases the terminal speed of the wind for slow or moderate rotational speeds. In particular, Friend & Abbott (1986) were unable to find a wind solution for rotational speeds higher than about |$75\hbox{ per cent}$| of the critical rotational speed. This was a strong indication that something unusual happens to the wind solution for rapidly rotating massive stars.

Among massive stars there are many rapid rotators, including some with rotational speeds near the critical rate. This latter condition has been generally assumed to explain the formation of circumstellar discs, for example in Be stars (Townsend, Owocki & Howarth 2004; Rivinius, Carciofi & Martayan 2013). A 2D model to link radiation-driven wind theory to the formation of a dense thin equatorial disc in Be stars was developed by Bjorkman & Cassinelli (1993) and termed the wind compressed disc (WCD) model. Then Owocki, Cranmer & Blondin (1994) implemented a time-dependent WCD model with purely radial forces and confirmed the basic WCD effect. However, Cranmer & Owocki (1995) included non-radial line forces in this model and showed that these forces inhibit the formation of an equatorial disc. This result was confirmed by Petrenz & Puls (2000).

Curé (2004) revisited the effects of stellar rotation in the theory of radiation-driven winds for 1D stationary models. He found a new type of wind solution when the rotational speed (⁠|$v$|_rot) is above |${\sim } 75\hbox{ per cent}$| of the critical rotational speed (⁠|$v$|_crit) and analysed its application to explain the equatorial disc of Be stars. This new solution, called the Ω-slow solution, where Ω = |$v$|_rot/|$v$|_crit, is slower and denser than the standard solution (hereafter the fast solution). This model combined with the bi-stability effect (Lamers & Pauldrach 1991) predicts equator-to-pole density contrasts of up to 10⁴, suggesting that the Ω-slow solution could play a key role in the formation of dense equatorial discs in B[e] stars (Curé, Rial & Cidale 2005). With the purpose of gaining a better understanding of this slow wind solution and its stability at high rotational speeds, Madura, Owocki & Feldmeier (2007) calculated time-dependent simulations and concluded that it was unlikely that the Ω-slow solution could reproduce the inferred densities of B[e]discs. In addition, they reported the presence of abrupt kink transitions in the wind velocity profiles at rotation rates of about 75–|$85\hbox{ per cent}$| of the critical rotational speed.

Despite the importance of 2D simulations, 1D calculations represent the best-case scenario. Thus, this paper aims to study the region where the fast solution switches to the Ω-slow solution (or vice versa) for a 1D equatorial model with a high rotational speed. Our calculations are performed using stationary and time-dependent hydrodynamic codes. We also study the dependence of a co-existence region on the stellar and line-force parameters (m-CAK formalism), and the stability of the solutions obtained within this region, using base-density wind perturbations.

The paper is organized 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 addition, a co-existence region is found where both types of solutions are present. The influence of the line-force parameters on this regions is also studied. Section 4 presents the solutions obtained from time-dependent simulations and a comparison between stationary and time-dependant results. In Section 5, we perform wind base-density perturbations, with the purpose of testing: (i) the stability of the two solutions inside the co-existence region, and (ii) the minimum conditions to induce regime-switching between wind solutions. Section 6 discusses the application of our results. Finally, we present our conclusions in Section 7.

2 HYDRODYNAMIC WIND MODEL

We adopt the time-dependent 1D rotating radiation-driven m-CAK model for the equator of a massive star. This model uses the following approximations: (i) an isothermal wind with temperature equal to the stellar effective temperature, and (ii) angular momentum conservation. This model neglects the effects of viscosity, heat conduction and magnetic fields. In addition, the distortion of the star's shape caused by its high rotational speed, gravity darkening, non-radial velocities and multiscattering are not considered.

2.1 Time-dependent equations

At the equatorial plane, the relevant 1D time-dependent hydrodynamic equations, namely continuity and radial momentum, are

\begin{eqnarray} \frac{\mathrm{\partial} \rho }{\mathrm{\partial} t} + \frac{1}{r^{2}}\frac{\mathrm{\partial} }{\mathrm{\partial} r}(\rho \, r^{2}\,v)=0 \end{eqnarray}

(1)

and

\begin{eqnarray} \frac{\mathrm{\partial} v}{\mathrm{\partial} t} + v\frac{\mathrm{\partial} v}{\mathrm{\partial} r}=-\frac{1}{\rho }\frac{\mathrm{\partial} p}{\mathrm{\partial} r}- \frac{G\,M (1-\Gamma _{\mathrm{E}})}{r^{2}}+\frac{v^{2}_{\phi }}{r} + g^{\mathrm{line}}, \end{eqnarray}

(2)

where t is the time, r is the radial coordinate, |$v$| is the fluid radial velocity, ρ is the mass density and g^line is the radiative line acceleration described in more detail below. The ideal gas pressure (p = a² ρ) is given in terms of the isothermal sound speed, a. The azimuthal speed |$v$|_ϕ = |$v$|_rot R_*/r is calculated assuming conservation of angular momentum, where |$v$|_rot is the star's surface rotational speed at the equator and R_* is the stellar radius. The classical Eddington factor, expressed as

\begin{eqnarray} \Gamma _{\mathrm{E}} = \frac{\sigma _{\mathrm{E}}\,L}{4\, \pi \,c \,G \,M}, \end{eqnarray}

(3)

is the ratio between the Thomson electron-scattering force and the gravitational force. The electron-scattering opacity per unit mass σ_E is given by

\begin{eqnarray} \sigma _{\mathrm{E}} = \sigma _{\mathrm{Th}} \frac{n_{\mathrm{E}}}{\rho }, \end{eqnarray}

(4)

with σ_Th the Thomson cross-section of electrons and n_E the electron number density. In addition,

\begin{eqnarray} \frac{n_{\mathrm{E}}}{\rho } = \, D\, = \frac{1}{m_{\mathrm{p}}} \left( \frac{1+Z_{\mathrm{He}}\, Y_{\mathrm{He}}}{1+4\, Y_{\mathrm{He}}} \right), \end{eqnarray}

(5)

where m_p is the mass of a proton, Y_He = n_He/n_H is the helium abundance relative to hydrogen, and Z_He is the number of free electrons provided by He.

Following the descriptions of Abbott (1982), Friend & Abbott (1986) and Pauldrach et al. (1986), the m-CAK standard parametrization for the radiative line acceleration term is

\begin{eqnarray} g^{\mathrm{line}} &=& \frac{\Gamma _{\mathrm{E}}\,G\,M\,k}{r^{2}} \left( \frac{1}{\sigma _{\mathrm{E}}\,v_{\mathrm{th}}} \right)^{\alpha } \left( \frac{1}{\rho } \frac{\mathrm{\partial} v}{\mathrm{\partial} r} \right)^{\alpha } \left( \frac{n_{\mathrm{E11}}}{W(r)} \right)^{\delta }\nonumber\\ &&\times\, f_{\mathrm{FD}}\left( r,v,\frac{\mathrm{\partial} v}{\mathrm{\partial} r}\right), \end{eqnarray}

(6)

where |$v$|_th is the mean thermal velocity of the protons, W(r) is the dilution factor and n_E11 is the electron number density in units of 10⁻¹¹ cm⁻³.

The m-CAK line-force parameters (Abbott 1982; Puls, Springmann & Lennon 2000) are α (the ratio between the line-force from optically thick lines and the total line-force), k (which is related to the number of lines contributing effectively to the driving of the wind) and δ (which accounts for changes in the ionization throughout the wind).

We include the effects of a finite-disc correction factor, f_FD, resulting from the finite size of the star (see equation 50 in Castor et al. 1975).

A different approach for the radiative line acceleration was considered by Gayley (1995), who introduced an alternative formalism whereby the line-force parameter k is replaced by the dimensionless line-strength parameter |$\bar{Q}$| and the cutoff parameter Q_o. These parameters are related by the following identity:

\begin{eqnarray} k= \frac{1}{1-\alpha } \left( \frac{v_{\mathrm{th}}}{c} \right)^{\alpha } \, \bar{Q}\, Q_{\mathrm{o}}^{-\alpha }. \end{eqnarray}

(7)

Based on the formalism of Gayley (1995), the g^line term can be expressed as

\begin{eqnarray} g^{\mathrm{line}} &=& \frac{\Gamma _{\mathrm{E}}\,G\,M\,\bar{Q}}{\left( 1 - \alpha \right) r^{2}} \left( \frac{1}{Q_{\mathrm{o}}\,\sigma _{\mathrm{E}}\,c} \right)^{\alpha } \left( \frac{1}{\rho } \frac{\mathrm{\partial} v}{\mathrm{\partial} r} \right)^{\alpha } \left( \frac{n_{\mathrm{E11}}}{W(r)} \right)^{\delta }\nonumber\\ &&\times f_{\mathrm{FD}}\left( r,v,\frac{\mathrm{\partial} v}{\mathrm{\partial} r}\right). \end{eqnarray}

(8)

Here we follow the original CAK formalism, but conversion to Gayley's method can be readily made. Owocki, Castor & Rybicki (1988) use a slightly different formalism by introducing a maximum opacity κ_max that limits the effects of strong driving lines in their line-deshadowing instability models.

2.2 Steady-state equations

For our stationary 1D spherical flow, we combine the continuity equation (1) and momentum equation (2) assuming |$\mathrm{\partial} /\mathrm{\partial} t \, \rightarrow 0$|⁠. This leads to our basic equation of motion:

\begin{eqnarray} \left( 1- \frac{a^{2}}{v^{2}} \right) v \frac{{\rm d}v}{{\rm d}r} =\frac{2\, a^{2}}{r} - \frac{G\,M (1-\Gamma _{\mathrm{E}})}{r^{2}}+\frac{v^{2}_{\phi }}{r} + g^{\mathrm{line}}. \end{eqnarray}

(9)

In order to facilitate a solution, the change of variables u = −R_*/r, |$w$| = |$v$|/a and w΄ = dw/du leads to (Castor et al. 1975; Curé 2004)

\begin{eqnarray} \nonumber F(u,w,w^{\prime }) & \equiv & \left( 1- \frac{1}{w^{2}} \right) w \frac{{\rm d}w}{{\rm d}u} + A + \frac{2}{u} + a^{2}_{\mathrm{rot}}\, u \\ & & - \, C^{\prime } \, f_{\mathrm{FD}} \, g(u)\, (w)^{-\delta } \left( w \frac{{\rm d}w}{{\rm d}u} \right)^{\alpha } =0, \end{eqnarray}

(10)

where a_rot = |$v$|_rot/a, |$A= \frac{G\, M \left[ 1- \Gamma _{\mathrm{E}}\right]}{a^{2}R_{*}}$|⁠, |$g(u)= \big( \frac{u^{2}}{1-\sqrt{1-u^{2}}} \big)^{\delta }$|

and

\begin{eqnarray} C^{\prime } = C \left( \frac{\dot{M}\, D}{2\, \pi } \frac{10^{-11}}{a\,R_{*}^{2}} \right)^{\delta } \left( a^{2}R_{*} \right)^{\alpha -1}. \end{eqnarray}

(11)

Here C΄ is the eigenvalue of the problem where |$C= \Gamma _{\mathrm{E}}\,G \,M \,k \left( {4\pi }/{\sigma _{\mathrm{E}}\, v_{\mathrm{th}}\, \dot{M}}\right)^{\alpha }$| and D is given by equation (5).

Based on this formalism, Curé (2004) demonstrated the existence of another physical solution for fast-rotating stars.

3 1D STEADY-STATE WIND SOLUTIONS FOR FAST-ROTATING STARS

When solving equation (10) for high rotational speeds, Curé (2004) found that the standard (or fast) solution ceases to exist and a new kind of solution is found from Ω ≳ 0.75, the so-called Ω-slow solution. Now we want to explore if there is an interval in the Ω-parameter space where the wind solution switches from the fast to the Ω-slow solution or vice versa. If this interval exists, it might be possible that the two solutions can co-exist.

In order to obtain numerical wind solutions, we used our stationary hydrodynamic code Hydwind described in Curé (2004). Hydwind solves equation (10) by the relaxation method, and thus requires an initial trial solution for the velocity profile and an estimate for the eigenvalue to begin the integration. Usually a β-law is used as the initial trial solution. For fast solutions, we typically use β ≈ 0.8 and |$v$|_∞ ≈ 1000 km s⁻¹ and for slow solutions β ≈ 3.5 and |$v$|_∞ ≈ 400 km s⁻¹. The type of solution obtained (fast or Ω-slow) depends solely on the initial trial solution for the velocity profile and not on the base wind density, because the latter is a fixed constant independent of the rotational speed.

We begin our analysis by calculating wind solutions for 0 ≤ Ω ≤ 0.95 using the stellar parameters for a typical main-sequence B2.5 star, following the work of Madura et al. (2007): T_eff =20 kK, R_* = 4 R_⊙ and logg = 4.11. The base density was fixed at ρ₀ = 8.709 × 10⁻¹³ g cm⁻³ and the line-force was parametrized using the formalism of Gayley (1995): α = 0.5, δ = 0.0 and |$\bar{Q}=1533$|¹ (or k = 0.6098, according to equation 7).

Fig. 1 shows some of our calculated stationary velocity and density profiles. In this figure, it can clearly be seen that the fast solution regime is characterized by high velocities and low densities, while the Ω-slow solution regime is characterized by lower velocities and higher densities. The switch of regime can be observed at Ω ∼ 0.74.

Figure 1.

Steady-state wind solutions at various rotational speeds Ω for a main-sequence B2.5 stellar model. The solutions were obtained with the stationary code Hydwind. Top and bottom panels correspond to the velocity and density profiles as a function of the inverse radial coordinate u, respectively. The solid and dashed lines correspond to fast and Ω-slow solutions, respectively.

In Fig. 2 we show the mass-loss rate (in units of the mass loss of the non-rotating case), the ratio between the terminal and the escape velocity, and the location of the critical point as a function of the rotational speed (the equation of motion topology and the integration procedure are explained in detail in Curé 2004). The differences between fast and Ω-slow solutions are quite clear. For Ω ∼ 0.74 there is a jump in the terminal velocity and the position of the critical point between these solutions, while the mass-loss rate displays a change in slope.

Figure 2.

Mass-loss rate (in units of the mass loss for the non-rotating case, top panel), terminal velocity (in terms of the escape velocity, middle panel) and location of the critical point (bottom panel) as a function of the rotational speed. These values were obtained using a hydrodynamical stationary solution with the code Hydwind using |$\dot{M}_{\Omega =0} =1.004\times 10^{-9}\, {\rm M}_{\odot }\mathrm{yr}^{-1}$|⁠.

The velocity profiles shown in Fig. 1 as well as the values of the mass-loss rates and |$v$|_∞/|$v$|_esc shown in Fig. 2 are similar to the analytical results given by Madura et al. (2007, see their figs 3 and 4). In that work, the authors described these solutions as a flow with a steep and a shallow acceleration for Ω ≤ 0.74 and Ω ≥ 0.75, respectively. In Section 4 we will discuss our time-dependent calculations and compare them with the results obtained by Madura et al. (2007).

3.1 Co-existence of fast and slow regimes

In this section, we focus on the interval in the Ω-space parameter domain where the fast solution ceases to exist and the Ω-slow solution emerges.

We apply our analysis to three models taken from the literature: a main-sequence B2.5 star (as in the previous section), a subgiant B0 star and a supergiant B3 star. The stellar and line-force parameters for these objects are given in Table 1. The B0 IV model parameters are from Sigut & Jones (2007) and the line-force parameters are derived from Pauldrach et al. (1986). The stellar and line-force parameters for the B3 I model are taken from the equatorial bi-stable wind model in Pelupessy, Lamers & Vink (2000).

Table 1.

Stellar and line-force parameters for our models.

Model	Teff	log g	R*	α	δ	k	|$\bar{Q}$|	ρ0
	[kK]	[dex]	[R⊙]					[g cm−3]
B0 IV	25.0	3.50	10.00	0.565	0.02	0.32	2792.0	5.0 × 10−11
B2.5 V	20.0	4.11	4.00	0.500	0.00	0.61	1533.0	8.7 × 10−13
B3 I	17.5	2.70	47.00	0.450	0.00	0.57	361.6	1.0 × 10−11

The wind solutions of these three models are calculated with Hydwind for various values of Ω. Their base densities are fixed to a certain value ρ₀, which has the same value in each model for both types of solution. These values are given in the last column of Table 1. The calculations begin with a model with Ω = 0 and, then, the rotational speed is gradually increased until the last fast solution is attained. For Ω-slow solutions, a β-law with β ≈ 3.5 and |$v$|_∞ ≈ 400 km s⁻¹ is used as the initial trial solution. In order to obtain the Ω-slow solutions, we start with a model rotating at Ω = 0.99, and then Ω is gradually decreased until the last Ω-slow solution is obtained. The calculated mass-loss rates and terminal velocities as a function of Ω are shown in Fig. 3. We can observe clearly the characteristic behaviour of the fast and Ω-slow solutions. Furthermore, it is found that there exists an interval (or region) where fast and Ω-slow solutions co-exist. The co-existence region of our models is enlarged and highlighted in Fig. 4. This co-existence region is different for each model and is located roughly in the interval 0.65 ≲ Ω ≲ 0.76. In Appendix A, we show our full topological analysis when Ω = 0.73. Table 2 summarizes the properties of the various solutions and lists the wind parameters calculated at the minimum (Ω₁) and maximum (Ω₂) values encompassed by the co-existence region. Overall, the terminal velocities of fast solutions are higher than those from the Ω-slow solutions by a factor of about 2. On the other hand, the mass-loss rates computed for the fast and Ω-slow solutions at Ω₁ are the same. However, the maximum mass-loss rate of the fast solution, at the upper limit of the interval, |$\dot{M}(\Omega _2$|⁠), is not more than 7 per cent of the value calculated for the Ω-slow solution.

Figure 3.

Mass-loss rates in units of the non-rotating case (top panel) and terminal velocities (bottom panel) as a function of Ω for the solutions obtained from three stellar models. The solutions are calculated with Hydwind. An overlap of fast and Ω-slow solutions is observed in the range 0.65 ≲ Ω ≲ 0.75.

Figure 4.

Similar to Fig. 3, but with a zoomed-in focus on the region where the two types of solutions co-exist. This region is highlighted for each stellar model.

Table 2.

Properties of the fast and Ω-slow solutions for our three stellar models. Wind parameters are calculated at values of Ω where the co-existence region begins and ends.

Model	Ω range	Solution type	Co-existence interval	|$\dot{M}(\Omega _1$|⁠)	|$v$|∞(Ω1)	rcrit(Ω1)	|$\dot{M}(\Omega _2$|⁠)	|$v$|∞(Ω2)	rcrit(Ω2)
			Ω1 − Ω2	[M⊙ yr−1]	[km s−1]	[R*]	[M⊙ yr−1]	[km s−1]	[R*]
B0 IV	0 − 0.760	Fast	0.749 − 0.760	2.44 × 10−7	880.4	1.056	2.52 × 10−7	844.5	1.058
	0.749 − 1	Ω-slow		2.44 × 10−7	426.7	16.32	2.45 × 10−7	422.3	16.62
B2.5 V	0 − 0.736	Fast	0.723 − 0.736	2.16 × 10−9	1109.3	1.027	2.26 × 10−9	1061.7	1.028
	0.723 − 1	Ω-slow		2.15 × 10−9	485.1	22.36	2.15 × 10−9	478.7	22.36
B3 I	0 − 0.687	Fast	0.666 − 0.687	1.38 × 10−6	616.1	1.048	1.47 × 10−6	578.9	1.050
	0.666 − 1	Ω-slow		1.38 × 10−6	305.1	13.83	1.38 × 10−6	299.0	14.05

Finally, we study the influence of the line-force parameters on the co-existence region. Figs 5 and 6 display mass-loss rates and terminal velocities as a function of Ω, as derived from the solutions obtained with Hydwind by systematically varying the values of α and δ for the stellar model B2.5 V. As the value of α increases there is a ‘shift’ of this region towards higher values of Ω. This is similar to the results obtained by Curé (2004), where for thin lines driving the wind (lower α), lower rotational speeds are needed to obtain a Ω-slow solution. Furthermore, the width of this region seems to be independent of the value of α. Regarding δ, the width of the co-existence region increases as δ increases below a certain threshold, in this case δ ≤ 0.04. For values larger than the threshold, a gap is observed (Venero et al. 2016), confirming the absence of monotonically increasing stationary solutions.

Figure 5.

Mass-loss rates in units of the non-rotating case (top panel) and terminal velocities (bottom panel) as a function of Ω for solutions obtained from the B2.5 V model, varying systematically the value of the line-force parameter α. A shift of the co-existence region towards larger rotational speeds can be seen as α increases.

Figure 6.

Similar to Fig. 5, but with the models calculated varying the value of the line-force parameter δ. A gap between fast and Ω-slow solutions can be seen for the highest values of δ.

4 1D TIME-DEPENDENT WIND SOLUTIONS FOR ROTATING STARS

In order to study the stability of the steady-state solutions presented in the previous section, we solve here the time-dependent m-CAK model in the equatorial plane. For this purpose, we employ the hydrodynamic code Zeus-3D (Clarke 1996). We have adapted this Eulerian finite-difference code to our specific problem, including the radiative driving terms, along with all our assumptions (see Section 2).

Our spatial mesh has 500 zones, from r₁ = R_* to r₅₀₀ = 100 R_*, distributed non-uniformly with higher concentrations near the stellar surface to account for steeper flow gradients. The zone spacing is defined by Δr_i + 1/Δr_i = 1.02, with Δr_i = r_i + 1 − r_i. Our test calculations have shown that this grid resolution is sufficient for our purpose and that a grid with higher resolution will produce very similar results.

At the inner boundary we set the base density to a factor of the sonic point density, ρ(R_*) = f ρ(R_s), in order to obtain a stable outflow. We found that f in the range 5 < f < 20 gives a stable solution. In addition, the inflow inner velocity boundary can be adjusted by linear extrapolation from the closest zones in the interior of the computational domain. This boundary also is limited to subsonic flow. At the outer boundary, we set the condition of an outflowing wind.

In order to have numerical stability in the integration process, we followed Madura et al. (2007), who truncate the radial velocity gradient to zero whenever it is negative; that is, d|$v$|/dr → max (d|$v$|/dr, 0). This comes from a simple compromise between the likely extremes of d|$v$|/dr: (i) a lower limit where a negative velocity gradient implies a prior line resonance that shadows radial photons from the star, and (ii) an upper limit where this shadowing can be weakened by forward scattering (for further details see Madura et al. 2007).

In order to begin the integration of the time-dependent hydrodynamic equations for rotating winds with Zeus-3D, we started with the stationary solutions obtained by Hydwind as initial or trial solutions. Moreover, we also performed calculations using various β-velocity field trial solutions, as we did for the stationary cases in the previous section.

Fig. 7 shows the velocity profiles for various values of Ω obtained with Zeus-3D for the model B2.5 V (see Table 1); we have overplotted the stationary solutions shown in Fig. 1. The agreement between the two sets of results, which used completely different numerical algorithms, is excellent. The relative error with respect to the stationary solution for terminal velocities is about |$1\hbox{ per cent}$|⁠, while the relative error in the mass-loss rate is less than |$12\hbox{ per cent}$|⁠. The disagreement in the mass-loss rate is caused mainly by the difference between the density profiles. This variance arises from the use of two codes (time-dependent and stationary), which employ different numerical approaches for different equations of motion.

Figure 7.

Velocity profiles obtained with the code Zeus-3D for various values of Ω for the B2.5 V model. As initial condition, a stationary solution for each value of Ω is used. Solid and dashed lines represent the fast and slow solutions, respectively. The time-dependent solutions are compared with the stationary solutions (black solid lines overlaid by the time-dependent solutions).

It is important to stress here that our results show that there are no kinks in the resulting velocity profiles calculated with time-dependent solutions (see Fig. 7), in contrast to the results obtained by Madura et al. (2007).

In order to reproduce the Madura et al. (2007) results, we calculated time-dependent solutions employing first an initial non-rotating trial solution as these authors had done. However, in our case we started from the solutions obtained by Hydwind. These results are shown in the top panel of Fig. 8, where kinks in the velocity profiles are now present for Ω = 0.74, 0.76, 0.80 and 0.82, in agreement with the calculations performed by Madura et al. (2007). These kinks occur far from the stellar surface owing to the overloading of the wind and become more prominent as the value of Ω increases. It is worth noting that these kink-solutions should be the Ω-slow solutions shown in Fig. 7 with dashed lines.

Figure 8.

Velocity profiles obtained with the code Zeus-3D at various values of Ω for the B2.5 V model. Top panel: As the initial condition we used a stationary solution with Ω = 0 from Hydwind. We obtained fast and Ω-slow solutions but, now, models with Ω = 0.74, 0.76, 0.80 and 0.82 present a kink, similar to the velocity profiles obtained by Madura et al. (2007). Bottom panel: As the initial condition we adopted a stationary solution with Ω = 0.95 from Hydwind. The Ω-slow solutions are recovered but all the fast solutions present a kink.

Second, we repeated these calculations, but now using as the initial trial solution a stationary solution with Ω = 0.95 from Hydwind. The bottom panel of Fig. 8 shows the different velocity profiles. Here it can be seen that the velocity profiles from the Ω-slow solutions do not present any kink, whereas the velocity profiles of all fast solutions do. These kinks follow the same behaviour as in our previous calculations; that is, the lower the rotational speed, the farther the kink from the stellar surface.

We note that the stationary solutions with a kink (attained with Zeus-3D) use the truncation condition d|$v$|/dr → max (d|$v$|/dr, 0) from the kink location up to the outer boundary. In contrast, models without a kink have a positive velocity gradient along the whole velocity profile so the truncation condition is not used.

We conclude that kinks arise when an initial (trial) solution is used that belongs to a different solution branch in the topology of the m-CAK model (see the topological analysis by Curé & Rial 2007).

4.1 Co-existence of fast and slow regimes

We now examine whether the time-dependent Zeus-3D code confirms the co-existence of the solutions found in Section 3.1. For this purpose, we calculated time-dependent solutions for the three models shown in Table 1. We found almost the same Ω-interval of co-existence, with a difference of the order of ∼10⁻³ for Ω₁ and Ω₂.

Fig. 9 (top and middle panels) shows the results for the B2.5 V model for Ω = 0.75. Both time-dependent velocity profiles match their corresponding stationary solutions (black dashed lines) very closely. The bottom panel illustrates the density ratio as a function of r/R_* − 1, for the time-dependent (green solid line) and stationary (black dashed line) solutions.

Figure 9.

Top panel: Velocity profiles as a function of the inverse radial coordinate u. Blue and red solid lines represent the time-dependent fast and Ω-slow solutions, respectively, for Ω = 0.75. Black dashed lines are the corresponding stationary solutions. Middle panel: As the top panel, but now as a function of r/R_* − 1. Bottom panel: Density ratio between the time-dependent fast and Ω-slow solutions as a function of r/R_* − 1. In all panels, the stationary solutions are also shown, with black dashed lines, and are overlaid with the time-dependent ones.

5 SWITCH BETWEEN THE FAST AND SLOW VELOCITY REGIMES

With the purpose of studying the stability of the solutions in the co-existence region and a possible switch of wind regime, we introduced density perturbations at the base of the wind using Zeus-3D.

We conducted simulations using three types of density perturbations as a function of time: (i) a random perturbation where the perturbed base density, ρ_pert, takes random values around the steady-state base density, ρ₀, and the amplitude of the perturbation is expressed in terms of ρ_pert/ρ₀; (ii) a square perturbation where the value of ρ₀ switches from a fixed upper (or lower) value ρ_pert; and (iii) a Gaussian perturbation where ρ_pert increases (decreases) with respect to ρ₀, following a Gaussian profile.

Fig. 10 illustrates the behaviour of the three types of perturbations of the base density with the same enhanced amplitude and time. Time is given in units of the flow time, t = R_*/|$v$|_∞. Random perturbations represent sudden change; square perturbations obey a strong and sudden change; and Gaussian perturbations correspond to a smooth and continuous change in density at every instant.

Figure 10.

Typical shapes of perturbations with an enhanced base density. The enhanced base density is in units of the base density ρ₀ and the time is in units of the flow time t = R_*/|$v$|_∞. The perturbation time is 1000 ks.

Based on these perturbations, we searched for the minimum amplitude (ρ_pert/ρ₀) and the minimum time that could lead to a switch from a fast to a slow regime, and vice versa.

From the numerical simulations we performed, we concluded that we need to increase the base density, ρ₀, to switch from a fast solution to a Ω-slow one and to reduce ρ₀ for the other case. Tables 3 and 4 show the results of our simulations for each switching case. These tables also show the minimum amplitudes and time required to switch regime. All these perturbations last a period of time that is typically about 500 flow times for transitions from fast to Ω-slow solutions and about 100 flow times for transitions from Ω-slow to fast solutions. These period of times are similar to the time relaxation from the initial trial solution to a steady state.

Table 3.

Minimum conditions required from a density-enhanced perturbation to switch from a fast to an Ω-slow solution within the co-existence region. The models with no data correspond to simulations where there is no switch of solutions under any conditions.

Star	Time	Amplitude	Type
	[s]	[ρpert/ρ0]
B2.5 V	6 × 106	1 × 104.8	Random
	–	–	Square
	–	–	Gaussian
B0 IV	1 × 107	1 × 102.4	Random
	4 × 106	1 × 106	Square
	4 × 106	1 × 106	Gaussian
B3 I	5 × 107	1 × 102.7	Random
	–	–	Square
	–	–	Gaussian

Table 4.

Similar to Table 3, but here the switch is from an Ω-slow to a fast solution with a decreased density perturbation.

Star	Time	Amplitude	Type
	[s]	[ρpert/ρ0]
B2.5 V	5 × 105	1 × 10−2.9	Random
	1 × 104	1 × 10−2.0	Square
	5 × 104	1 × 10−2.0	Gaussian
B0 IV	–	–	Random
	2 × 105	1 × 10−2.0	Square
	9 × 105	1 × 10−2.0	Gaussian
B3 I	5 × 106	1 × 10−5.5	Random
	3 × 106	1 × 10−2.0	Square
	3 × 106	1 × 10−2.0	Gaussian

In the case of random perturbations, the model B2.5 V requires a much higher density amplitude Furthermore, for models B2.5 V and B3 I, square and Gaussian density perturbations do not switch at all, and the B0 IV model requires an unrealistic amplitude perturbation. When the regimes do not switch at all, a kink structure occurs in most of the studied cases.

When the Ω-slow solution is perturbed, the switch to the fast solution is attained by lowering the density by a factor ≲1/100 for square and Gaussian perturbations for all models. For random perturbations, the B0 IV model never switches to the fast solution, and the other two models require very high amplitudes; for example, the B2.5 V model requires an amplitude factor ≈1/1000 to switch. Based on the physically unlikely high values required to switch from one solution to another, both type of solutions seem very stable. However, the switch from a Ω-slow to a fast solution is more likely than one in the opposite direction.

The perturbation times required to switch regimes are quite different. To obtain an Ω-slow solution, starting from a fast solution, requires a perturbation time of the order of 4 × 10⁶–5 × 10⁷ s or 500–600 flow times. For the other switching case, the time is of the order of 10⁴–5 × 10⁶ s or a few–100 flow times. This difference is explained by the characteristic flow times of the two regimes to achieve a steady state from an initial trial solution.

Density and velocity contour-plots, starting from the fast regime, are shown in Fig. 11 for the B0 IV model, while those starting from the Ω-slow regime are shown in Fig. 12 for the B2.5 V model. In each case we show the results when a random perturbation was applied to the base density under the conditions given in Tables 3 and 4. In both cases, the solution switches its regime.

Figure 11.

Velocity contours as a function of time for the B0 IV model. Time is in units of the flow time, t = R_*/|$v$|_∞. The simulation begins from a fast solution with Ω = 0.75. Then a random enhanced density perturbation with amplitude 10^2.4 is applied and a switch to a Ω-slow solution is attained. Dashed lines indicate the initial and final times of the perturbation process.

Figure 12.

As Fig. 11, but this simulation corresponds to a B2.5 V model, starting with a Ω-slow solution rotating at Ω = 0.73. A random decreased density perturbation with amplitude 10^−2.9 is applied, and a switch to a fast solution is attained.

6 DISCUSSION

In the previous sections we studied various density perturbations that could yield a switch of wind regime in the co-existence region. In this section we discuss the implications of this regime-switching in relation to the variability of massive stars and their winds.

There is a consensus (Puls et al. 2008) that line-profile variability is induced by the interplay between disturbances in the photosphere, such as non-radial pulsations, stellar spots, etc. Thus, any mechanism (or a combination thereof) that induces a significant change in the photosphere density when the stellar rotational speed is in the co-existence range (Ω₁ ≤ Ω ≤ Ω₂) might have a dramatic impact in the equatorial plane, triggering a switch of wind regime. This could lead to the formation of a circumstellar disc in the equatorial plane of these stars, where the mass injection mechanism is given by radiation pressure in a Ω-slow wind regime. These discs might appear at rotational speeds of the order of Ω ∼ 0.65–0.75, values that are much lower than the critical rotational speed. Moreover, these rotational velocities are in better agreement with the observed mean value from the distribution of rotational speeds of Be stars (Zorec et al. 2016, and references therein). If subsequent perturbations do not have sufficient amplitudes, this scenario is stable in the frame of this 1D model. Otherwise, if the amplitudes of perturbations are high enough, there could be another switch to the fast regime, leading to dissipation of the disc. This new configuration could also be stable for a long period of time.

This regime-switching might also explain the change of phase between B and Be stars, and also episodic mass-loss events in Be stars (Balona 2010).

We can link this regime-switch with stellar evolution: as the internal structure of a star evolves, the critical rotational speed decreases (Ω increases, Langer 1998) and Ω might reach the co-existence region, creating or dissipating a disc. In a later evolutionary scenario characterized by Ω > Ω₂, only the Ω-slow regime wind is present at the equatorial plane (Curé 2004), leading to a different evolutionary phase with a higher mass-loss rate. Araya et al. (2017) demonstrated that whenever an Ω-slow regime is established, a stationary wind density structure is able to explain the observed Hα emission lines in Be stars.

In contrast to some radiation-driven viscous disc models (Kee, Owocki & Sundqvist 2016; Krtička, Owocki & Meynet 2011) that need near-critical rotational speed to create/dissipate a disc, our model requires under-critical rotation.

A current weakness of this m-CAK model is that it does not consider the role of viscosity and its influence on angular momentum transport, a mechanism that might explain a Keplerian disc. Furthermore, it is important to emphasize that the models presented here are purely 1D, and solutions in multiple dimensions might be quite different from what we present here. Thus, a more extensive study considering a 2D/3D model with non-radial forces (and the effects of stellar distortion and gravity darkening) is required to confirm the co-existence region.

7 CONCLUSIONS

The topology of the non-linear m-CAK differential equation predicts two wind solutions (from different solution branches) as a function of Ω for a given set of line-force parameters. The fast solution ceases to exist when Ω is at about a certain threshold value (0.65 ≲ Ω ≲ 0.76) that depends on the stellar and line-force parameters, while the Ω-slow solution begins to exist from this interval up to Ω ≲ 1.

We investigated the region in the Ω-space where these solutions are simultaneously present by solving 1D stationary and 1D time-dependent hydrodynamic equations at the equatorial plane. From the steady-state study we found a co-existence region inside a small interval of Ω, where both solutions are simultaneously present, both satisfying the same boundary conditions at the stellar surface (ρ(R_*) = ρ₀).

We demonstrated that time-dependent solutions are very sensitive to the initial solution topology. If the topology of the initial velocity profile is far from the asymptotic steady state (belonging to a different solution branch), the calculation leads to non-monotonic ‘kink’ solutions. Therefore, to obtain either fast or Ω-slow solutions with a globally monotonic acceleration, we need to use an initial solution representative of the fast or slow regime, respectively. Our calculations confirmed the stability of these known stationary solutions, and, when using the proper initial conditions, these solutions have no kinks.

The location and width of the co-existence region depend on stellar and line-force parameters. Higher values of α shift the location of this region towards higher values of Ω, while the width of the region is almost unchanged. For some particular values of δ, we find a gap in the steady-state solutions, confirming the results of Venero et al. (2016).

For the co-existence region, we determined which minimum conditions can induce a switch of regime by means of density perturbations at the base of the wind. We used random, square and Gaussian perturbations with different amplitudes and perturbation times. Our results indicate that the switching processes are not identical. While a switch from the fast to the Ω-slow regime requires a perturbation amplitude of some hundreds times the wind base density, the reverse case requires a hundredth of this base density.

Our main interpretation of this work is that the switching process can trigger the formation and/or dissipation of an outflowing equatorial disc without the need for a star to be rotating at almost the critical rotational speed.

ACKNOWLEDGEMENTS

The authors would like to thank the referee, Achim Feldmeier, for his thoughtful comments and suggestions. We also thank Diego Rial for his helpful comments concerning the topological analysis. IA acknowledges support from Fondo Institucional de Becas FIB-UV. MC and IA acknowledge support from Centro de Astrofísica de Valparaíso. AuD acknowledges support from NASA through Chandra Award numbers GO5-16005X, AR6-17002C, G06-17007B and TM7-18001X issued by the Chandra X-ray Observatory Center, which is operated by the Smithsonian Astrophysical Observatory for and on behalf of NASA under contract NAS8- 03060. AuD was supported by the FRS-FNRS at the University of Liège (Belgium) for a research stay. LC acknowledges financial support from CONICET (PIP 0177), La Agencia Nacional de Promoción Científica y Tecnológica (PICT 2016-1971) and the Programa de Incentivos (G11/137) of the Universidad Nacional de La Plata (UNLP), Argentina. LC and MC are grateful for support from the project CONICYT + PAI/Atracción de capital humano avanzado del extranjero (folio PAI80160057).

Footnotes

1

Madura et al. (2007) assume the formalism of Gayley (1995) with |$\bar{Q} \approx Q_{\mathrm{o}}$|⁠, contrary to what Puls et al. (2000) states for T_eff < 35 000 K. Independently of this, we use the same assumption with the purpose of comparing results.

APPENDIX A: TOPOLOGICAL ANALYSIS

In this section we present a topological analysis of a rotational radiation-driven wind. This analysis is based on the work of Curé & Rial (2007), where we presented the topological theory of radiation-driven winds in detail. Here we present a brief summary and refer the reader to Curé & Rial (2007) for more details.

A1 Topology of singular points

We define the logarithmic variables η and ζ in terms of u, w and w΄ as follows:

\begin{eqnarray} \eta = \ln (2 w^{\prime } / w) \end{eqnarray}

(A1)

and

\begin{eqnarray} \zeta = \ln (w w^{\prime }). \end{eqnarray}

(A2)

Then, from the equation of motion (equation 10) plus the regularity and singularity conditions at the singular point, we are able to solve for ζ = ζ(u, η, C΄) and then define two new functions, H(u, η, C΄) and R(u, η) (see section 3.4 in Curé & Rial 2007). A singular point occurs when the following conditions are simultaneously satisfied:

\begin{eqnarray} H(u,\eta ,C^{\prime })\,=\, R(u,\eta )\,=\,0. \end{eqnarray}

(A3)

Fig. A1 shows the H(u, η, C΄) and R(u, η) functions in terms of u, η and C΄. Dark-grey continuous lines show the contour plots H(u, η, C΄) = 0 for various values of the eigenvalue C΄ (labelled squares). In addition, the dashed contour lines are the H(u, η, C΄) = 0 function for C΄ = 70.37, which is the eigenvalue of the fast solution. The contour line corresponding to the eigenvalue C΄ = 70.37 has two components: one on the left of left component of the contour line of C΄ = 70 and the other inside the closed contour of C΄ = 70 at the right border of Fig. A1. The dotted closed line corresponds to the level H(u, η, C΄) = 0 for C΄ = 71.32, which is the eigenvalue of the Ω-slow solution. The function R(u, η) = 0 is also plotted in this figure, with solid thick black lines. The functions H(u, η, C΄) and R(u, η) intersect at five locations for both eigenvalues, three for the fast-solution eigenvalue and two for the Ω-slow-solution eigenvalue. All singular points are labelled from A to E clockwise, starting from the fast-solution eigenvalues as shown in Table A1. This table also show the location of each singular point in the u, η-plane and the value of the determinant of matrix B, which gives information concerning the topology type of a singular point (see the definition of matrix B in section 5.2 of Curé & Rial 2007). In all these cases, the determinant values were negative, and thus all singular points are X-type.

Figure A1.

H(u, η, C΄) and R(u, η) functions in terms of u, η and C΄; see text for details.

Table A1.

Locations and topology of the singular points.

Label	C΄	u	η	Det(B)	Topology
A (fast)	70.37	−0.97327	4.2516	<0	X-type
B	70.37	−0.02579	4.2221	<0	X-type
C	71.32	−0.04045	3.5097	<0	X-type
D (slow)	71.32	−0.04396	1.1126	<0	X-type
E	70.37	−0.03142	−0.1433	<0	X-type

A2 Integration from singular points

Consequently, we start to integrate from each singular point outwards and inwards to obtain the velocity profile as a function of u. Fig. A2 shows the velocity profiles of each of the integrations.

Figure A2.

Velocity profile as a function of u. We start from any singular point and integrate up- and downstream. Only fast and Ω-slow solutions (solid lines) reach the stellar surface. All other solutions have no physical meaning because they never reach the stellar surface.

Fast (A) and Ω-slow (D) solutions are plotted with black solid lines. The similarity with Fig. 9 is remarkable. The integration from singular points B, C and E are plotted with dashed, dotted and dot–dashed lines, respectively. None of these velocity profiles reaches the stellar surface, so we conclude that these are non-physical solutions of our rotating radiation-driven wind.

Finally, we note that although both physical solutions (fast and Ω-slow) found by Hydwind and Zeus-3D have the same initial condition at the stellar surface, the rotational speed, stellar and line-force parameter and the stationary solutions have different locations of the singular point and different eigenvalues.