A&A 428, 545-554 (2004)
DOI: 10.1051/0004-6361:20040325

The influence of rotation in radiation driven winds from hot stars

II. CAK topological analysis

M. Curé ¹ - D. F. Rial ²

1 - Departamento de Física y Meteorología, Facultad de Ciencias, Universidad de Valparaíso, Valparaíso, Chile
2 - Departamento de Matemáticas, Facultad de Ciencias Exactas y Naturales, Universidad de Buenos Aires, Argentina

Received 24 February 2004 / Accepted 26 August 2004

Abstract
The topological analysis from Bjorkman (1995) for the standard model that describes the winds from hot stars by Castor et al. (1975) has been extended to include the effect of stellar rotation and changes in the ionization of the wind. The differential equation for the momentum of the wind is non-linear and transcendental for the velocity gradient. Due to this non-linearity the number of solutions that this equation possess is not known. After a change of variables and the introduction of a new physically meaningless independent variable, we manage to replace the non-linear momentum differential equation by a system of differential equations where all the derivatives are explicitely given. We then use this system of equations to study the topology of the rotating-CAK model. For the particular case when the wind is frozen in ionization (${\delta =0}$) only one physical solution is found, the standard CAK solution, with a X-type singular point. For the more general case ( ${\delta \neq 0}$), besides the standard CAK singular point, we find a second singular point which is focal-type (or attractor). We find also, that the wind does not adopt the maximal mass-loss rate but almost the minimal.

Key words: hydrodynamics - methods: analytical - stars: early-type - stars: mass-loss - stars: rotation - stars: winds, outflows

1 Introduction

Since the launch of the first satellite with a telescope on board, it has been established the widespread presence of stellar winds from hot stars. These winds are driven by the transfer of momentum of the radiation field to the gas by scattering of radiation in spectral lines (Lucy & Solomon 1970). The theory of radiation driven stellar winds is the standard tool to describe the observed properties of the winds from these stars. Castor et al. (1975, hereafter "CAK'') obtained an analytical hydrodynamic model for these winds, based in the Sobolev approximation. The CAK model has been improved by Friend & Abbott (1986, "FA'') and Pauldrach et al. (1986, "PPK''), giving a general agreement with the observations. For a extended review see Kudritzki & Puls (2000, "KP'') and references therein.

This agreement with the observations led to the development of a new method to determine galactic distances using Supergiants as targets, namely the Wind Momentum Luminosity relationship ("WML'', Kudritzki et al. 1999; KP, and references therein).

More detailed studies from Puls et al. (1996) and Lamers & Leitherer (1993) came to the conclusion that the line-driven wind theory shows a systematic discrepancy with the observations. Lamers & Leitherer (1993) suggest that this discrepancy may arise due to an inadequate treatment of multiple scattering. Abbott & Lucy (1985), Puls (1987) and Gayley et al. (1995) have shown that multiple scattering can provide an enhancement of the wind momentum over that from single scattering only by a factor of two - three for O stars (Abbott & Lucy 1985, found a factor of 3.3 for the wind of $\zeta$ Pup).

Vink et al. (2000) calculate, including the multiple scattering effects, mass-loss rates for a grid of wind models that covers a wide range of stellar parameters. They found a much better agreement between theory and observation, concluding that the inclusion of multiple scattering increases the confidence of the WML relationship to derive extragalactic distances. In all the calculations involved in the WML relationship, the solution of the improved (or modified) CAK wind (hereafter m-CAK) is not used. Instead an ad hoc $\beta$-field velocity profile is utilized (see KP; Vink et al. 2000).

The unsatisfactory results of the velocity field obtained from the m-CAK model when applied in the WML relationship could come from the complex structure of this non-linear transcendental equation for the velocity gradient and its solution schema. Due to this non-linearity in the momentum differential equation, there exist many solution branches in the integration domain. A physical solution that describes the observed winds must start at the stellar photosphere, satisfying certain boundary condition and reach infinity. There is no solution branch that covers the whole integration domain, thus a solution must pass through a singular point in order to match a second solution branch. Therefore, the solution in this second solution branch reaches infinity. To find the location of singular points is one of the most difficult aspects of topological analysis of non-linear differential equations.

Bjorkman (1995) performed a topological analysis of the CAK differential equation. He showed that the original solution from CAK, which passes through a X-type critical point and has a monotonically increasing velocity field, is the only physical solution that satisfies the condition of zero pressure at infinite radius. In this study Bjorkman did not include the influence of the star's rotational speed.

Although it is known that the line-force parameters (see below) are not constant through in the wind (Abbott 1982), the standard m-CAK model still uses these parameters as constant. For the particular case of extreme low metalicities, Kudritzki (2002) introduced a new treatment of the line-force with depth dependent radiative force multipliers. As a test, he applied this new treatment for the most massive and most luminous O stars in the Galaxy and in the Magellanic Clouds (due to the lower metalicity) finding an acceptable agreement between theory and observations. Then, it was used to understand the influence of stellar winds on the evolution of very massive stars in the early universe and on the interstellar medium in the early stages of galaxy formation.

On the other hand, it is known from observations, that all early type stars have moderate to large rotational speeds (Hutchings et al. 1979; Abt et al. 2002) and for Oe and Be stars, their rotational speed is a large fraction of their break-up speed (Slettebak 1976; Chauville et al. 2001). The incorporation of rotation in the CAK and m-CAK models has been studied by Castor (1979), Marlborough & Zamir (1984), FA and PPK, concluding that the effect of the centrifugal force results in a downstream-shift of the position of the singular point, a slightly lower terminal velocity and a slightly larger mass loss rate as a result of an increasing in the star's rotational speed. Maeder (2001) studied the influence of the stellar rotation in the WML relationship, finding just a very small effect on it.

A revision of the influence of the stellar rotation in radiation driven winds has been done by Curé (2004) finding that there exists a second singular point in these winds. He studied the case when the stellar rotational speed is high and found numerical solutions, that pass through this second singular point, and which are denser and slower than the standard m-CAK solution.

In view of these results, it is crucial to understand the solution topology of the standard model, forall when one wants to incorporate other physical processes into the theory.

The purpose of this article is to study the topology of the rotating-CAK model. In Sect. 2 we give a brief exposition of the radiation driven winds theory and the non-linear differential equation for the momentum, including rotation, is shown. In Sect. 3 after a coordinate transformation, we develop a general method that allows to replace the non-linear momentum equation in a simple and straightforward manner by a system of ordinary differential equations, where all the derivatives are explicitely given. In Sect. 4 a general condition for the eigenvalue of the problem is developed. This condition allows to classify the topology of the singular point (Saddle or Focal) and constrains the location of it in the integration domain. Section 5 is devoted to the application of the criteria developed in Sect. 4 for the rotating-CAK model. In Sect. 6 we show numerical results of this topological analysis, first for a wind frozen in ionization (setting the $\delta$ parameter of the line-force to zero) and compare our results with Bjorkman`s (1995) for the non-rotating CAK model and the rotating-CAK model from Marlborough & Zamir (1984). Furthermore, Sect. 6 shows the results of the influence of changes in the wind`s ionization structure ( ${\delta \neq 0}$) on the topology and discuss the rotating-CAK wind model. Conclusion are in Sect. 7.

   
2 The non-linear differential equation

The standard stationary model for radiation driven stellar winds treats an one-component isothermal radial flow, ignoring the influence of heat conduction, viscosity and magnetic fields (see e.g., Kudritzki et al. 1989, "KPPA'').

For a star with mass M, radius $R_{\ast}$, effective temperature T and luminosity L, the momentum equation with the inclusion of the centrifugal force due to star's rotation, reads:

$$v\frac{{\rm d}v}{{\rm d}r}=-\frac{1}{\rho }\frac{{\rm d\wp}}{... ... \frac{v_{\phi }^{2}(r)}{r}+g^{\rm line}(\rho, v', n_{\rm E})$$ (1)

where v is the fluid velocity, $v' = {\rm d}v/{\rm d}r$ is the velocity gradient, $\rho$ is the mass density, $\dot{M}$ is the star's mass loss rate, $\wp$ is the fluid pressure, $v_{\phi } = v_{\rm rot} ~ R_{\ast} / r$, where $v_{\rm rot}$ is the star's rotational speed at the equator, $\Gamma$ represents the ratio of the radiative acceleration due to the continuous flux mean opacity, $\sigma_{\rm e}$, relative to the gravitational acceleration, i.e., $\Gamma= \sigma_{\rm e} L / 4 \pi c G M$ and the last term  $g^{\rm line}$ represents the acceleration due to the lines.

The standard parameterization of the line-force (Abbott 1982) reads:

$$g^{\rm line}=\frac{C}{r^{2}}\;f_{\rm D}(r,v,v')\;\left( r^{2}v\; v'\right) ^{\alpha }\;\left( {n_{E}}/{W(r)}\right) ^{\delta }$$ (2)

W(r) is the dilution factor and $f_{\rm D}(r,v,v')$ is the finite disk correction factor. The line force parameters are: $\alpha$, $\delta$ and k (the last has been incorporated in the constant C), typical values of these parameter from LTE and non-LTE calculations are summarized in Lamers & Cassinelli (1999, Chap. 8, "LC''). We have not used the absolute value of the velocity gradient in the line-force term, because we are interested in monotonic velocity laws.

The constant C represents the eigenvalue of the problem (see below) and is given by:

$$C=\Gamma G M k \left( \frac{4\pi }{\sigma _{\rm E} v_{\rm th} \dot{M}}\right) ^{\alpha }\; ,$$ (3)

where $v_{\rm th}=({2k_{\rm bol}T/m_{\rm H}})^{1/2}$ is the hydrogen thermal speed, $n_{\rm E}$ is the electron number density in units of $10^{-11}~\rm cm^{-3}$ (Abbott 1982), while the meaning of all other quantities is standard (see, e.g., LC).

Together with the momentum Eq. (1), the continuity equation reads:

$$\; 4 \pi r^{2}\rho \;v = \dot{M}.$$ (4)

Introducing the following change of variables:

$$u=-R_{\ast}/r ,\\$$ (5)

$$w= \; v/a ,\\$$ (6)

$$w'=\; {\rm d}w/{\rm d}u ,$$ (7)

where a is the isothermal sound speed, i.e., $\wp=a^{2}\rho$. Replacing the density $\rho$ from Eq. (4), the m-CAK momentum Eq. (1) with the line force (2), we obtain:

 

$$F(u,w,w^\prime)\equiv \left( 1-\frac{1}{w^{2}} \right) w \; w'+ \... ...verline{C} \;f_{\rm D}\;g(u)(w)^{-\delta }\left( w \;w'\right)^{\alpha}\ = \; 0$$ (8)

here the constants A, $\overline{C}$ and $\Omega$ are:

$$A = \frac{GM(1-\Gamma )}{a^{2}R_{\ast }}=\frac{v_{\rm esc}^{2}}{2a^{2}} =\frac{v_{\rm break\mbox{-}up}^{2}}{a^{2}} ,$$ (9)

$$\overline{C} = C\;\left( \frac{\dot{M} D}{2\pi}\frac{10^{-11}... ...\right)^{\delta }\;\left(a^{2}R_{\ast }\right)^{(\alpha -1)} ,$$ (10)

$$\Omega=\frac{v_{\rm rot}}{v_{\rm break\mbox{-}up}}= \frac{a_{\rm rot}}{\sqrt{A}} ,$$ (11)

where $a_{\rm rot}$ is defined by:

$$a_{\rm rot}=\frac{v_{\rm rot}}{a} ,$$ (12)

here $v_{\rm esc}$ is the escape velocity and $v_{\rm break\mbox{-}up}$ is the "break-up'' velocity. The function g(u) is defined as:

$$g(u)=\left( {1+\sqrt{1-u^{2}}}\right)^{\delta } ,$$ (13)

and the constant D is:

$$D= \frac{1}{m_{\rm H}} \frac{(1+Z_{\rm He} Y_{\rm He})}{(1+A_{\rm He} Y_{\rm He})} ,$$ (14)

$Y_{\rm He}$ is the helium abundance relative to the hydrogen, $Z_{\rm He}$ is the amount of free electrons provided by helium, $A_{\rm He}$ is the atomic mass number of helium and $m_{\rm H}$ is the mass of the proton.

The standard solution, from this non-linear differential Eq. (8), starts at the stellar surface and after crossing the singular point reaches infinity. At the stellar surface the differential equation must satisfy a boundary condition, namely the monochromatic optical depth integral (see Kudritzki 2002, Eq. (48)):

$$\tau_{\rm Phot}= \int_{R_{\ast}}^{\infty} \sigma_{\rm E} ~\rho ~(r) ~{\rm d}r ~ =~ \frac{2}{3}\cdot$$ (15)

A numerically equivalent boundary condition is to set the density at the stellar surface to a specific value,

$$\rho~(R_{\ast}) = \rho_{\ast}.$$ (16)

When the singularity condition,

$$\frac{\partial}{\partial w'}F(u,w,w')= F_{w'} = 0$$ (17)

is satisfied, its location ( $u=u_{\rm c}$) corresponds to a singular (or critical) point, and in order to get a physical solution, the regularity condition must be imposed, namely:

$$\frac{\rm d}{{\rm d}u} F(u,w,w')=\frac{\partial F}{\partial u}+\frac{\partial F}{\partial w}w'= F_{u}+ F_{w}~w'=0 ,$$ (18)

hereafter, all partial derivatives are written in a shorthand form, i.e. $F_{u}={\partial F}/{\partial u}$.

In order to satisfy simultaneously Eqs. (8), (17), (18) and (15) or (16), the value of the constant $\overline{C}$ is not arbitrary, i.e., the constant $\overline{C}$ is the eigenvalue of this non-linear problem.

  
3 Coordinate transformation and system of equations

In this section we will apply another coordinate transformation and introduce a new independent variable, $\tau$, without physical meaning. This will allow us to transform the non-linear differential equation for the momentum (8) into a system of coupled differential equations, which is numerically integrable.

Defining

$$y = \frac{1}{2} \;w^{2} ,$$ (19)

and

$$p =w\;w'=\frac{{\rm d}y}{{\rm d}u} ,$$ (20)

the momentum Eq. (8) can be written in a general form as:

F(u, y, p) = 0. (21)

Differentiating this function and using ${\rm d}y=p~{\rm d}u$, we obtain:

$${\rm d}F=\left(F_{u}~+ p~F_{y}\right) ~{\rm d}u + F_{p}~{\rm d}p = 0.$$ (22)

We introduce now a new independent variable, $\tau$, defined implicitely by:

$${\rm d}u=F_{p}~{\rm d}\tau .$$ (23)

Because u and $\tau$ are independent variables, they have to be related between them. While $F_{p}\neq 0$, we can write $\tau$ as function of u as:

$$\frac{{\rm d}y}{{\rm d}u}=\frac{{\rm d}y}{\rm d\tau} \left(\frac{{\rm d}u}{\rm d\tau}\right)^{-1}~ = \;p.$$ (24)

We can transform from Eq. (22) to the following system of ordinary differential equations:

$$\frac{{\rm d}u}{\rm d\tau} \equiv U= F_{p} ,$$ (25)

$$\frac{{\rm d}y}{\rm d\tau} \equiv Y=p~F_{p} ,$$ (26)

$$\frac{{\rm d}p}{\rm d\tau} \equiv P=-\left( F_{u}~+p~F_{y}\right).$$ (27)

A solution of this system of differential equations is also a solution of the original momentum equation, since if any initial condition (u₀,y₀,p₀) satisfies F(u₀,y₀,p₀) =0, then a solution of Eqs. (25)-(27) verifies that $F(u(\tau),y(\tau),p(\tau))=0$.

An advantage of this equation system (25)-(27) over the CAK momentum differential Eq. (8) is that all the derivatives are explicitely given, therefore there is no need to use root-finding algorithms to find the value of the velocity-gradient. Also, standard numerical methods (e.g., Runge-Kutta) can be used to integrate this system.

   
4 Linearization and eigenvalue criteria

All critical points of the system (25)-(27) satisfy simultaneously F=0, U=0 and P=0. Thus, in order to study the behavior of the solution in the neighborhood of a singular point we linearise this system of differential equations, using the Groebman-Hartman theorem ("GH'', see appendix; For more details see Amann 1990), we obtain:

$${\rm d}U =\left( U_{u}~+pU_{y}\right) ~{\rm d}u+U_{p}~{\rm d}p,$$ (28)

$${\rm d}P =\left( P_{u}~+pP_{y}\right) ~{\rm d}u+P_{p}~{\rm d}p.$$ (29)

We have not included ${\rm d}Y$ because U and Y are dependent (Y=pU, Eq. (26)). The GH theorem indicates that the eigenvalues (and the eigenvectors) of the partial derivative matrix B, defined by:

$$B=\left( \begin{array}[c]{cc} U_{u}+p~U_{y} & U_{p}\\ P_{u}+p~P_{y} & P_{p}% \end{array}\right) ,$$ (30)

provide the information concerning the topology structure of the critical (singular) points.

On the other hand, if we consider the eigenvalue $\bar{C}$ as a free paramenter, the critical points are the solution of the following system of equations:

$$F\left(u,y,p,\bar{C}\right) = 0,$$ (31)

$$U\left(u,y,p,\bar{C}\right) = 0,$$ (32)

$$P\left(u,y,p,\bar{C}\right) = 0.$$ (33)

In this case the number of incognits is greater than the number of equations, then we can only solve for the incognits in terms of one of them (implicit function theorem).

If $\left(u_{\rm c},y_{\rm c},p_{\rm c},\bar{C}\right)$ satisfies the previous system and furthermore, $\Delta \neq 0$ at the singular point, we can solve for $y_{\rm c}$, $p_{\rm c}$ and $\bar{C}$ and its derivatives as a function of $u_{\rm c}$. We obtain for the gradient of $\bar{C}$:

$$\frac{{\rm d}\bar{C}}{{\rm d}u_{\rm c}}=-\frac{1}{\Delta}\det... ... U_{p} & U_{u}\\ P_{y} & P_{p} & P_{u}% \end{array}\right), %$$ (34)

where the determinant $\Delta$ is defined by:

$$\Delta=\det\left( \begin{array}[c]{ccc}% F_{y} & F_{p} & F_{\... ...U_{\bar{C}}\\ P_{y} & P_{p} & P_{\bar{C}} \end{array}\right),$$ (35)

Considering that at the critical point: F_p=0 and F_u=-pF_y, Eq. (34) transforms to:

$$\frac{{\rm d}\bar{C}}{{\rm d}u_{\rm c}}=-\frac{1}{\Delta}\det... ... \end{array}\right) =\frac{F_{y}}{\Delta}\det\left( B\right) .$$ (36)

The GH theorem establishes that the singular point topology is determined by the sign of the eigenvalues of the matrix B. A critical point is of X-type (Saddle) when the eigenvalues are both real and have opposite sign, or equivalently,

$$\det(B) ~ < ~0 .$$ (37)

Applying this to the rotating-CAK model (see next section for details) it is verified that,

$$\frac{F_{y}}{\Delta} ~ > ~0 ,$$ (38)

then a X-type singular point corresponds to the condition:

$$\frac{{\rm d}\bar{C}}{{\rm d}u_{\rm c}} ~< ~0 .$$ (39)

Therefore any X-type physically relevant singular point, is related to the behavior of the eigenvalue $\bar{C}$, specifically when condition (39) holds.

   
5 The topology of the rotating-CAK model

The analysis we have shown in the last section (up to Eq. (37)) is valid for the general radiation driven stellar winds theory, i.e., including the finite-disk correction factor. In the remainder of this paper we will study the original CAK model, i.e. $f_{\rm D}=1$. The topological study of the m-CAK model will be the scope of a future article.

The non-linear momentum Eq. (8) for the rotational-CAK model (including $\delta$), in (u, y, p) coordinates reads:

 

$$F(u,y,p) =\left( 1-\frac{1}{2y}\right) p+\frac{2}{u}+A\left(1+\Omega^{2}u\right) -\bar{C}~g\left( u\right) ~y^{-\delta/2}p^{\alpha} = 0 ~.$$ (40)

From this equation, after a straightforward calculation, the corresponding system of differential equations U and P, (Eqs. (25) and (27), respectively) are:

$$U\left( u,y,p\right) = 1-\frac{1}{2~y}-\alpha~\bar{C}~g(u)~y^{-\delta/2}~p^{-1+\alpha} ,$$ (41)

 

$$P\left(u,y,p\right) =\frac{2}{u^{2}}-\Omega^{2}A-\frac{p^{2}}{2~y... ...lta/2}~p^{1+\alpha} +\bar{C}~g\left( u\right) ~h(u)~~y^{-\delta/2}~~p^{\alpha},$$ (42)

where

$$h\left( u\right) =-~\delta~u\left(1-u^{2}+\sqrt{1-u^{2}}\right)^{-1} ,$$ (43)

F_y and $\Delta$ are:

$$F_{y}=\frac{p}{2~y^{2}}+\frac{\delta}{2}\bar{C}~g\left( u\right) ~~y^{-(1+\delta/2)}~~p^{\alpha} ,%$$ (44)

$$\Delta (1-\alpha)~g(u)~y^{-\delta/2}~p^{\alpha} \;\;$$ =

$$\times \bigg[ \frac{p}{2y^{4}} +\frac{1}{4}~\bar{C}~g(u)~y^{-(3+\delta)}~p^{\alpha-1} \;\;$$

$$\times \left( (4\alpha+\delta)~y^{\delta/2}~ + 2~\alpha~\delta~\bar{C}~g(u)~y\;p^{\alpha } \right) \bigg].$$ (45)

It is easy to verify that both, $F_{y}\; > \; 0$ and $\Delta\; > \; 0$.

Once the location of the critical point, $u_{\rm c}$, is known, we can solve $y_{\rm c}$, $p_{\rm c}$ and $\bar{C}$ from Eqs. (40)-(42), obtaining:

$$y_{\rm c}=\frac{1}{2}+\frac{\alpha}{\left( 1-\alpha\right) }\... ...rac{2}{u_{\rm c}}+A\left(1+\Omega^{2}u_{\rm c}\right)\right) ,$$ (46)

$$p_{\rm c}=\frac{1}{2}q_{\rm c}+\frac{\alpha}{1-\alpha}\left( \frac{2}{u_{\rm c} }+A\left(1+\Omega^{2}u_{\rm c}\right)\right) ,$$ (47)

$$\bar{C}=\frac{~y_{\rm c}^{\delta/2}}{\left( 1-\alpha\right) g... ...rac{2}{u_{\rm c}}+A\left(1+\Omega^{2}u_{\rm c}\right)\right) ,$$ (48)

where $q_{\rm c}$ is the positive solution of the quadratic equation

$$q^{2}+\frac{\delta}{1-\alpha}\left( \frac{2}{u_{\rm c}}+A\lef... ...a^{2}u_{\rm c}\right)\right) q=\gamma\left( u_{\rm c}\right) ,$$ (49)

with

 

$$\gamma\left( u_{\rm c}\right)=-2 A \Omega^{2} + \frac{4}{u_{\rm c... ...1-\alpha}\left( \frac{2}{u_{\rm c}}+A\left(1+\Omega^{2}u_{\rm c}\right)\right).$$ (50)

Equations (46)-(48) are generalizations of the Eqs. (49)-(51) from Marlborough & Zamir (1984), including now the effect of the line-force term $(n_{\rm E}/W)^{\delta}$.

Since the eigenvalue $\bar{C}$ must be positive, we have:

$$\frac{2}{u_{\rm c}} + A\left(1+\Omega^{2}u_{\rm c}\right)~>~0 .$$ (51)

This last inequality imposes a restriction in the position (furthest from the stellar surface in the radial coordinate r) of the singular point, namely:

$$u_{\rm c}~ < ~ u_{\max} ~\equiv~ -\frac{4}{A\left( 1+\sqrt{1-8~\Omega^{2}/A}\right)} \cdot$$ (52)

   
6 Topology of the rotating CAK wind

In this section we show the results of our topological analysis. Following Bjorkman (1995), we choose a typical B2  V star with stellar parameters summarized in Table 1. We have adopted the line-force multiplier parameters k and $\alpha$ from Abbott (1982) and the $\delta$ parameter is from LC, these values are summarized in Table 2.

Table 1: B2 V stellar parameters.

Table 2: Line-force parameters.

6.1 The frozen-in ionization ( ${\delta =0}$)

The factor $(n_{\rm E}/W)^{\delta}$ in the line-force takes into account the changes in the ionization of the wind. As a first step to understand the topology of the rotating-CAK model we set ${\delta =0}$.

6.1.1 The critical point interval

This case is simpler because we can obtain analytically from Eqs. (46)-(48) the variable y as function of u and p, i.e.:

$$y=\frac{p~}{2~\left( 2/u+A(1+\Omega^{2}~u)~+p~-C~p^{\alpha}~\right) }\cdot$$ (53)

Bjorkman (1995) obtained the same result (see his Eq. (13)) for the non-rotational case, $\Omega =0$. Moreover, $q_{\rm c}$ can be expressed by:

$$q_{\rm c}^2=\gamma(u_c)= 4 u_{\rm c}^{-2} - 2 A ~ \Omega^{2} .$$ (54)

As we mentioned in previous sections, q_c must be positive. This gives an additional restriction for the location of the singular point (nearest to the stellar surface in the radial coordinate r). The value of $u_{\rm min}$ is defined when $\gamma(u_{\rm min})=0$ and is given by:

$$u_{\min} =\max\left\{ -1,-\frac{\sqrt{2/A}}{\Omega}\right\} \cdot$$ (55)

Therefore, $u_{\rm c}$ is restricted to the interval $(u_{\rm min}, u_{\rm max})$. A similar result which restricts the location of the critical point to an interval has been found by Marlborough & Zamir (1984) see their Eq. (65)).

Figure 1a shows the behavior of $\gamma$ (Eq. (54)) versus $r/R_{\ast }-1$ for different values of the rotational parameter $\Omega$. For the non-rotational case $\Omega =0$, the function $\gamma$ is positive for almost the whole integration domain, therefore the singular point can be placed anywhere. Thus, is the lower boundary condition which fixes the position of the singular point. For the rotational case, the second term in the RHS of Eq. (54) is the dominant term for almost any value of $\Omega$. Therefore the larger is $\Omega$ the larger is the value of  $r_{\rm min}$ ( $\equiv -1/u_{\rm min}$) as the different curves in Fig. 1a show.

On the other hand, the value of the term $8 \Omega^2/A$ in Eq. (52) is almost negligible and consequently the value of $r_{\rm max}(\equiv\! -1/u_{\rm max})$ is almost constant at $r=A R_{\ast}/ 2$. Table 3 shows the interval $(r_{\rm min},r_{\rm max})$ for different values of the parameter $\Omega$. It is clear from this table that from a very low value of $\Omega=\sqrt{2/A}$ (=0.038 for our test star) the location of $r_{\rm min}~>~ R_{\ast}$ and the critical point, $r_{\rm crit}$, is strongly shifted downstream in the wind. Once  $r_{\rm crit}$ (or $u_{\rm c}$) is fixed, the value of $\gamma(r_{\rm crit})$ is inserted in Eq. (54) and $q_{\rm c}$ and $\bar{C}$ are obtained. Figure 1b shows $\bar{C}$ from Eq. (48) against $r/R_{\ast }-1$ for the same values of $\Omega$ of Fig. 1a (same type of lines too). The value of the derivative, ${\rm d}\bar{C}/{\rm d}u$ is always negative indicating that the critical point is an X-type.

Figure 1: a) The function $\gamma$ versus $r/R_{\ast }-1$ (Eq. (54)) for different values of the rotational parameter $\Omega$. $\Omega =0$, continuous-line; $\Omega =0.25$, dashed-line; $\Omega =0.8$, dotted-line. The horizontal line is for $\gamma =0$. b) The eigenvalue from Eq. (48) as a function of $r/R_{\ast }-1$ for the same values of $\Omega$ as a). Note that ${\rm d}\bar{C}/{\rm d}u$ is always negative in the interval  $(r_{\rm min},r_{\rm max})$.

Open with DEXTER

Table 3: Analytical approximation and numerical results for the rotational-CAK model with ${\delta =0}$. Note: the mass loss rate is given in units of $10^{-9}~ M_{\odot} / \rm year$ and the terminal velocity is in km s^-1.

6.1.2 Linearization of the B-matrix: Eigenvalues and eigenvectors

The matrix of linearization B (Eq. (30)) is given by:

$$\left. B~\right\vert _{\left( u_{\rm c},y_{\rm c},p_{\rm c}\r... ...\eta\right) }{u_{\rm c}^{2}} & -2\theta^{2} \end{array}\right)$$ (56)

where $\theta$ and $\eta$ are given by:

$$\theta=\sqrt{1-\Omega^{2}~A~u_{\rm c}^{2}/2} ,$$ (57)

and

$$\eta=-\left( 2+A~u_{\rm c}~\left(1+\Omega^{2}~u_{\rm c}\right)\right) /\left( 1-\alpha\right) .$$ (58)

The eigenvalues of the B-matrix are:

$$\mu_{\pm}=\frac{1}{u_{\rm c}^{2}p_{\rm c}}\left( -\theta^{2}\... ...ft( 1-\alpha\right) \left( 1+2\theta^{3}\right) \eta }\right).$$ (59)

It is easy to check that the product of the eigenvalues $\mu_{+}\times \mu_{-}$ is negative, i.e., the topology of the singular point is saddle or X-type as the GH theorem establishes. We can write the associated eigenvectors as $\left( 1,\nu_{\pm}\right)$, where

$$\nu_{\pm}=\frac{\left( -3\theta^{2}\pm\sqrt{9\theta^{4}+4\alp... ...a\right) }{\alpha\left( 1-\alpha\right) u_{\rm c}% ^{2}\eta} ,$$ (60)

where $\nu_{+}$ ($\nu_{-}$) corresponds to the unstable (stable) manifold. The maximum value of $\bar{C}$, that accounts for the minimum mass loss rate, occurs when $u_{\rm c}=u_{\min}$, Eq. (48) becomes:

$$\begin{array}[l]{lll} \bar{C}_{\max} =&\frac{A~(1-\Omega^{2})... ...ight) ^{1-\alpha} &, \;if\; \Omega ~> ~ \sqrt{2/A}. \end{array}$$ (61)

6.1.3 Phase diagram

In order to obtain a numerical solution for the wind, we need an approximation for the location of the critical point. We use the value of $u_{\rm c}=u_{\rm min}$ as a first guess of the critical point and use this value of $u_{\rm c}$ to calculate $\bar{C}_{\rm max}$ as our first guest for the eigenvalue $\bar{C}$. Table 3 shows the values of  $r_{\rm min}$, $r_{\rm crit}$ and $\bar{C}_{\rm max}$, $\bar{C}$ confirming that this is a very good approximation. Now using standard numerical algorithms we integrate our system of equations (Eqs. (25)-(27)) from the singular point up and downstream to obtain the numerical solution.

As Bjorkman (1995) pointed out, it is insightfull to study the solution topology in a p versus $(r/R_{\ast }-1)$ diagram. Figure 2 shows this phase diagram.

	Figure 2: The topology of the freeze in ionization rotating-CAK model (${\delta =0}$), p versus $(r/R_{\ast }-1)$ for $\Omega =0.25$. The unique curve that starts at the stellar surface and reaches infinity is the CAK original solution (continuous-line).
Open with DEXTER

If we start to integrate at the singular point, we cannot leave this point because all the equations, U=0, Y=0 and P=0 are simultaneously satisfied at this point. Therefore we have to move slightly up and downstream along the direction of the unstable manifold. After this, we can integrate obtaining the different solutions showed in Fig. 2. From this figure, it is clear that the only solution that reaches the stellar surface ( $\tau \rightarrow +\infty$, in the direction of $-(1,\nu_{+})$) and also reaches infinite ( $\tau \rightarrow +\infty$, in the direction of $(1,\nu_{+})$) is the original CAK solution (continuous-line). The results from Table 3 for the non-rotational case are the same one obtained by Bjorkman (1995).

	Figure 3: The velocity profile v(r) (in km s-1) as function of $r/R_{\ast }-1$ for the non rotating case (continuous-line) and for the rotational cases: $\Omega =0.25$ dashed-line; $\Omega =0.8$ dotted-line. The location of the singular points is indicated by a larger-dot.
Open with DEXTER

Figure 3 shows the velocity profile, v (in km s^-1) versus $r/R_{\ast }-1$ for our B2  V test star. We have chosen this value of $\Omega =0.25$ from the study of Abt et al. (2002), that concluded that B-stars rotate at a $25\%$ of their break-up speed. The values $\Omega =0.8$ accounts for a fast rotator, e.g., a typical Be-Star (Chauville et al. 2001). We see from this figure that neglects the rotational speed always overpredicts the value of the terminal velocity.

We conclude that the rotational speed shifts the location of the critical point downstream and reduces the terminal velocity, but has almost no influence on the value of the eigenvalue (mass loss rate). Furthermore, we can see from our approximate and numerical results summarized in Table 3 that the CAK wind do not have the maximum mass-loss rate as Feldmeier et al. (2002) and Owocki & ud-Doula (2004) concluded for a non-rotating CAK model with zero sound speed. Contrary to expectation, the rotating-CAK wind critical solution corresponds to an almost minimum mass-loss rate (maximum eigenvalue).

6.2 The rotating CAK model ( ${\delta \neq 0}$)

Abbott (1982) studied the indirect influence of the density in the line-force through the dependence of the ionization balance in the electron density. He found that this dependency modifies the force-multiplier and therefore the line-force by a factor $(n_{\rm E}/W)^{\delta}$, where $n_{\rm E}$ is the electron density (in units of $10^{-11} \rm ~gr~/ cm^{3}$) and W is the dilution factor.

Although this is a weak influence, because $\delta$ ranges between 0.0 and 0.2, it is important to study how its inclusion in the momentum equation modifies the topology of the rotating CAK model.

As we pointed out in Sect. 5, the existence of $q_{\rm c}~>~0$ implies $\gamma\left(u_{\rm c}\right)~>~0$.

	Figure 4: Same as Figs. 1a and b but for $\delta =0.02$. See text for discussion.
Open with DEXTER

The behavior of $\gamma$ (Eq. (50)) versus $r/R_{\ast }-1$ is shown in Fig. 4a for the same values of $\Omega$ of Fig. 1. We can interpret from this figure, that $\gamma$ in $(R_{\ast},r_{\rm max})$ is a decreasing function in the neighborhood of $r \rightarrow R_{\ast}^{+}$ and an increasing function in the neighborhood of $r \rightarrow r_{\rm max}$. Furthermore, the function $\gamma$ posses only one minimum in the integration domain located at $r \equiv r_{\rm min}$. The value of $\gamma(r_{\rm min})$ decreases as $\Omega$ increases. For small values of $\Omega$, $\gamma(r_{\rm min})~>~0$ and the location of the singular point  $r_{\rm crit}$ can be anywhere in the integration domain $(R_{\ast},r_{\rm max})$. For larger values of $\Omega$, $\gamma(r_{\rm min})~<~0$ and $\gamma(r)=0$ has two roots, at $r \equiv r_{\pm}$ where r_- <  r₊ (or u_- <  u₊). Thus, $r_{\rm crit}$ can be located now in two different intervals, i.e., $r_{\rm crit} \in (R_{\ast},r_{-})$ or $r_{\rm crit} \in (r_{+},r_{\rm max})$. For the particular case where $\gamma(r_{\rm min})~=~0$, the value of the rotational speed parameter is $\Omega~\equiv~\Omega_{\rm bif}$, where the subscript "bif'' accounts for the bifurcation in the wind solution topology.

Figure 4b shows $\bar{C}$ (Eq. (48)) against $r/R_{\ast }-1$ for different values of $\Omega$. When $\gamma(r)$ is positive, $\bar{C}$ exhibits now a different behavior compared with the ${\delta =0}$ case. As long as $\Omega < \Omega_{\rm bif}$, ${\rm d}\bar{C}/{\rm d}u$ is positive in the interval $(R_{\ast},r_{\rm min})$, i.e., any singular point in this interval is an attractor. Furthermore, ${\rm d}\bar{C}/{\rm d}u$ reaches its maximum, when $\gamma$ is minimum, i.e., in the neighborhood of $r=r_{\rm min}$, from this point up to $r_{\rm max}$, ${\rm d}\bar{C}/{\rm d}u ~< ~0$ and any singular point in $(r_{\rm min},r_{\rm max})$ is X-type. When $\Omega > \Omega_{\rm bif}$, the minimum of $\gamma(r)$ is negative and the intervals are reduced to: $(R_{\ast},r_{-})$ for the attractor type singular point and $(r_{+}, r_{\rm max})$ for the X-type singular point.

The behavior of $y_{\rm c}$, $p_{\rm c}$ and $\bar{C}$ (Eqs. (46)-(48)) in the neighborhood of $r~\rightarrow r_{\pm}$, in the inverse radial coordinate u, is as follows:

 

$$y_{\rm c} \sim \left\vert u-u_{\pm}\right\vert ^{-1}~\rightarrow ~ \infty ,$$ (62)

$$p_{\rm c} \rightarrow \frac{\alpha}{1-\alpha}\left( \frac{2}{u_{\pm}} +A (1+\Omega^{2}u_{\pm})\right) ~ ,$$ (63)

$$\bar{C} \sim \left\vert u-u_{\pm}\right\vert ^{-\delta/2}~\rightarrow ~+\infty.$$ (64)

This divergent behavior of $\bar{C}$ in the neighborhood of $r_{\pm }$, is shown in the curves for $\Omega \ne 0$ in Fig. 4. The almost constant value of $\bar{C}$ explains why a slightly change in the eigenvalue cause an enormous change in the location in the singular point.

Tables 4 and 5 summarise the numerical calculation for our test star with $\delta =0.02$ and $\delta =0.1$ respectively. The data of the $\dot{M}$ column (see also the $\bar{C}$ column) show that the effect of the rotation on $\dot{M}$ ($\bar{C}$) is almost negligible. Figure 5 shows the velocity profile for three different values of $\Omega$ ( 0.0;0.25;0.8), panel a) for $\delta =0.02$ and panel b) for $\delta =0.1$. A large dot shows the respective positions of the critical points. The effect of shifting the position of the critical point is stronger for low rotational speeds and decreases when $\Omega$ increases as a comparison between Figs. 3 and 5 clearly shows. The terminal velocity, is a decreasing function of the rotational speed and has almost the same behavior as in the ${\delta =0}$  case. But for high rotational speeds, the influence of $\delta$ in  $v_{\infty}$ is negligible/small as a comparison between Tables 3 and 4/5 shows.

We conclude from that the factor $(n_{\rm E}/W)^{\delta}$ strongly shifts outwards the location of the critical point and produces a bifurcation in the solution topology.

Table 4: Numerical results for the rotation-CAK model with $\delta =0.02$. Note: the mass loss rate is given in units of $10^{-9}~ M_{\odot} / \rm year$ and the terminal velocity is in km s^-1.

Table 5: Same as Table 4 but for $\delta =0.1$.

	Figure 5: The velocity profile v(r) (in km s-1) as function of $r/R_{\ast }-1$ for the non rotating case (continuous-line) and for the rotational cases: $\Omega =0.25$ dashed-line; $\Omega =0.8$ dotted-line.
Open with DEXTER

   
6.2.1 Bifurcation rotational speed

In order to have an analytical approximation for the value of  $\Omega_{\rm bif}$, we can approximate for the function h(u), Eq. (43):

$$h(u) \simeq -\frac{\delta }{2 }u ,$$ (65)

then the minimum of $\gamma(u)$ is achieved at:

$$u_{\rm bif} \simeq -2\left( \frac{1-\alpha }{\delta A}\right) ^{1/3} ,$$ (66)

and the minimum value of $\gamma(u_{\rm bif})$ is

$$\gamma(u_{\rm bif}) \simeq -2~\left( 1-\alpha \right) ~\Omega^{2}... ...^{-2/3}\right) +3~~\left( 1-\alpha \right) ^{1/3}\delta ^{2/3}A^{2/3}-2~\delta.$$ (67)

From $\gamma(u_{\rm bif}) = 0$, we can obtain the bifurcation value of $\Omega$ given by

$$\Omega_{\rm bif} \simeq \sqrt{\frac{3~~\left( 1-\alpha \right... ...t( 1-\alpha \right) ^{-1/3}\delta ^{1/3}A^{-2/3}\right)}}\cdot$$ (68)

	Figure 6: The curves $r_{\rm min}$ and $r_{\pm }$ are shown with continuous lines as function of $\Omega$. These curves represent the boundaries for the type of the singular point topology. The numerical results for the location of the singular point is also shown in dashed-line. a) is for $\delta =0.02$ and b) for $\delta =0.1$. See text for details.
Open with DEXTER

Figure 6 show the curves $r_{\rm min}$ and $r_{\pm }$ as function of $\Omega$. The intersection point for all the curves (in continuous-line) is at $\Omega=\Omega_{\rm bif}$. Critical points can not be located between the curves r₊ and r_- (filled region). In addition, we show curves for the location of the critical point from numerical calculations (dashed-lines), with the lower boundary condition, $\tau_{\rm Phot}=2/3$.

We can clearly see from this figure, that the position of the critical point is shifted outwards from the stellar surface and the greater is $\delta$ the further is the position of this critical point. It can be inferred from this figure, that the location of the singular point remains almost constant as long as $\Omega ~ \le ~ \Omega_{\rm bif}$, but from values of $\Omega ~ > ~ \Omega_{\rm bif}$ the position of  $r_{\rm crit}$ grows almost linear with $\Omega$.

This behavior of the solution topology can be applied for the winds of Be-Stars. At polar latitudes, i.e., slow rotational speed, the wind behaves as the standard CAK wind, but as the latitude approaches to the equator, the rotational speed is larger than  $\Omega_{\rm bif}$ and the wind is slower and denser. This transition from polar to equatorial latitudes seems to have a similar behavior described by Curé (2004) for the more general rotating m-CAK wind. The study of the influence of this bifurcation in the winds of Be stars will be the scope of a forthcoming article.

   
6.2.2 Phase diagram with ${\delta \neq 0}$

Figure 7 shows the phase diagram p versus $r/R_{\ast }-1$ for $\Omega =0.25$ and $\delta =0.02$ for our test star. The solution topology seems to be similar to the ${\delta =0}$ case. Here the shallow and steep curves are in continuous line and dashed line respectively. The shallow solution is the CAK solution while the steep solution correspond to accretion flows or for radiation driven disk winds (Feldmeier et al. 2002). As we mentioned in previous section, we move slightly from the singular point in the unstable manifold in the directions $\pm(1,\nu_{+})$ and then integrate outwards and inwards, obtaining the solution topology of Fig. 7. Figure 8 shows both shallow and steep solutions, but here we have started from the singular point but with different values of the eigenvalue $\bar{C}$. The almost horizontal dotted curves correspond to $p_{\rm c}$ from Eq. (47) and is not a continuous curve because of the bifurcation in the solution topology. The left-most shallow curve that start at a X-type singular point do not reach stellar surface, but is trapped by the attractor singular point. This result implies that there is a shorter window of locations of singular points, that fixes the eigenvalues $\bar{C}$ and can reach infinity passing through the X-type critical point. This feature is not present when ${\delta =0}$.

	Figure 7: The topology of the rotating-CAK model, p versus $r/R_{\ast }-1$ for $\Omega =0.25$ and $\delta =0.02$. The unique curve that starting at the stellar surface and reaches infinity is the CAK original solution (continuous-line).
Open with DEXTER

	Figure 8: The topology of the rotating-CAK model, p versus $(r/R_{\ast }-1)$ for $\Omega =0.25$. The different curves correspond to different values of the eigenvalues. The left-most curve starts at the "normal'' critical point and is trapped by the attractor. See text for details.
Open with DEXTER

   
7 Conclusions

In this article we have examined the topology of the rotating-CAK wind non-linear differential equation. After the introduction of an additional physically meaningless independent variable ($\tau$) we transform the momentum equation to a set of equations where all the derivatives are explicitly given. This formalism permitted us to linearise the equations in the neighborhood of the critical points. This linearization let us to define a condition for the derivative of the eigenvalue that defines the topology of the critical point, i.e, X-type or attractor.

We have applied our results to the case of a point star (CAK) for a frozen in ionization rotating wind, recovering and generalizing the results of previous studies (Bjorkman 1995; Marlborough & Zamir 1984). The most significant result (with ${\delta =0}$) is that the wind does not assume the maximum mass-loss rate but almost the minimum.

For the more general case, where changes in the wind ionization are taken into account, our analysis shows the existence of a bifurcation in the solution topology, where two critical points exist. The first critical point (closer to the star's surface) is an attractor while the second is the standard CAK critical point. Besides the known fact that the rotational speed shifts the location of the critical point outwards in the wind, the inclusion of the term $(n_{\rm E}/W)^{\delta}$ produces the same effect, reinforcing this displacement.

The bifurcation topology seems to explain the results from Curé (2004) that there exist two regions in the wind of a fast rotating hot star, one where the wind is the one from the standard solution (fast wind) and the other with a new solution that is slower and denser. This result shows us the necessity to perform a topological analysis of the rotating m-CAK wind. This study is currently underway.

Acknowledgements

This work has been possible thanks to the research cooperation agreement UBA/UV and DIUV project 15/2003. M.C. wants to thank the hospitality of the colleges from the mathematics department from the UBA during his stay in Buenos Aires.

Appendix A: Linearization in the neighborhood of a singular point

Considering the system of differential equations given by:

$$\dot{x}= {\rm d}x/{\rm d}t =X\left( x,y\right) ,$$

$$\dot{y}= {\rm d}y/{\rm d}t =Y\left( x,y\right) ,$$ (A.1)

where x and y are the dependent variables and t is the independent variable. The point $\left( x_{\rm c},y_{\rm c}\right)$ is singular (critical) if and only if verifies that:

$$\begin{eqnarray*}X\left( x_{\rm c},y_{\rm c}\right)& =0 \nonumber\\ Y\left( x_{\rm c},y_{\rm c}\right)& =0. \nonumber \end{eqnarray*}$$

In order to understand the topological behavior of the solutions close to a singular point, we expand this system using Taylor series at $\left( x_{\rm c},y_{\rm c}\right)$, obtaining:

$$\begin{eqnarray*}\dot{x} =X\vert _{\left( x_{\rm c},y_{\rm c}\right) }+X_{x}\ver... ...,y_{\rm c}\right) }\left( y-y_{\rm c}\right) +o\left(x,y\right). \end{eqnarray*}$$

Neglecting superior order terms and using that $\left( x_{\rm c},y_{\rm c}\right)$is a singular point, we have, in matrix form:

$$\left( \begin{array}[c]{c}% \dot{x}\\ \dot{y}% \end{array}\r... ...n{array}[c]{c}% x-x_{\rm c}\\ y-y_{\rm c}% \end{array}\right)$$ (A.2)

where B is the Jacobian matrix at $\left( x_{\rm c},y_{\rm c}\right)$, i.e.,

$$B=\left. \left( \begin{array}[c]{cc}% X_{x} & X_{y}\\ Y_{x} ... ...ray}\right) \right\vert _{\left( x_{\rm c},y_{\rm c}\right) }.$$ (A.3)

For a 3-dimensional system case:

 

$$\dot{x} =X\left( x,y,z\right) ,$$

$$\dot{y} =Y\left( x,y,z\right) ,$$

$$\dot{z} =Z\left( x,y,z\right),$$ (A.4)

with a constraint $\phi\left( x,y,z\right)$ satisfying:

$$\phi_{x}X+\phi_{y}Y+\phi_{z}Z = 0 ,$$ (A.5)

the surface defined by:

$$S=\left\{ \left( x,y,z\right) :\phi\left( x,y,z\right) =0\right\} ,$$ (A.6)

is invariant for the differential equation system evolution.

From the Eq. (A.6), we can solve for z (using the implicit function theorem) obtaining $z=z\left( x,y\right)$ and reduce the system (A.4) to:

 

$$\dot{x} =X\left( x,y,z\left( x,y\right) \right) ,$$

$$\dot{y} =Y\left( x,y,z\left( x,y\right) \right) .$$ (A.7)

If $\left( x_{\rm c},y_{\rm c},z_{\rm c}\right)$ is a singular point of Eq. (A.7), the Jacobian matrix B is given by:

$$B=\left. \left( \begin{array}[c]{cc}% X_{x}+z_{x}X_{z} & X_{y... ...y}\right) \right\vert _{\left( x_{\rm c},y_{\rm c}\right) } .%$$ (A.8)

Since $z_{x}=-\phi_{x}/\phi_{z}$ and $z_{y}=-\phi_{y}/\phi_{z}$, we have

$$B=\frac{1}{\phi_{z}}\left. \left( \begin{array}[c]{cc}% \phi_... ...y}\right) \right\vert _{\left( x_{\rm c},y_{\rm c}\right) } .%$$ (A.9)

For $\phi=F$, X=U and Y=P, we obtain Eq. (30).

References

• Abbott, D. C. 1982, ApJ, 259, 282 [NASA ADS] [CrossRef] (In the text)
• Abbott, D. C., & Lucy, L. B. 1985, ApJ, 288, 679 [NASA ADS] [CrossRef] (In the text)
• Abt, H. A., Levato, H., & Grosso, M. 2002, ApJ, 573, 359 [NASA ADS] [CrossRef] (In the text)
• Amann, H. 1990, Ordinary differential equations: an introduction to nonlinear analysis. Gruyter studies in mathematics 13 (Berlin: Ed. Walter de Gruyter)
• Bjorkman, J. E. 1995, ApJ, 453, 369 [NASA ADS] [CrossRef] (In the text)
• Castor, J. I., Abbott, D. C., & Klein, R. 1975, ApJ, 195, 157 [NASA ADS] [CrossRef] (In the text)
• Castor, J. I. 1979, Mass Loss and Evolution of O-Type Stars, ed. P. S. Conti, & C. W. de Loore (Dordrecht: Reidel), IAU Symp., 83, 175 (In the text)
• Chauville, J., Zorec, J., Ballereau, D., et al. 2001, A&A, 378, 861 [EDP Sciences] [NASA ADS] [CrossRef] (In the text)
• Curé, M. 2004, ApJ, 614, 929 [NASA ADS] [CrossRef] (In the text)
• Friend, D., & Abbott, D. C. 1986, ApJ, 311, 701 [NASA ADS] [CrossRef] (In the text)
• Feldmeier, A., Shlosman, I., & Hamann, W.-R. 2002, A&A, 566, 392 [NASA ADS] (In the text)
• Gayley, K. G., Owocki, S. P., & Cranmer, S. R. 1995, ApJ, 442, 296 [NASA ADS] [CrossRef] (In the text)
• Hutchings, J. B., Nemec, J. M., & Cassidy, J. 1979, PASP, 91, 313 [NASA ADS] [CrossRef] (In the text)
• Kudritzki, R. P., Pauldrach, A., Puls, J., & Abbott, D. C. 1989, A&A, 219, 205 [NASA ADS] (In the text)
• Kudritzki, R. P., Puls, J., Lennon, D. J., et al. 1999, A&A, 350, 970 [NASA ADS] (In the text)
• Kudritzki, R. P., & Puls, J. 2000, ARA&A, 38, 613 [NASA ADS] (In the text)
• Kudritzki, R. P. 2002, ApJ, 577, 389 [NASA ADS] [CrossRef] (In the text)
• Lamers, H. J. G. L. M., & Leitherer, C. 1993, ApJ, 412, 771 [NASA ADS] [CrossRef] (In the text)
• Lamers, H. J. G. L. M., & Cassinelli, J. 1999, Introduction to stellar winds (Cambridge: Cambridge University Press) (In the text)
• Lucy, L. B., & Solomon, P. M. 1970, ApJ, 159, 879 [NASA ADS] [CrossRef] (In the text)
• Maeder, A. 2001, A&A, 373, 122 [EDP Sciences] [NASA ADS] [CrossRef] (In the text)
• Marlborough, J. M., & Zamir, M. 1984, ApJ, 276, 706 [NASA ADS] [CrossRef] (In the text)
• Owocki, S. P., & ud-Doula, A. 2004, ApJ, 600, 1004 [NASA ADS] [CrossRef] (In the text)
• Pauldrach, A., Puls, J., & Kudritzki, R. P. 1986, A&A, 164, 86 [NASA ADS] (In the text)
• Puls, J. 1987, A&A, 184, 227 [NASA ADS] (In the text)
• Puls, J., Kudritzki, R.-P., Herrero, A., et al. 1996, A&A, 305, 171 [NASA ADS] (In the text)
• Slettebak, A. 1976, in Be and Shell Stars, ed. Arne Slettebak (Dordrecht: D. Reidel Pub. Co.), IAU Symp., 70, 123 (In the text)
• Vink, J. S., de Koter, A., & Lamers, H. J. G. L. M. 2000, A&A, 362, 295 [NASA ADS] (In the text)