© ESO, 2014

1. Introduction

In recent years considerable progress has been made in our theoretical modelling of rotating massive stars (Maeder & Meynet 2012; Langer 2012) and our basic understanding of spherical radiation-driven winds (Puls et al. 2008). However, in order to get a grasp on the non-spherical 2D outflows of rotating massive stars, involving B[e] supergiants (Zickgraf et al. 1985), classical Be stars (Porter & Rivinius 2003) as well as luminous blue variable (LBV) outflows (Groh et al. 2006), it is paramount to combine the intrinsically 2D nature of rotation and mass loss (e.g. Lovekin 2011; Espinosa & Rieutord 2013). This is not only required for a basic understanding of massive star evolution, but also for linking the oftentimes non-spherical supernova (SN) data with their progenitors (e.g. Maund et al. 2007; Hoffman et al. 2008).

We need to develop 2D wind models in order to obtain a physical understanding of how rotation might affect both the strength and latitudinal dependence of their outflows. In turn winds may be able to remove significant quantities of angular momentum, potentially down to masses as low as 10–15 M_⊙ (Vink et al. 2010). Whether the mass loss originates from the pole or the equator remains currently unknown. Yet, is of paramount importance for understanding whether rapid rotation is maintained or leads to stellar spin loss (Meynet & Maeder 2007), highly relevant for our understanding of the progenitor evolution of long-duration gamma ray bursts (GRBs).

Previous models of the winds from rotating stars have mostly been 1D models of the equatorial flow versus the polar flow, although one 2D numerical calculation has been performed by Poe (1987). The 1D model of Friend & Abbott (1986, hereafter FA) concerned the influence of stellar rotation on the hydrodynamics of a stellar wind, involving a solution of the fluid equations in the equatorial plane, which included centrifugal forces. They used a form of the radiation force after Castor et al. (1975, hereafter CAK), but corrected for the finite angular size of the stellar “disk” (see also Pauldrach et al. 1986).

Bjorkman & Cassinelli (1993, hereafter BC) provided an analytical approximation for the axisymmetric 2D supersonic solution (i.e. for the velocity field and density distribution) of a rotating radiation-driven wind, obtained from the FA 1D model of the equatorial flow. In order to find the streamline trajectories, they rotated the FA 1D solution (in the equatorial plane) up to the initial co-latitude θ₀ of the streamline at the stellar surface, adjusting the equatorial rotation velocity of the central star v_rot by v_rotsinθ₀. The supersonic solutions obtained this way provided the velocity and density as a function of θ₀ and radius r, i.e. in a non-explicit form: given a location (r,θ), one needs to find θ₀ of the streamline that passes through that location, by iteratively solving additional equations. As a result, they explained how rotation can lead to the production of a dense equatorial disk around, e.g. Be stars, by means of their wind-compressed disk (WCD) model.

Refinements of the BC model have been also made for simulating the density structure of rotating O-star winds (e.g. Petrenz & Puls 1996). Owocki et al. (1996) showed that the inclusion of non-radial line forces leads to a small polarwards v_θ component, which may inhibit disk formation in Be stars. Moreover, Maeder (1999) showed that gravity darkening as a result of the Von Zeipel (1924) effect will generally lead to a polar wind. It should be noted that the hydrodynamical wind models of Owocki and colleagues that lead to polar winds employ an approximated line driving formula from CAK for 1D.

A generally prolate wind structure was however confirmed by the sectorial 1.5D wind models of Pelupessy et al. (2000) that employed 1D detailed Monte Carlo line acceleration computations of Vink et al. (1999). In their models for B[e] supergiants, Pelupessy et al. also showed that when models are in close proximity to the bi-stability jump it is possible to overcome the polar enhancement due to the Von Zeipel effect, and drive equatorial enhancements, as originally suggested by Lamers & Pauldrach (1991) for Be and B[e] supergiant disk formation. Cure et al. (2005) and Madura et al. (2007) derived 1D hydrodynamical models for very rapid rotators (above 75% of the critical rate) finding a slow solution to the classical CAK theory, which may enable disk formation in Be stars and B[e] supergiants, when accounting for the wind bi-stability effects of Vink et al. (1999). However, again, these models employ a simplified treatment of the line acceleration due to the 2 (or 3) parameter power-law approximation due to CAK. What is eventually required in order to resolve the intricate problem of stellar rotation with mass loss are 3D Monte Carlo radiative transfer calculations in combination with a full hydrodynamic solution. Most published models have necessarily made significant assumptions with respect to either the line-force calculations or the wind hydrodynamics.

We suggested a new parametrisation of the line acceleration (Müller & Vink 2008, hereafter Paper I), expressing it as a function of radius rather than of the velocity gradient as in CAK theory. The implementation of this formalism allowed for local dynamical consistency as we were able to determine the momentum transfer at each location in the wind through the use of Monte Carlo simulations. In Muijres et al. (2012), we tested our hydrodynamic wind solutions and velocity laws by additional explicit numerical integrations of our fluid equation, also accounting for a temperature stratification. These results were in excellent agreement with both our full and our approximated solutions from Paper I. We here build on those results, now deriving analytical solutions for the 2D case concerning both the velocity and density structure in an axisymmetric mass outflow (or inflow) scenario. Furthermore, we extend our iterative method from Paper I for the simultaneous solution of the mass-loss rate and velocity field to the 2D case of a rotating non-spherical stellar wind.

We obtain the velocity field fully analytically without any previous fits to numerical solutions of the fluid equation of FA, if we neglect the polar velocity v_θ in our model. We are justified in doing so as long as the stellar rotation speeds are well below the critical value where disks are formed, that is, for example for O-star winds, where non-spherical outflows have only been detected in a small (~20%) minority only involving “special” O-type sub-groups (Oe, Onfp) from linear spectropolarimetry observations (Harries et al. 2002; Vink et al. 2009).

However, the non-spherical problem for the flow in the equatorial plane including the case of the outflow at the pole (where v_θ ≡ 0 for symmetrical reasons) are as well fully analytically solved by our improved 2D wind model without any restrictions to the velocity components and the rotational speed of the central star. For the specific case of a non-spherical radiation-driven wind, we do not rely on the CAK expression for the radiation force, rather we describe the line acceleration as a function of stellar radius g(r,θ) at a given constant co-latitude θ. In addition, the critical point of our stellar wind is the sonic point (depending on latitude) and not the CAK critical point. The calculation of g(r,θ) is performed through Monte Carlo simulations accounting for multi-line transfer, and the wind parameters are solved simultaneously – in an iterative way – for each latitude of interest.

The set-up of the paper is as follows. In Sects. 2.2–2.5, the hydrodynamic equations for a non-spherical axisymmetric steady flow are introduced including a derivation of the mathematical description of the radiative line acceleration as a function of radius for the case of a rotating radiation-driven wind. The process for obtaining the fully analytical 2D solutions is described and discussed in Sect. 2.6. Here, the velocity field for the entire family of solutions is provided in an explicit general expression from which the solutions for a rotating radiation-driven wind or mass accretion flux (e.g. collapsing protostellar cloud) follow as unique trans-sonic solutions through the critical point. Moreover, an approximate analytical solution for the supersonic flow is presented. In Sect. 3, we describe our numerical computation obtaining the radiative acceleration in our stellar wind models. Furthermore, our iterative method for the determination of the consistent solution for the mass-loss rate in case of a spherical wind is being extended and applied to the wind from a rotating star. In Sect. 4, we present the application of our models to a differentially rotating stellar wind from a typical 40 M_⊙ O5–V-star, including the effects of oblateness and gravity darkening, and from a rotating 60 M_⊙ O-giant star. We discuss the results, before we summarise and discuss our findings in Sect. 5.

2. Radiation hydrodynamics of rotating and expanding or collapsing systems

2.1. The velocity field

The velocity field of the differentially rotating system at location r = (r,   θ,   φ) can generally be described by its spherical components ${v}\hspace{0.17em}\mathrm{\left(}{r}\mathrm{\right)}\mathrm{=}{\mathit{v}}_{\mathit{r}}\hspace{0.17em}\mathrm{\left(}{r}\mathrm{\right)}\hspace{0.17em}{{e}}_{\mathit{r}}\hspace{0.17em}\mathrm{\left(}{r}\mathrm{\right)}\mathrm{+}{\mathit{v}}_{\mathit{\phi }}\hspace{0.17em}\mathrm{\left(}{r}\mathrm{\right)}\hspace{0.17em}{{e}}_{{\phi }}\hspace{0.17em}\mathrm{\left(}{r}\mathrm{\right)}\mathrm{+}{\mathit{v}}_{\mathit{\theta }}\hspace{0.17em}\mathrm{\left(}{r}\mathrm{\right)}\hspace{0.17em}{{e}}_{{\theta }}\hspace{0.17em}\mathrm{\left(}{r}\mathrm{\right)}$(1)with the unit vectors e_r, e_θ, e_φ, in radial, polar and azimuthal direction, respectively, where the following two presentations of the unit vectors ${{e}}_{\mathit{r}}\hspace{0.17em}\mathrm{\left(}{r}\mathrm{\right)}\mathrm{=}\left(\begin{array}{c}\\ \sin \mathit{\theta }\hspace{0.17em}\cos \mathit{\phi }\\ \sin \mathit{\theta }\hspace{0.17em}\sin \mathit{\phi }\\ \cos \mathit{\theta }\end{array}\right) \mathrm{and} {{e}}_{{\phi }}\hspace{0.17em}\mathrm{\left(}{r}\mathrm{\right)}\mathrm{=}\left(\begin{array}{c}\\ \mathrm{-}\sin \mathit{\phi }\\ \cos \mathit{\phi }\\ \mathrm{0}\end{array}\right)$(2)will be useful¹.

2.2. Basic equations of hydrodynamics

Considering the particular chosen system as a non-viscous², i.e., ideal fluid, the momentum equation $\mathit{\rho }\hspace{0.17em}\frac{\mathrm{D}\hspace{0.17em}{v}}{\mathrm{D}\hspace{0.17em}\mathit{t}}\mathrm{=}{f}\mathrm{-}{\nabla }\hspace{0.17em}\mathit{p}$(3)is valid (see, e.g., Mihalas & Weibel Mihalas 1984), where D/D t is the covariant Lagrangean or co-moving time derivative in the fluid-frame of a material element and v is its velocity, f is the total external body force per volume acting on a mass element of fluid, and ∇   p is the divergence term of a diagonal isotropic stress tensor ∇·T, in which T =  −p   I and p is the hydrostatic pressure.

One also needs to consider the equation of continuity $\frac{\mathit{\partial \rho }}{\mathit{\partial }\mathit{t}}\mathrm{+}{\nabla }\mathrm{·}\left(\mathit{\rho }\hspace{0.17em}{v}\right)\mathrm{=}\mathrm{0}\mathit{,}$(4)with the covariant divergence ∇·v of the velocity vector.

2.3. Simplifying assumptions

Besides the assumption of an inviscid flow and to account for the axisymmetric and stationary problem, we make the further following simplifying assumptions to solve the hydrodynamic equations analytically:

• 1.

  The stellar wind (flow) is first assumed to be isothermal to derivethe analytical expressions for the hydrodynamic solutions. In thiscase, the equation of state$\mathit{p}\mathrm{=}{\mathit{a}}^{\mathrm{2}}\mathit{\rho }$(5)is valid, where a is the isothermal speed of sound and ρ is the density of the wind (system). This assumption, however, will be relaxed later by applying our iterative method (described in Sect. 3), to compute the mass-loss rates and the parameters in our analytical wind solutions consistently with the radiation field and ionisation/excitation state of the gas of an outflowing stellar model atmosphere in non local thermodynamic equilibrium, assuming radiative equilibrium, by the use of the ISA-Wind code allowing a temperature stratification³. Additionally, in Muijres et al. (2012), we verified our hydrodynamic solutions, especially for the case of a spherical wind without rotation, by explicit numerical integrations of our fluid equation also accounting for a temperature distribution.

• 2.

  We assume a stationary axisymmetric flow, since we are interested in a rotating star (system) where rotational effects dominate the flow, i.e. we set $\frac{\mathit{\partial }}{\mathit{\partial }\mathit{t}}\mathrm{\equiv }\mathrm{0}\mathit{,} \frac{\mathit{\partial }}{\mathit{\partial }\mathit{\phi }}\mathrm{\equiv }\mathrm{0}\mathit{,}$(6)and therefore disregard shocks as well. Furthermore, we exclude the presence of clumps.

• 3.

  For symmetrical reasons, the polar component of the flow velocity should be ${\mathit{v}}_{\mathit{\theta }}\mathrm{=}\mathrm{0}\mathit{,}$(7)as well as the polar and azimuthal force components, ${\mathit{f}}_{\mathit{\theta }}\mathrm{=}\mathrm{0}$(8)and ${\mathit{f}}_{\mathit{\phi }}\mathrm{=}\mathrm{0}\mathit{,}$(9)respectively, in the equatorial plane and at the pole for an axisymmetric flow. At the pole, the hydrodynamic solutions must pass over into those for the spherical case without rotation. Eqs. (7) to (9) allow us then the seperation of the radial motion for any individual latitude. However, for intermediate latitudes, these approximations do generally not hold: for instance, neglecting the meridional velocity there, restricts the application of our model only to those cases where matter exchange between layers of different latitude (and therefore also the occurrance of a friction between) can be neglected. We therefore apply our model, in the case of stellar winds, only to rotating O-stars with rotational speeds below the critical value where disks are formed. Moreover, the additional involvement of the distortion of a central star due to its rotation (and consequently the effect of gravity darkening as well) in the course of our investigations, may result in a non-vanishing radiative force f_θ in θ-direction at all mid-latitudes. Therefore, for winds from those rotating O-stars, we employ our formalism and method (in Sect. 4) exlusively only on the equatorial plane and the pole, to compute the more accurate wind parameters and mass-loss rates there, where the above assumptions (Eqs. (7)–(9)) are satisfied best⁴. Then, at latitudes between the pole and equator, all corresponding hydrodynamic quantities adopt values which lie in-between these two extreme values (constraints).

• 4.

  In the case of a wind from a luminous early-type star, the wind is primarily driven by continuum plus line radiation forces, where the radial acceleration on a mass element is $\frac{{\mathit{f}}_{\mathit{r}}}{\mathit{\rho }}\mathrm{=}\mathrm{-}\frac{\mathit{G}\hspace{0.17em}\mathit{M}}{{\mathit{r}}^{\mathrm{2}}}\hspace{0.17em}\mathrm{(}\mathrm{1}\mathrm{-}{\mathrm{\Gamma }}^{\mathrm{)}}\mathrm{+}{\mathit{g}}_{\mathrm{rad}}^{\mathrm{line}}$(10)with $\mathrm{\Gamma :}\mathrm{=}\frac{{\mathit{g}}_{\mathrm{rad}}^{\mathrm{cont}}}{\mathit{g}}\mathit{,}$(11)the force ratio between the radiative acceleration ${\mathit{g}}_{\mathrm{rad}}^{\mathrm{cont}}$ due to the continuum opacity (dominated by electron scattering) and the inward acceleration of gravity g. Γ is supposed to be independent of radius r, may however vary with polar angle θ, and ${\mathit{g}}_{\mathrm{rad}}^{\mathrm{line}}\hspace{0.17em}\mathrm{\left(}\mathit{r,\theta }\mathrm{\right)}$ is the outward radiative acceleration due to spectral lines. M is the mass of the central star.

2.4. Simplified hydrodynamic equations

2.4.1. Wind density and mass-loss rate

If we use the covariant derivative (see Mihalas & Weibel Mihalas 1984), for spherical coordinates and apply assumptions 2 and 3 to the equation of continuity (4), we find $\frac{\mathit{\partial }}{\mathit{\partial }\mathit{r}}\hspace{0.17em}\mathrm{(}{\mathit{r}}^{\mathrm{2}}\hspace{0.17em}\mathit{\rho }\hspace{0.17em}{\mathit{v}}_{\mathit{r}}\mathrm{)}\mathrm{=}\mathrm{0}$(12)for a two-dimensional axisymmetric and steady flow.

This equation looks the same as that one, we would get for a one-dimensional, spherically symmetric and steady flow, however, this expression is quite more general and only fulfilled for a given constant polar angle θ. But similar to a spherically symmetric flow, the integration of Eq. (12) yields $\mathrm{4}\hspace{0.17em}\mathit{\pi }\hspace{0.17em}{\mathit{r}}^{\mathrm{2}}\hspace{0.17em}\mathit{\rho }\hspace{0.17em}\mathrm{\left(}\mathit{r,\theta }\mathrm{\right)}\hspace{0.17em}{\mathit{v}}_{\mathit{r}}\hspace{0.17em}\mathrm{\left(}\mathit{r,\theta }\mathrm{\right)}\mathrm{=}\mathrm{const}\mathrm{.}\mathrm{=}\mathrm{:}\mathit{Ṁ}\hspace{0.17em}\mathrm{\left(}\mathit{\theta }\mathrm{\right)}\mathit{,}$(13)in which here Ṁ(θ) is not the total mass flux through a spherical shell surrounding the star, but its mass flux (or mass-loss rate) at co-latitude θ through the star surface multiplied by 4   π   R² (cf. definition in BC) that is conserved for each angle θ. The value of Ṁ   (θ) can later be determined by the given values of velocity v_r   (R,θ) and density ρ   (R) at the stellar radius R (or at the surface of an inner core).

The total mass-loss rate must then be $\mathrm{ℳ̇}\mathrm{\equiv }{\mathrm{\int }}_{\mathrm{△}\hspace{0.17em}\mathit{\pi }}\frac{\mathrm{ℳ̇}\hspace{0.17em}\mathrm{\left(}\mathit{\theta }\mathrm{\right)}}{\mathrm{△}\hspace{0.17em}\mathit{\pi }}\hspace{0.17em}\mathrm{d}\mathrm{\otimes }\mathit{,}$(14)integrated over the solid angle Ω. Note that in case of a collapsing system (e.g. protostellar cloud) this mass-loss rate is negative, because the inner core gains mass.

Finally, by Eq. (13), we obtain the 2D density distribution $\mathit{\rho }\hspace{0.17em}\mathrm{\left(}\mathit{r,\theta }\mathrm{\right)}\mathrm{=}\frac{\mathit{Ṁ}\hspace{0.17em}\mathrm{\left(}\mathit{\theta }\mathrm{\right)}}{\mathrm{4}\hspace{0.17em}\mathit{\pi }\hspace{0.17em}{\mathit{r}}^{\mathrm{2}}\hspace{0.17em}{\mathit{v}}_{\mathit{r}}\hspace{0.17em}\mathrm{\left(}\mathit{r,\theta }\mathrm{\right)}}\mathrm{=}\frac{\mathit{F}\hspace{0.17em}\mathrm{\left(}\mathit{\theta }\mathrm{\right)}}{{\mathit{r̂}}^{\mathrm{2}}\hspace{0.17em}{\mathit{v}}_{\mathit{r}}\hspace{0.17em}\mathrm{\left(}\mathit{r̂}\mathit{,\theta }\mathrm{\right)}}$(15)at location (r,θ), with the defined flux F = Ṁ   (θ)/4   π   R² through the star’s surface at radius R and the dimensionless radius .

Please note that all formulae derived in this Sect. 2 are expressed in terms of referring to the radius R, which is (throughout this section) the stellar (i.e. photospheric) radius of the central star (or the inner core radius of any central object, respectively). However, all formulae can generally also be applied with respect to the reference radius R to be the inner boundary radius R_in, from where the numerical computations of the (stellar wind) model start.

2.4.2. The azimuthal velocity component and the correction for oblateness

By using the correct contravariant components of acceleration (D   v_i/D   t) in Eq. (3), for spherical coordinates, and replacing them by their equivalent physical components (see again, e.g., Mihalas & Weibel Mihalas 1984), and applying assumptions 1–4, we obtain the simplified r- and φ-component of the momentum equation

$\begin{array}{ccc}& & {\mathit{v}}_{\mathit{r}}\hspace{0.17em}\frac{\mathit{\partial }}{\mathit{\partial }\mathit{r}}\hspace{0.17em}{\mathit{v}}_{\mathit{r}}\mathrm{-}\frac{{\mathit{v}}_{\mathit{\phi }}^{\mathrm{2}}}{\mathit{r}}\mathrm{=}\frac{{\mathit{f}}_{\mathit{r}}}{\mathit{\rho }}\mathrm{-}\frac{{\mathit{a}}^{\mathrm{2}}}{\mathit{\rho }}\hspace{0.17em}\frac{\mathit{\partial }\mathit{\rho }}{\mathit{\partial }\mathit{r}}\mathit{,}\\ & & {\mathit{v}}_{\mathit{r}}\hspace{0.17em}\frac{\mathit{\partial }}{\mathit{\partial }\mathit{r}}\hspace{0.17em}{\mathit{v}}_{\mathit{\phi }}\mathrm{+}\frac{{\mathit{v}}_{\mathit{r}}\hspace{0.17em}{\mathit{v}}_{\mathit{\phi }}}{\mathit{r}}\mathrm{=}\mathrm{0}\mathit{,} \mathrm{orequivalently}\\ & & \frac{{\mathit{v}}_{\mathit{r}}}{\mathit{r}}\hspace{0.17em}\frac{\mathit{\partial }}{\mathit{\partial }\mathit{r}}\hspace{0.17em}\mathrm{(}\mathit{r}\hspace{0.17em}{\mathit{v}}_{\mathit{\phi }}\mathrm{)}\mathrm{=}\mathrm{0}\mathit{,}\end{array}$with the external radial force per unit mass f_r, i.e. the radial acceleration on the mass element in Eq. (10), in case of a stellar wind.

The φ-component of the momentum Eq. (17) is nothing more than the conservation of angular momentum per unit mass $\begin{array}{ccc}\mathit{r}\hspace{0.17em}{\mathit{v}}_{\mathit{\phi }}\hspace{0.17em}\mathrm{\left(}\mathit{r,}\hspace{0.17em}\mathit{\theta }\mathrm{\right)}\mathrm{=}\mathrm{const}\mathrm{.}\mathit{,}& & \end{array}$as one would expect for external central forces and axisymmetry, what we have, of course, supposed before. Then, the last differential Eq. (17) can be solved (i.e. integrated) immediately and separately from Eq. (16), to obtain the unknown velocity component v_φ$\begin{array}{ccc}\mathit{r}\hspace{0.17em}{\mathit{v}}_{\mathit{\phi }}\hspace{0.17em}\mathrm{\left(}\mathit{r,}\hspace{0.17em}\mathit{\theta }\mathrm{\right)}\stackrel{\mathrm{!}}{\mathrm{=}}\mathit{R}\hspace{0.17em}{\mathit{v}}_{\mathit{\phi }}\hspace{0.17em}\mathrm{\left(}\mathit{R,}\hspace{0.17em}\mathit{\theta }\mathrm{\right)}& & \end{array}$by choosing an adequate boundary (initial) condition ${\mathit{v}}_{\mathrm{rot}}\hspace{0.17em}\mathrm{\left(}\mathit{\theta }\mathrm{\right)}\mathrm{:}\mathrm{=}{\mathit{v}}_{\mathit{\phi }}\hspace{0.17em}\mathrm{(}\mathit{R,}\hspace{0.17em}\mathit{\theta }\mathrm{)}\mathit{,}$(18)i.e. rotational speed of the star (inner core) surface at co-latitude θ.

Hence, the φ-component of the velocity of a particle at distance in its orbit, originating from the stellar surface at co-latitude θ and ejected with v_rot   (θ), remaining on the cone surface of constant angle θ, becomes ${\mathit{v}}_{\mathit{\phi }}\hspace{0.17em}\mathrm{\left(}\mathit{r̂}\mathit{,\theta }\mathrm{\right)}\mathrm{=}\frac{\mathrm{1}}{\mathit{r̂}}\hspace{0.17em}{\mathit{v}}_{\mathrm{rot}}\hspace{0.17em}\mathrm{\left(}\mathit{\theta }\mathrm{\right)}\mathit{.}$(19)Assuming that the central star surface (or inner core surface) behaves like a rotating rigid sphere (R = const.), the rotational velocity at co-latitude θ would then be described by ${\mathit{v}}_{\mathrm{rot}}\hspace{0.17em}\mathrm{\left(}\mathit{\theta }\mathrm{\right)}\mathrm{=}{\mathit{V}}_{\mathrm{rot}}\hspace{0.17em}\sin \mathit{\theta }$(20)with the equatorial rotation speed V_rot.

However, due to rapid rotation, the central star (core) can become oblate from the centrifugal forces and can then be described as a rotating rigid ellipsoid (R = R(θ)) with rotational velocity ${\mathit{v}}_{\mathrm{rot}}\hspace{0.17em}\mathrm{\left(}\mathit{\theta }\mathrm{\right)}\mathrm{=}\frac{\mathit{R}\mathrm{\left(}\mathit{\theta }\mathrm{\right)}}{{\mathit{R}}_{\mathrm{eq}}}\hspace{0.17em}{\mathit{V}}_{\mathrm{rot}}\hspace{0.17em}\sin \mathit{\theta }$(21)at co-latitude θ, where R_eq is the radius R(θ = π/2) at the equator.

Here, the stellar (core) radius R(θ), depending on latitude and rotation speed is given by (Cranmer & Owocki 1995, hereafter CO; Petrenz & Puls 1996, hereafter PP) $\mathit{R}\mathrm{\left(}\mathit{\theta }\mathrm{\right)}\mathrm{=}\frac{\mathrm{3}{\mathit{R}}_{\mathrm{p}}}{\mathrm{\Omega }\sin \mathit{\theta }}\cos \left(\frac{\mathit{\pi }\mathrm{+}\arccos \left(\mathrm{\Omega }\sin \mathit{\theta }\right)}{\mathrm{3}}\right)$(22)where R_p is the polar radius and assumed to be independent of the rotational velocity, i.e. used as stellar input parameter, and Ω is the normalised stellar angular velocity and defined by $\mathrm{\Omega }\mathrm{=}\frac{\mathit{\omega }}{{\mathit{\omega }}_{\mathrm{crit}}}\mathrm{=}\frac{\mathrm{1}}{{\mathit{\omega }}_{\mathrm{crit}}}\frac{{\mathit{V}}_{\mathrm{rot}}}{{\mathit{R}}_{\mathrm{eq}}}$(23)with the critical angular velocity ${\mathit{\omega }}_{\mathrm{crit}}\mathrm{=}\sqrt{\frac{\mathrm{8}\hspace{0.17em}\mathit{GM}\mathrm{(}\mathrm{1}\mathrm{-}{\mathrm{\Gamma }}^{\mathrm{)}}}{\mathrm{27}\hspace{0.17em}{\mathit{R}}_{\mathrm{p}}^{\mathrm{3}}}}$(24)and equatorial radius ${\mathit{R}}_{\mathrm{eq}}\mathrm{=}\frac{{\mathit{R}}_{\mathrm{p}}}{\mathrm{1}\mathrm{-}{\mathit{V}}_{\mathrm{rot}}^{\mathrm{2}}{\mathit{R}}_{\mathrm{p}}\mathit{/}\mathrm{(}\mathrm{2}\hspace{0.17em}\mathit{GM}{\mathrm{(}\mathrm{1}\mathrm{-}{\mathrm{\Gamma }}^{\mathrm{)}}}^{\mathrm{)}}}\mathrm{·}$(25)Note that the above definition of ω_crit differs from the break-up velocity ω_{crit,spherical} = v_crit/R_p introduced in the following Sect. 2.4.4, where the rotational distortion of the stellar surface is neglected.

2.4.3. The effect of gravity darkening

If one considers the distortion of the central star due to its rotation, one also has to account for the effect of gravity darkening caused by its oblateness. The early work of von Zeipel (1924) for distorted stars states that the radiative flux F   (θ) emerging from the surface at co-latitude θ is proportional to the local effective gravity $\mathit{F}\hspace{0.17em}\mathrm{\left(}\mathit{\theta }\mathrm{\right)}\mathrm{\propto }{\mathit{g}}_{\mathrm{\perp }}\hspace{0.17em}\mathrm{\left(}\mathit{\theta }\mathrm{\right)}\mathit{.}$(26)Therein the normal component of gravity is (see Collins 1965; or PP) $\begin{array}{ccc}{\mathit{g}}_{\mathrm{\perp }}\hspace{0.17em}\mathrm{\left(}\mathrm{\Omega }\mathrm{\right(}{\mathit{V}}_{\mathrm{rot}}\mathrm{\left)}\mathit{,\theta }\mathrm{\right)}& \mathrm{=}& \frac{\mathit{GM}}{{\mathit{R}}_{\mathrm{p}}^{\mathrm{2}}}\hspace{0.17em}\frac{\mathrm{8}}{\mathrm{27}}[{\left(\frac{\mathrm{27}}{\mathrm{8}}\hspace{0.17em}{\left(\frac{{\mathit{R}}_{\mathrm{p}}}{\mathit{R}\mathrm{\left(}\mathit{\theta }\mathrm{\right)}}\right)}^{\mathrm{2}}\mathrm{-}\frac{\mathit{R}\mathrm{\left(}\mathit{\theta }\mathrm{\right)}}{{\mathit{R}}_{\mathrm{p}}}\hspace{0.17em}{\mathrm{\Omega }}^{\mathrm{2}}\hspace{0.17em}{\sin }^{\mathrm{2}}\mathit{\theta }\right)}^{\mathrm{2}}\\ & & \mathrm{+}{{\mathrm{\Omega }}^{\mathrm{4}}\hspace{0.17em}{\left(\frac{\mathit{R}\mathrm{\left(}\mathit{\theta }\mathrm{\right)}}{{\mathit{R}}_{\mathrm{p}}}\right)}^{\mathrm{2}}\hspace{0.17em}{\sin }^{\mathrm{2}}\mathit{\theta }\hspace{0.17em}{\cos }^{\mathrm{2}}\mathit{\theta }]}^{\mathrm{1}\mathit{/}\mathrm{2}}\mathit{,}\end{array}$(27)with the stellar radius R(θ) in Eq. (22), and the proportional constant $\mathit{C}\hspace{0.17em}\mathrm{\left(}\mathrm{\Omega }\mathrm{\right)}\mathrm{=}\frac{\mathit{L}}{\mathrm{\Sigma }\hspace{0.17em}\mathrm{\left(}\mathrm{\Omega }\mathrm{\right)}}$(28)is given by the constraint that the surface-integrated flux is equal to the total luminosity L of the (oblate) star, where $\begin{array}{ccc}\mathrm{\Sigma }\hspace{0.17em}\mathrm{\left(}\mathrm{\Omega }\mathrm{\right)}& \mathrm{:}\mathrm{=}& 􏽉\mathit{A}{\mathit{g}}_{\mathrm{\perp }}\hspace{0.17em}\mathrm{\left(}\mathrm{\Omega }\mathit{,\theta }\mathrm{\right)}\hspace{0.17em}\mathrm{d}\mathit{A}\\ & \mathrm{\approx }& \mathrm{4}\mathit{\pi }\hspace{0.17em}\mathit{G}\hspace{0.17em}\mathit{M}\mathrm{(}\mathrm{1}\mathrm{-}\mathrm{0.19696}\hspace{0.17em}{\mathrm{\Omega }}^{\mathrm{2}}\mathrm{-}\mathrm{0.094292}\hspace{0.17em}{\mathrm{\Omega }}^{\mathrm{4}}\\ & & \hspace{0.17em}\mathrm{+}\hspace{0.17em}\mathrm{0.33812}\hspace{0.17em}{\mathrm{\Omega }}^{\mathrm{6}}\mathrm{-}\mathrm{1.3066}\hspace{0.17em}{\mathrm{\Omega }}^{\mathrm{8}}\mathrm{+}\mathrm{1.8286}\hspace{0.17em}{\mathrm{\Omega }}^{\mathrm{10}}\\ & & \hspace{0.17em}\mathrm{-}\hspace{0.17em}\mathrm{0.92714}\hspace{0.17em}{\mathrm{\Omega }}^{\mathrm{12}}\mathrm{)}\end{array}$(29)is the surface-integrated gravity (see the power series in CO).

Together with the use of the Stefan-Boltzmann law for the flux emitted at co-latitude θ$\mathit{F}\hspace{0.17em}\mathrm{\left(}\mathit{\theta }\mathrm{\right)}\mathrm{=}{\mathit{\sigma }}_{\mathrm{B}}\hspace{0.17em}{\mathit{T}}_{\mathrm{eff}}^{\mathrm{4}}\hspace{0.17em}\mathrm{\left(}\mathit{\theta }\mathrm{\right)}\mathit{,}$(30)where σ_B is the Boltzmann constant, we obtain finally the following equation for the local effective temperature ${\mathit{T}}_{\mathrm{eff}}\hspace{0.17em}\mathrm{\left(}{\mathit{V}}_{\mathrm{rot}}\mathit{,\theta }\mathrm{\right)}\mathrm{=}{\left[\frac{\mathit{L}}{{\mathit{\sigma }}_{\mathit{B}}\hspace{0.17em}\mathrm{\Sigma }\hspace{0.17em}\mathrm{\left(}{\mathit{V}}_{\mathrm{rot}}\mathrm{\right)}}\hspace{0.17em}{\mathit{g}}_{\mathrm{\perp }}\hspace{0.17em}\mathrm{\left(}{\mathit{V}}_{\mathrm{rot}}\mathit{,\theta }\mathrm{\right)}\right]}^{\mathrm{1}\mathit{/}\mathrm{4}}\mathit{.}$(31)Through the angular velocity Ω (cf. Eq. (23)), the local effective temperature in Eq. (31) and the stellar radius R(θ) in Eq. (22) depends also on the continuum Eddington factor Γ (cf. Eq. (11)) which is, in case of an homogeneous spherical star, $\mathrm{\Gamma }\mathrm{=}\frac{{\mathit{\sigma }}_{\mathrm{e}}\hspace{0.17em}\mathit{L}}{\mathrm{4}\mathit{\pi }\hspace{0.17em}\mathit{c}\hspace{0.17em}\mathit{G}\hspace{0.17em}\mathit{M}}\mathit{,}$(32)where σ_e is the electron scattering cross-section.

Then, this Eq. (32) is also used in our case to calculate a mean value of from the prescribed value of L of the non-spherical star, to be able to evaluate Eq. (22) and Eq. (31) for the stellar parameters R(θ) and T_eff   (θ) at a given co-latitude θ of interest.

However, since the effective gravity and therefore the flux vary over the surface of the rotating star, we still need to determine the local value of Γ   (θ) of the non-spherical star that can be defined as $\mathrm{\Gamma }\hspace{0.17em}\mathrm{\left(}\mathit{\theta }\mathrm{\right)}\mathrm{:}\mathrm{=}\frac{{\mathit{\sigma }}_{\mathrm{e}}\hspace{0.17em}\mathit{L}\hspace{0.17em}\mathrm{\left(}\mathit{\theta }\mathrm{\right)}}{\mathrm{4}\mathit{\pi }\hspace{0.17em}\mathit{c}\hspace{0.17em}\mathit{G}\hspace{0.17em}\mathit{M}}\mathrm{=}\frac{{\mathit{\sigma }}_{\mathrm{e}}\hspace{0.17em}{\mathit{\sigma }}_{\mathrm{B}}}{\mathit{c}\hspace{0.17em}\mathit{G}\hspace{0.17em}\mathit{M}}\hspace{0.17em}{\mathit{R}}^{\mathrm{2}}\hspace{0.17em}\mathrm{\left(}\mathit{\theta }\mathrm{\right)}\hspace{0.17em}{\mathit{T}}_{\mathrm{eff}}^{\mathrm{4}}\hspace{0.17em}\mathrm{\left(}\mathit{\theta }\mathrm{\right)}$(33)by means of the definition of the latitude-dependent luminosity $\mathit{L}\hspace{0.17em}\mathrm{\left(}\mathit{\theta }\mathrm{\right)}\mathrm{:}\mathrm{=}\mathrm{4}\hspace{0.17em}\mathit{\pi }\hspace{0.17em}{\mathit{\sigma }}_{\mathit{B}}\hspace{0.17em}{\mathit{R}}^{\mathrm{2}}\hspace{0.17em}\mathrm{\left(}\mathit{\theta }\mathrm{\right)}\hspace{0.17em}{\mathit{T}}_{\mathrm{eff}}^{\mathrm{4}}\hspace{0.17em}\mathrm{\left(}\mathit{\theta }\mathrm{\right)}\mathit{,}$(34)which is the luminosity of a corresponding spherical star with radius R(θ) and effective temperature T_eff   (θ) of the considered non-spherical star at co-latitude θ.

2.4.4. The equation of motion

Next, we wish to solve the r-component of the momentum Eq. (16), i.e. find an expression for the radial velocity component v_r of the non-spherical axisymmetric steady flow. Equation (16) can then be rewritten $\mathit{v̂}\mathit{r}\hspace{0.17em}\frac{\mathit{\partial }}{\mathit{\partial }\mathit{r̂}}\hspace{0.17em}\mathit{v̂}\mathit{r}\mathrm{-}\frac{\mathit{v̂}\begin{array}{c}\mathrm{2}\\ \mathit{\phi }\end{array}}{\mathit{r̂}}\mathrm{=}\mathrm{-}\frac{{\mathit{v̂}\mathrm{crit}}^{\mathrm{2}}}{{\mathit{r̂}}^{\mathrm{2}}}\mathrm{+}\mathit{ĝ}\begin{array}{c}\mathrm{line}\\ \mathrm{rad}\end{array}\mathrm{-}\frac{\mathrm{1}}{\mathit{\rho }}\hspace{0.17em}\frac{\mathit{\partial \rho }}{\mathit{\partial }\mathit{r̂}}$(35)in non-dimensional form. In which the following dimensionless velocities (in units of the isothermal sound speed a) $\mathit{v̂}\mathit{r}\mathrm{:}\mathrm{=}\frac{{\mathit{v}}_{\mathit{r}}}{\mathit{a}}\mathit{,} \mathit{v̂}\mathit{\phi }\mathrm{:}\mathrm{=}\frac{{\mathit{v}}_{\mathit{\phi }}}{\mathit{a}}\mathit{,} \mathit{v̂}\mathrm{crit}\hspace{0.17em}\mathrm{\left(}\mathit{\theta }\mathrm{\right)}\mathrm{:}\mathrm{=}\frac{\mathrm{1}}{\mathit{a}}\hspace{0.17em}\sqrt{\frac{\mathit{GM}}{\mathit{R}\hspace{0.17em}\mathrm{\left(}\mathit{\theta }\mathrm{\right)}}\hspace{0.17em}\mathrm{(}\mathrm{1}\mathrm{-}\mathrm{\Gamma ̅}\mathrm{)}}\mathit{,}$(36)and dimensionless line acceleration $\mathit{ĝ}\begin{array}{c}\mathrm{line}\\ \mathrm{rad}\end{array}\hspace{0.17em}\mathrm{\left(}\mathit{\theta }\mathrm{\right)}\mathrm{:}\mathrm{=}\frac{\mathit{R}\hspace{0.17em}\mathrm{\left(}\mathit{\theta }\mathrm{\right)}}{{\mathit{a}}^{\mathrm{2}}}\hspace{0.17em}{\mathit{g}}_{\mathrm{rad}}^{\mathrm{line}}\hspace{0.17em}\mathrm{\left(}\mathit{\theta }\mathrm{\right)}$(37)are used, where v_crit   (θ = 0) equals the break-up velocity of the rotating central object (usually without consideration of radiative line acceleration terms and rotational distortion). By means of Eq. (15) and applying the chain rule to the function , we obtain $\begin{array}{ccc}\frac{\mathit{\partial \rho }}{\mathit{\partial }\mathit{r̂}}& \mathrm{=}& \left(\mathrm{-}\frac{\mathrm{2}}{{\mathit{r̂}}^{\mathrm{3}}}\hspace{0.17em}\frac{\mathrm{1}}{{\mathit{v}}_{\mathit{r}}\hspace{0.17em}\mathrm{\left(}\mathit{r̂}\mathrm{\right)}}\mathrm{-}\frac{\mathrm{1}}{{\mathit{r̂}}^{\mathrm{2}}}\hspace{0.17em}\frac{\mathrm{1}}{{{\mathit{v}}_{\mathit{r}}\hspace{0.17em}\mathrm{\left(}\mathit{r̂}\mathrm{\right)}}^{\mathrm{2}}}\hspace{0.17em}\frac{\mathit{\partial }{\mathit{v}}_{\mathit{r}}\hspace{0.17em}\mathrm{\left(}\mathit{r̂}\mathrm{\right)}}{\mathit{\partial }\mathit{r̂}}\right)\hspace{0.17em}\mathit{F}\\ & \mathrm{\equiv }& \mathrm{-}\mathit{\rho }\hspace{0.17em}\mathrm{\left(}\mathit{r̂}\mathit{,\theta }\mathrm{\right)}\hspace{0.17em}\left(\frac{\mathrm{2}}{\mathit{r̂}}\mathrm{+}\frac{\mathrm{1}}{{\mathit{v}}_{\mathit{r}}\hspace{0.17em}\mathrm{\left(}\mathit{r̂}\mathit{,\theta }\mathrm{\right)}}\hspace{0.17em}\frac{\mathit{\partial }{\mathit{v}}_{\mathit{r}}\hspace{0.17em}\mathrm{\left(}\mathit{r̂}\mathit{,\theta }\mathrm{\right)}}{\mathit{\partial }\mathit{r̂}}\right)\mathrm{·}\end{array}$Using this expression for in Eq. (35) together with our relation for the azimuthal velocity, Eq. (19), we finally find the dimensionless differential equation of motion (EOM) for the radial velocity at constant co-latitude θ$\left(\mathit{v̂}\mathit{r}\mathrm{-}\frac{\mathrm{1}}{\mathit{v̂}\mathit{r}}\right)\hspace{0.17em}\frac{\mathit{\partial }}{\mathit{\partial }\hspace{0.17em}\mathit{r̂}}\hspace{0.17em}\mathit{v̂}\mathit{r}\mathrm{=}\frac{{\mathit{v̂}\mathrm{rot}}^{\mathrm{2}}\hspace{0.17em}\mathrm{\left(}\mathit{\theta }\mathrm{\right)}}{{\mathit{r̂}}^{\mathrm{3}}}\mathrm{-}\frac{{\mathit{v̂}\mathrm{crit}}^{\mathrm{2}}\hspace{0.17em}\mathrm{\left(}\mathit{\theta }\mathrm{\right)}}{{\mathit{r̂}}^{\mathrm{2}}}\mathrm{+}\frac{\mathrm{2}}{\mathit{r̂}}\mathrm{+}\mathit{ĝ}\begin{array}{c}\mathrm{line}\\ \mathrm{rad}\end{array}\hspace{0.17em}\mathrm{\left(}\mathit{\theta }\mathrm{\right)}\mathit{,}$(38)that is now independent of ρ.

2.5. The line acceleration term and the final equation of motion

To derive a sophisticated mathematical expression for the radiative line acceleration $\mathit{ĝ}\begin{array}{c}\mathrm{line}\\ \mathrm{rad}\end{array}\hspace{0.17em}\mathrm{\left(}\mathit{r̂}\mathit{,\theta }\mathrm{\right)}$ at a constant co-latitude θ of interest as a function of radius r only, we demand the same physically motivated mathematical properties as described in Paper I (as for the case of a radiation-driven spherical wind): $\mathit{ĝ}\begin{array}{c}\mathrm{line}\\ \mathrm{rad}\end{array}\mathrm{\left(}\mathit{r̂}\mathit{,\theta }\mathrm{\right)}\mathrm{=}\frac{\mathit{ĝ}\mathrm{0}}{\mathit{r̂}\mathrm{1}\mathrm{+}\mathit{\delta }}\hspace{0.17em}{\left(\mathrm{1}\mathrm{-}\frac{\mathit{r̂}\mathrm{0}}{{\mathit{r̂}}^{\mathit{\delta }}}\right)}^{\mathit{\gamma }}\mathrm{\equiv }\frac{\mathit{ĝ}\mathrm{0}}{\mathit{r̂}\mathrm{1}\mathrm{+}\mathit{\delta }\hspace{0.17em}\left(\mathrm{1}\mathrm{+}\mathit{\gamma }\right)}\hspace{0.17em}{\mathrm{(}{\mathit{r̂}}^{\mathit{\delta }}\mathrm{-}{\mathit{r̂}\mathrm{0}}^{\mathrm{)}}}^{\mathit{\gamma }}\mathit{.}$(39)This function is independent of $\mathit{v̂}{\mathit{r}}_{}$ and () and dependent on only, at constant co-latitude θ. Note that, herein, the set of four parameters all depend on latitude: ĝ₀ = ĝ₀(θ), , γ = γ(θ), δ = δ(θ), due to the rotation of the central star.

The line force is zero at radius near the stellar photosphere () and everywhere else positive for . To guarantee the decrease of $\mathit{ĝ}\begin{array}{c}\mathrm{line}\\ \mathrm{rad}\end{array}\hspace{0.17em}\mathrm{\left(}\mathit{r̂}\mathrm{\right)}$ as with increasing radial distance from the central star at intermediate radii (at the right side of the $\mathit{ĝ}\begin{array}{c}\mathrm{line}\\ \mathrm{rad}\end{array}$ peak) in particular, we had to introduce the parameter δ in addition to γ (where 0 < γ ≲ 1 and 0 < δ ≲ 1).

Hence, the equation of motion (38) for each latitude becomes $\begin{array}{ccc}\left(\mathit{v̂}\mathit{r}\mathrm{-}\frac{\mathrm{1}}{\mathit{v̂}\mathit{r}}\right)\hspace{0.17em}\frac{\mathit{\partial }}{\mathit{\partial }\hspace{0.17em}\mathit{r̂}}\hspace{0.17em}\mathit{v̂}\mathit{r}\mathrm{=}\frac{{\mathit{v̂}\mathrm{rot}}^{\mathrm{2}}}{{\mathit{r̂}}^{\mathrm{3}}}\mathrm{-}\frac{{\mathit{v̂}\mathrm{crit}}^{\mathrm{2}}}{{\mathit{r̂}}^{\mathrm{2}}}\mathrm{+}\frac{\mathrm{2}}{\mathit{r̂}}\mathrm{+}\frac{\mathit{ĝ}\mathrm{0}}{\mathit{r̂}\mathrm{1}\hspace{0.17em}\mathrm{+}\hspace{0.17em}\mathit{\delta }\hspace{0.17em}\left(\mathrm{1}\hspace{0.17em}\mathrm{+}\hspace{0.17em}\mathit{\gamma }\right)}\hspace{0.17em}{\mathrm{(}{\mathit{r̂}}^{\mathit{\delta }}\mathrm{-}{\mathit{r̂}\mathrm{0}}^{\mathrm{)}}}^{\mathit{\gamma }}\mathit{.}& & \end{array}$(40)v_rot is the azimuthal velocity of the surface of the inner rotating star (object) at co-latitude θ of interest. This equation is fully solvable analytically as in the spherical case without rotation (cf. Paper I).

2.6. Analytical solutions of the equation of motion

2.6.1. The critical point and critical solutions

The EOM (Eqs. (38) and (40)) yields several families of solutions that have quite different mathematical behaviour and physical significance (cf. Fig. 2).

The left hand side of Eq. (40) vanishes for at the critical radius , where . That is, the critical point velocity here is equal to the isothermal sound speed , and the critical radius is just the sonic radius $\mathit{r̂}\mathrm{c}\mathrm{\equiv }\mathit{r̂}\mathrm{s}\mathit{,}$(41)as is also the case for thermal winds or mass accretion events (where $\mathit{ĝ}\begin{array}{c}\mathrm{line}\\ \mathrm{rad}\end{array}\mathrm{\approx }\mathrm{0}$).

We are now interested in under which conditions one can obtain a continuous and smooth trans-sonic flow through the critical point of Eq. (40). For the case of a stellar wind, this means how to obtain a smooth transition from subsonic and subcritical flow () at small to supercritical and supersonic flow () at large , when this critical solution has a finite positive slope ( at (cf. solid curve 1 in Fig. 2)? ⁵ Then, it is evident from the left hand side of Eq. (40) that one can obtain such a trans-sonic wind, if the right hand side (1) vanishes at the critical radius , (2) is negative for , and (3) is positive for .

The opposite situation occurs for the case of mass accretion in e.g. a collapsing cloud. If (, we obtain the second unique trans-sonic and critical solution in which is monotonically decreasing from supersonic speeds for , e.g. near the protostar, to subsonic speeds for at the outer edge of the cloud (see also the second solid line 2, in Fig. 2, for the case of a corresponding accretion flow with a star as the central object).

Here we are mainly interested in the critical wind solution of Eq. (40). The right hand side of Eq. (40) vanishes at the critical radius that solves the equation $\mathrm{2}\hspace{0.17em}{\mathit{r̂}}^{\mathit{\delta }\hspace{0.17em}\left(\mathrm{1}\hspace{0.17em}\mathrm{+}\hspace{0.17em}\mathit{\gamma }\right)\hspace{0.17em}\mathrm{+}\hspace{0.17em}\mathrm{1}}\mathrm{-}{\mathit{v̂}}_{\mathrm{crit}}^{\mathrm{2}}\hspace{0.17em}{\mathit{r̂}}^{\mathit{\delta }\hspace{0.17em}\left(\mathrm{1}\hspace{0.17em}\mathrm{+}\hspace{0.17em}\mathit{\gamma }\right)}\mathrm{+}\mathit{ĝ}\mathrm{0}\hspace{0.17em}\mathit{r̂}\hspace{0.17em}{\mathrm{(}{\mathit{r̂}}^{\mathit{\delta }}\mathrm{-}{\mathit{r̂}\mathrm{0}}^{\mathrm{)}}}^{\mathit{\gamma }}\mathrm{+}{\mathit{v̂}\mathrm{rot}}^{\mathrm{2}}\hspace{0.17em}{\mathit{r̂}}^{\mathit{\delta }\hspace{0.17em}\left(\mathrm{1}\hspace{0.17em}\mathrm{+}\hspace{0.17em}\mathit{\gamma }\right)\hspace{0.17em}\mathrm{-}\hspace{0.17em}\mathrm{1}}\mathrm{=}\mathrm{0.}$(42)Therefore, the critical radius (for each latitude) has to be determined numerically by means of the above equation and the line term parameters $\mathit{ĝ}{\mathrm{0}}_{}$, γ, δ and $\mathit{r̂}{\mathrm{0}}_{}$, depending on the rotational speed V_rot of the central object. However, if one assumes values of γ and δ close to 1, one can provide a good analytical approximation (for the solution of Eq. (42)) for the critical radius ${\mathit{r̂}}_{\mathrm{c}}\mathrm{\approx }\frac{\mathrm{1}}{\mathrm{4}}\hspace{0.17em}\left(\mathrm{(}{\mathit{v̂}}_{\mathrm{crit}}^{\mathrm{2}}\mathrm{-}{\mathit{ĝ}\mathrm{0}}^{\mathrm{)}}\mathrm{+}\sqrt{\mathrm{(}{\mathit{v̂}}_{\mathrm{crit}}^{\mathrm{2}}\mathrm{-}{\mathit{ĝ}\mathrm{0}}^{\mathrm{)}}{}^{\mathrm{2}}\mathrm{+}\mathrm{8}\hspace{0.17em}\left(\mathit{ĝ}\mathrm{0}\hspace{0.17em}\mathit{r̂}\mathrm{0}\mathrm{-}{\mathit{v̂}\mathrm{rot}}^{\mathrm{2}}\right)}\right)\mathit{,}$(43)if the rotation speed at co-latitude θ fulfills the condition $\begin{array}{ccc}\mathit{v̂}\mathrm{rot}\mathit{<}\sqrt{\frac{\mathrm{1}}{\mathrm{8}}\hspace{0.17em}\mathrm{(}{\mathit{v̂}}_{\mathrm{crit}}^{\mathrm{2}}\mathrm{-}{\mathit{ĝ}\mathrm{0}}^{\mathrm{)}}{}^{\mathrm{2}}\mathrm{+}\mathit{ĝ}\mathrm{0}\hspace{0.17em}\mathit{r̂}\mathrm{0}}\mathit{.}& & \end{array}$For the simpler case of a thermal wind (or any other non line-driven mass flow), where $\mathit{ĝ}\begin{array}{c}\mathrm{line}\\ \mathrm{rad}\end{array}$ can be set to zero (in Eq. (42)), we obtain the analytical solution ${\mathit{r̂}}_{\mathrm{c}}\mathrm{=}\frac{\mathrm{1}}{\mathrm{4}}\hspace{0.17em}\left({\mathit{v̂}}_{\mathrm{crit}}^{\mathrm{2}}\mathrm{+}\sqrt{{\mathit{v̂}}_{\mathrm{crit}}^{\mathrm{4}}\mathrm{-}\mathrm{8}\hspace{0.17em}{\mathit{v̂}\mathrm{rot}}^{\mathrm{2}}}\right)\mathit{.}$(44)Equation (44), and Eq. (43) under the assumption of a constant remaining line force (i.e. $\mathit{ĝ}{\mathrm{0}}_{}$ value), imply that the critical radius moves closer to the inner core radius (at any latitude besides of that of the pole) with increasing rotational speed V_rot for these particular cases of a trans-sonic flow.

2.6.2. Solving the equation of motion

The equation of motion (40) can be solved by first integrating the left hand side over $\mathit{v̂}{\mathit{r}}_{}$, and then integrating the right hand side over , separately, which yields $\begin{array}{ccc}{\mathit{v̂}\mathit{r}}^{\mathrm{2}}\mathrm{-}\ln {\mathit{v̂}\mathit{r}}^{\mathrm{2}}& \mathrm{=}& \mathrm{2}\frac{{\mathit{v̂}}_{\mathrm{crit}}^{\mathrm{2}}}{\mathit{r̂}}\mathrm{+}\mathrm{4}\hspace{0.17em}\ln \mathit{r̂}\\ & & \mathrm{+}\frac{\mathrm{2}}{\mathit{r̂}\mathrm{0}}\hspace{0.17em}\frac{\mathit{ĝ}\mathrm{0}}{\mathit{\delta }\hspace{0.17em}\left(\mathrm{1}\mathrm{+}\mathit{\gamma }\right)}{\left(\mathrm{1}\mathrm{-}\frac{\mathit{r̂}\mathrm{0}}{{\mathit{r̂}}^{\mathit{\delta }}}\right)}^{\mathrm{1}\mathrm{+}\mathit{\gamma }}\mathrm{-}\frac{{\mathit{v̂}}_{\mathrm{rot}}^{\mathrm{2}}}{{\mathit{r̂}}^{\mathrm{2}}}\mathrm{+}\mathit{C,}\end{array}$(45)with the right hand side of Eq. (45) denoted as the function $\begin{array}{ccc}\mathit{f}\hspace{0.17em}\mathrm{\left(}\mathit{r̂}\mathit{,\theta }\mathrm{;}{\mathit{r̂}}^{\mathrm{\prime }}\mathit{,}{\mathit{v̂}\mathit{r}}^{\mathrm{\prime }}\mathrm{\right)}& \mathrm{:}\mathrm{=}& \mathrm{2}\frac{{\mathit{v̂}}_{\mathrm{crit}}^{\mathrm{2}}}{\mathit{r̂}}\mathrm{+}\mathrm{4}\ln \mathit{r̂}\\ & \mathrm{+}& \frac{\mathrm{2}}{\mathit{r̂}\mathrm{0}}\frac{\mathit{ĝ}\mathrm{0}}{\mathit{\delta }\left(\mathrm{1}\mathrm{+}\mathit{\gamma }\right)}{\left(\mathrm{1}\mathrm{-}\frac{\mathit{r̂}\mathrm{0}}{{\mathit{r̂}}^{\mathit{\delta }}}\right)}^{\mathrm{1}\mathrm{+}\mathit{\gamma }}\mathrm{-}\frac{{\mathit{v̂}}_{\mathrm{rot}}^{\mathrm{2}}}{{\mathit{r̂}}^{\mathrm{2}}}\mathrm{+}\mathit{C}\hspace{0.17em}\mathrm{\left(}{\mathit{r̂}}^{\mathrm{\prime }}\mathit{,\theta }\mathit{,}{\mathit{v̂}\mathit{r}}^{\mathrm{\prime }}\mathrm{\right)}\end{array}$(46)with the constant of integration C, that is determined by the boundary condition of the radial velocity ${\mathit{v̂}{\mathit{r}}_{}}^{\mathrm{\prime }}$ at a given location .

From Eq. (45), we can determine for the solution that passes through the particular point , $\begin{array}{ccc}\mathit{C}\hspace{0.17em}\mathrm{\left(}{\mathit{r̂}}^{\mathrm{\prime }}\mathit{,\theta }\mathit{,}{\mathit{v̂}\mathit{r}}^{\mathrm{\prime }}\mathrm{\right)}& \mathrm{=}& {\mathit{v̂}\mathit{r}}^{\mathrm{\prime }\mathrm{2}}\mathrm{-}\ln {\mathit{v̂}\mathit{r}}^{\mathrm{\prime }\mathrm{2}}\mathrm{-}\mathrm{2}\frac{{\mathit{v̂}}_{\mathrm{crit}}^{\mathrm{2}}}{\mathit{r̂}\mathrm{\prime }}\mathrm{-}\mathrm{4}\hspace{0.17em}\ln {\mathit{r̂}}^{\mathrm{\prime }}\\ & & \mathrm{-}\hspace{0.17em}\frac{\mathrm{2}}{\mathit{r̂}\mathrm{0}}\hspace{0.17em}\frac{\mathit{ĝ}\mathrm{0}}{\mathit{\delta }\left(\mathrm{1}\mathrm{+}\mathit{\gamma }\right)}{\left(\mathrm{1}\mathrm{-}\frac{\mathit{r̂}\mathrm{0}}{{\mathit{r̂}}^{\mathrm{\prime }\mathit{\delta }}}\right)}^{\mathrm{1}\mathrm{+}\mathit{\gamma }}\mathrm{+}\frac{{\mathit{v̂}}_{\mathrm{rot}}^{\mathrm{2}}}{{\mathit{r̂}}^{\mathrm{\prime }\mathrm{2}}}\mathrm{·}\end{array}$(47)Therefore the function f in Eq. (46) becomes $\begin{array}{ccc}\mathit{f}\hspace{0.17em}\mathrm{\left(}\mathit{r̂}\mathit{,\theta }\mathrm{;}{\mathit{r̂}}^{\mathrm{\prime }}\mathit{,}{\mathit{v̂}\mathit{r}}^{\mathrm{\prime }}\mathrm{\right)}& \mathrm{=}& {\mathit{v̂}\mathit{r}}^{\mathrm{\prime }\mathrm{2}}\mathrm{-}\ln {\mathit{v̂}\mathit{r}}^{\mathrm{\prime }\mathrm{2}}\mathrm{+}\mathrm{2}\hspace{0.17em}{\mathit{v̂}}_{\mathrm{crit}}^{\mathrm{2}}\left(\frac{\mathrm{1}}{\mathit{r̂}}\mathrm{-}\frac{\mathrm{1}}{\mathit{r̂}\mathrm{\prime }}\right)\mathrm{+}\mathrm{4}\hspace{0.17em}\ln \left(\frac{\mathit{r̂}}{\mathit{r̂}\mathrm{\prime }}\right)\\ & & \mathrm{+}\hspace{0.17em}\frac{\mathrm{2}}{\mathit{r̂}\mathrm{0}}\hspace{0.17em}\frac{\mathit{ĝ}\mathrm{0}}{\mathit{\delta }\left(\mathrm{1}\mathrm{+}\mathit{\gamma }\right)}\left[{\left(\mathrm{1}\mathrm{-}\frac{\mathit{r̂}\mathrm{0}}{{\mathit{r̂}}^{\mathit{\delta }}}\right)}^{\mathrm{1}\mathrm{+}\mathit{\gamma }}\mathrm{-}{\left(\mathrm{1}\mathrm{-}\frac{\mathit{r̂}\mathrm{0}}{{\mathit{r̂}}^{\mathrm{\prime }\mathit{\delta }}}\right)}^{\mathrm{1}\mathrm{+}\mathit{\gamma }}\right]\\ & & \mathrm{+}\hspace{0.17em}{\mathit{v̂}}_{\mathrm{rot}}^{\mathrm{2}}\hspace{0.17em}\left(\frac{\mathrm{1}}{{\mathit{r̂}}^{\mathrm{\prime }\mathrm{2}}}\mathrm{-}\frac{\mathrm{1}}{{\mathit{r̂}}^{\mathrm{2}}}\right)\mathrm{·}\end{array}$(48)And from this, Eq. (45) now reads $\begin{array}{ccc}{\mathit{v̂}\mathit{r}}^{\mathrm{2}}\mathrm{-}\ln {\mathit{v̂}\mathit{r}}^{\mathrm{2}}\mathrm{=}\mathit{f}\hspace{0.17em}\mathrm{\left(}\mathit{r̂}\mathit{,\theta }\mathrm{;}{\mathit{r̂}}^{\mathrm{\prime }}\mathit{,}{\mathit{v̂}\mathit{r}}^{\mathrm{\prime }}\mathrm{\right)}& & \end{array}$or, equivalently, $\mathrm{-}{\mathit{v̂}\mathit{r}}^{\mathrm{2}}\hspace{0.17em}{\mathrm{e}}^{\mathrm{-}{\mathit{v̂}\mathit{r}}^{\mathrm{2}}}\mathrm{=}\mathrm{-}{\mathrm{e}}^{\mathrm{-}\mathit{f}\hspace{0.17em}\mathrm{\left(}\mathit{r̂}\mathit{,\theta }\mathrm{;}{\mathit{r̂}}^{\mathrm{\prime }}\mathit{,}{\mathit{v̂}\mathit{r}}^{\mathrm{\prime }}\mathrm{\right)}}\mathit{,}$(49)which is solved explicitly and fully analytically in terms of the Lambert W function (cf. Corless et al. 1993, 1996).

2.6.3. The solution(s) of the equation of motion

It is now possible to provide an explicit analytical expression for the solution $\mathit{v̂}{\mathit{r}}_{}$ of the equation of motion (40), or Eq. (49), by means of the W function (cf. Paper I). If we compare Eq. (49) with the defining equation of the Lambert W function ${\mathit{W}}_{\mathit{k}}\mathrm{\left(}\mathit{x}\mathrm{\right)}\hspace{0.17em}{\mathrm{e}}^{{\mathit{W}}_{\mathit{k}}\mathrm{\left(}\mathit{x}\mathrm{\right)}}\mathrm{=}\mathit{x,}$(50)we find that $\begin{array}{ccc}\mathrm{-}{\mathit{v̂}\mathit{r}}^{\mathrm{2}}\mathrm{=}{\mathit{W}}_{\mathit{k}}\mathrm{\left(}\mathit{x}\mathrm{\right)}& & \end{array}$or $\mathit{v̂}\mathit{r}\mathrm{=}\mathrm{\pm }\sqrt{\mathrm{-}{\mathit{W}}_{\mathit{k}}\mathrm{\left(}\mathit{x}\mathrm{\right)}}$(51)is the general solution of the equation of motion that passes through the point , with the argument function of the W function $\mathit{x}\hspace{0.17em}\mathrm{\left(}\mathit{r̂}\mathit{,\theta }\mathrm{\right)}\mathrm{=}\mathrm{-}{\mathrm{e}}^{\mathrm{-}\mathit{f}\hspace{0.17em}\mathrm{\left(}\mathit{r̂}\mathit{,\theta }\mathrm{;}{\mathit{r̂}}^{\mathrm{\prime }}\mathit{,}{\mathit{v̂}\mathit{r}}^{\mathrm{\prime }}\mathrm{\right)}}\mathit{.}$(52)Since the argument of the W function in Eq. (51) is always real and negative, it is guaranteed that the argument of the square root never becomes negative, and hence the solution is always real.

Inserting from Eq. (48) into Eq. (52) yields $\begin{array}{ccc}\mathit{x}\hspace{0.17em}\mathrm{\left(}\mathit{r̂}\mathit{,\theta }\mathrm{;}{\mathit{r̂}}^{\mathrm{\prime }}\mathit{,}{\mathit{v̂}\mathit{r}}^{\mathrm{\prime }}\mathrm{\right)}& \mathrm{=}& \mathrm{-}{\left(\frac{\mathit{r̂}\mathrm{\prime }}{\mathit{r̂}}\right)}^{\mathrm{4}}{\mathit{v̂}\mathit{r}}^{\mathrm{\prime }\mathrm{2}}\exp [\mathrm{-}{\mathit{v̂}}_{\mathrm{rot}}^{\mathrm{2}}\hspace{0.17em}\left(\frac{\mathrm{1}}{{\mathit{r̂}}^{\mathrm{\prime }\mathrm{2}}}\mathrm{-}\frac{\mathrm{1}}{{\mathit{r̂}}^{\mathrm{2}}}\right)\mathrm{-}\mathrm{2}\hspace{0.17em}{\mathit{v̂}}_{\mathrm{crit}}^{\mathrm{2}}\hspace{0.17em}\left(\frac{\mathrm{1}}{\mathit{r̂}}\mathrm{-}\frac{\mathrm{1}}{\mathit{r̂}\mathrm{\prime }}\right)\\ & & \mathrm{-}\frac{\mathrm{2}}{\mathit{r̂}\mathrm{0}}\frac{\mathit{ĝ}\mathrm{0}}{\mathit{\delta }\left(\mathrm{1}\mathrm{+}\mathit{\gamma }\right)}\left({\left(\mathrm{1}\mathrm{-}\frac{\mathit{r̂}\mathrm{0}}{{\mathit{r̂}}^{\mathit{\delta }}}\right)}^{\mathrm{1}\mathrm{+}\mathit{\gamma }}\mathrm{-}{\left(\mathrm{1}\mathrm{-}\frac{\mathit{r̂}\mathrm{0}}{{\mathit{r̂}}^{\mathrm{\prime }\mathit{\delta }}}\right)}^{\mathrm{1}\mathrm{+}\mathit{\gamma }}\right)\mathrm{-}{\mathit{v̂}\mathit{r}}^{\mathrm{\prime }\mathrm{2}}]\end{array}$(53)the general expression for the argument function x depending on the parameters .

Thus, for the trans-sonic case of a stellar wind or general accretion flow, where and , the analytical solution is $\mathit{v̂}\mathit{r}\hspace{0.17em}\mathrm{\left(}\mathit{r̂}\mathit{,\theta }\mathrm{\right)}\mathrm{=}\mathrm{\pm }\sqrt{\mathrm{-}{\mathit{W}}_{\mathit{k}}\mathrm{\left(}\mathit{x}\hspace{0.17em}\mathrm{\right(}\mathit{r̂}\mathit{,\theta }\mathrm{\left)}\mathrm{\right)}}$(54)with $\begin{array}{ccc}\mathit{x}\hspace{0.17em}\mathrm{\left(}\mathit{r̂}\mathit{,\theta }\mathrm{\right)}& \mathrm{=}& \mathrm{-}{\left(\frac{\mathit{r̂}\mathrm{c}}{\mathit{r̂}}\right)}^{\mathrm{4}}\hspace{0.17em}\exp [\mathrm{-}{\mathit{v̂}}_{\mathrm{rot}}^{\mathrm{2}}\hspace{0.17em}\left(\frac{\mathrm{1}}{{\mathit{r̂}\mathrm{c}}^{\mathrm{2}}}\mathrm{-}\frac{\mathrm{1}}{{\mathit{r̂}}^{\mathrm{2}}}\right)\mathrm{-}\mathrm{2}\hspace{0.17em}{\mathit{v̂}}_{\mathrm{crit}}^{\mathrm{2}}\hspace{0.17em}\left(\frac{\mathrm{1}}{\mathit{r̂}}\mathrm{-}\frac{\mathrm{1}}{\mathit{r̂}\mathrm{c}}\right)\\ & & \mathrm{-}\frac{\mathrm{2}}{\mathit{r̂}\mathrm{0}}\frac{\mathit{ĝ}\mathrm{0}}{\mathit{\delta }\left(\mathrm{1}\mathrm{+}\mathit{\gamma }\right)}\hspace{0.17em}\left({\left(\mathrm{1}\mathrm{-}\frac{\mathit{r̂}\mathrm{0}}{{\mathit{r̂}}^{\mathit{\delta }}}\right)}^{\mathrm{1}\hspace{0.17em}\mathrm{+}\hspace{0.17em}\mathit{\gamma }}\mathrm{-}{\left(\mathrm{1}\mathrm{-}\frac{\mathit{r̂}\mathrm{0}}{{\mathit{r̂}\mathrm{c}}^{\mathit{\delta }}}\right)}^{\mathrm{1}\hspace{0.17em}\mathrm{+}\hspace{0.17em}\mathit{\gamma }}\right)\mathrm{-}\mathrm{1}]\mathit{.}\end{array}$(55)We are only interested in the possible two real values of W(x), the k = 0, − 1–branches in Eq. (51) or Eq. (54), where x is real and −1/e ≤ x < 0. The branch point at x =  −1/e, where these two branches meet, corresponds to the critical point , where the velocity in Eq. (54) becomes ≡ 1. Depending on which branch of W one is approaching this point x =  −1/e, one obtains a different shape of the –curve, i.e. a stellar wind or a collapsing system.

However, to determine which of the two branches to choose at a certain range of radius between and as to guarantee a continuous, monotonically increasing, and smooth trans-sonic flow (as, e.g., in the case of a stellar wind), one needs to examine the behaviour of the argument function of W, in Eq. (55), with radius at a given co-latitude θ. Then, the argument function decreases monotonically from the stellar radius (with value of nearly zero) to its minimum at with x =  −1/e to afterwards increase monotonically again.

Finally, a detailed investigation (cf. Paper I) yields the following amount of the radial velocity component (i.e. the axisymmetric two-dimensional trans-sonic analytical solution of our equation of motion for a rotating and expanding or collapsing system) for a given latitude, i.e. polar angle θ

• (a)

  for the case of a stellar wind, and Eq. (54) (and thepositive sign in front of the root) with the argument function inEq. (55), choosing the branch$\mathit{k}\mathrm{=}\{\begin{array}{c}\\ \hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\mathrm{0}& \mathrm{for}& \mathrm{1}\mathrm{\le }\mathit{r̂}\mathrm{\le }\mathit{r̂}\mathrm{c}\hspace{0.17em}\mathrm{\left(}\mathit{\theta }\mathrm{\right)}\\ -1& \mathrm{for}& \mathit{r̂}\mathit{>}\mathit{r̂}\mathrm{c}\hspace{0.17em}\mathrm{\left(}\mathit{\theta }\mathrm{\right)}\end{array}$(56)of the W function at a certain radius , and

• (b)

  in case of a general accretion flux, as well, by Eq. (54) (but now with the negative sign in front of the root) and the argument function in Eq. (55), choosing the branch $\mathit{k}\mathrm{=}\{\begin{array}{c}\\ -1& \mathrm{for}& \mathrm{1}\mathrm{\le }\mathit{r̂}\mathrm{\le }\mathit{r̂}\mathrm{c}\hspace{0.17em}\mathrm{\left(}\mathit{\theta }\mathrm{\right)}\\ \hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\mathrm{0}& \mathrm{for}& \mathit{r̂}\mathit{>}\mathit{r̂}\mathrm{c}\hspace{0.17em}\mathrm{\left(}\mathit{\theta }\mathrm{\right)}\end{array}$(57)depending on the location , where

• (c)

  in the special cases of a thermal wind or collapsing system like a collapsing protostellar cloud, the argument function simplifies to $\mathit{x}\hspace{0.17em}\mathrm{\left(}\mathit{r̂}\mathit{,\theta }\mathrm{\right)}\mathrm{=}\mathrm{-}{\left(\frac{\mathit{r̂}\mathrm{c}}{\mathit{r̂}}\right)}^{\mathrm{4}}\exp \left[\mathrm{-}{\mathit{v̂}}_{\mathrm{rot}}^{\mathrm{2}}\left(\frac{\mathrm{1}}{{\mathit{r̂}\mathrm{c}}^{\mathrm{2}}}\mathrm{-}\frac{\mathrm{1}}{{\mathit{r̂}}^{\mathrm{2}}}\right)\mathrm{-}\mathrm{2}\hspace{0.17em}{\mathit{v̂}}_{\mathrm{crit}}^{\mathrm{2}}\left(\frac{\mathrm{1}}{\mathit{r̂}}\mathrm{-}\frac{\mathrm{1}}{\mathit{r̂}\mathrm{c}}\right)\mathrm{-}\mathrm{1}\right]$(58)with $\mathit{r̂}{\mathrm{c}}_{}$ given by Eq. (44), while choosing the same branches of W as mentioned above in case (a), or (b) respectively, and the appropriate sign in Eq. (54) (in front of the root).

Accordingly, one obtains different expressions for the density distribution, by Eq. (15) that depends on , which is therefore also piece-wise defined for these different ranges of radius .

In addition to these two critical solutions (type 1 and 2, cf. numbering in Fig. 2), already discussed, that pass through the critical point (i.e. sonic point), all the other four types of solutions were obtained from our general velocity law, Eq. (51) with Eq. (53), choosing the following branches of the W function and values of (), for the point we demand the solution to go through:

• Type-3: Subsonic (subcritical) solutions$\begin{array}{ccc}\mathit{k}\mathrm{=}\mathrm{0}\mathit{,} {\mathit{r̂}}^{\mathrm{\prime }}\mathrm{=}{\mathit{r̂}}_{\mathrm{c}}\hspace{0.17em}\mathrm{\left(}\mathit{\theta }\mathrm{\right)}\mathit{,} {\mathit{v̂}\mathit{r}}^{\mathrm{\prime }}\mathit{<}\mathrm{1}& & \end{array}$

• Type-4: Supersonic (supercritical) solutions $\begin{array}{ccc}\mathit{k}\mathrm{=}\mathrm{-}\mathrm{1}\mathit{,} {\mathit{r̂}}^{\mathrm{\prime }}\mathrm{=}{\mathit{r̂}}_{\mathrm{c}}\hspace{0.17em}\mathrm{\left(}\mathit{\theta }\mathrm{\right)}\mathit{,} {\mathit{v̂}\mathit{r}}^{\mathrm{\prime }}\mathit{>}\mathrm{1}& & \end{array}$

• Type-5 and 6: Double-valued solutions $\begin{array}{ccc}\mathit{k}\mathrm{=}\mathrm{0} \mathrm{and} \mathit{k}\mathrm{=}\mathrm{-}\mathrm{1} \mathrm{for}& & \end{array}$$\begin{array}{ccc}{\mathit{r̂}}^{\mathrm{\prime }}\mathit{,}{\mathit{v̂}\mathit{r}}^{\mathrm{\prime }}\mathrm{=}\mathrm{arbitrarily,}\hspace{0.17em}\mathrm{where} {\mathit{r̂}}^{\mathrm{\prime }}\ne {\mathit{r̂}}_{\mathrm{c}}\hspace{0.17em}\mathrm{\left(}\mathit{\theta }\mathrm{\right)}\mathit{,} {\mathit{v̂}\mathit{r}}^{\mathrm{\prime }}\ne \mathrm{1.}& & \end{array}$

Type-3 solutions are everywhere subsonic (choosing only the principal branch, k = 0, of the W function). Those of type 4 are everywhere supersonic (selecting only the k =  −1-branch in the velocity law), and those of type 5 and 6 are double-valued, composed of both the k = 0 and k =  −1-branch, below and above the sonic line, respectively. In this connection, the two sub- and supersonic pairs of curves of this last mentioned types, subdivided into (5a, 6a) and (5b, 6b) in Fig. 2, belong together. They are fixed, not only by the same chosen branch of W in Eq. (51), but also by the same selected parameters for the solution through the identical given point ().

Subsequently, we derive an analytical expression for the wind solution in the supersonic approximation (that is only valid in the supersonic region and is not supposed to be applied to the subsonic region, where it even becomes imaginary, particularly in our wind model in the range of ). The reasons for the necessity of deriving this approximated solution are given in our previous paper (Paper I) for the solution of a spherical wind.

2.6.4. Approximated solution of the equation of motion

By neglecting the pressure term in the equation of motion (35), which is a good approximation for the stellar wind solution in the supersonic region with , Eq. (40) becomes $\mathit{v̂}\mathit{r}\hspace{0.17em}\frac{\mathit{\partial }}{\mathit{\partial }\hspace{0.17em}\mathit{r̂}}\hspace{0.17em}\mathit{v̂}\mathit{r}\mathrm{=}\frac{{\mathit{v̂}\mathrm{rot}}^{\mathrm{2}}}{{\mathit{r̂}}^{\mathrm{3}}}\mathrm{-}\frac{{\mathit{v̂}\mathrm{crit}}^{\mathrm{2}}}{{\mathit{r̂}}^{\mathrm{2}}}\mathrm{+}\frac{\mathit{ĝ}\mathrm{0}}{\mathit{r̂}\mathrm{1}\hspace{0.17em}\mathrm{+}\hspace{0.17em}\mathit{\delta }\hspace{0.17em}\left(\mathrm{1}\hspace{0.17em}\mathrm{+}\hspace{0.17em}\mathit{\gamma }\right)}\hspace{0.17em}\mathrm{(}{\mathit{r̂}}^{\mathit{\delta }}\mathrm{-}{\mathit{r̂}\mathrm{0}}^{\mathrm{)}}\gamma$(59)at the given latitude of interest. This simplified equation of motion can again be solved by first integrating the left hand side over $\mathit{v̂}{\mathit{r}}_{}$, and then integrating the right hand side over , separately, which yields ${\mathit{v̂}\mathit{r}}^{\mathrm{2}}\mathrm{=}\mathrm{-}\frac{{\mathit{v̂}\mathrm{rot}}^{\mathrm{2}}}{{\mathit{r̂}}^{\mathrm{2}}}\mathrm{+}\mathrm{2}\frac{{\mathit{v̂}}_{\mathrm{crit}}^{\mathrm{2}}}{\mathit{r̂}}\mathrm{+}\frac{\mathrm{2}}{\mathit{r̂}\mathrm{0}}\hspace{0.17em}\frac{\mathit{ĝ}\mathrm{0}}{\mathit{\delta }\hspace{0.17em}\left(\mathrm{1}\mathrm{+}\mathit{\gamma }\right)}{\left(\mathrm{1}\mathrm{-}\frac{\mathit{r̂}\mathrm{0}}{{\mathit{r̂}}^{\mathit{\delta }}}\right)}^{\mathrm{1}\hspace{0.17em}\mathrm{+}\hspace{0.17em}\mathit{\gamma }}\mathrm{+}\mathit{C}\mathit{.}$(60)To determine the integration constant C, we assume a boundary condition for the wind velocity at radius and polar angle θ, close to the stellar photosphere, i.e. $\mathit{C}\hspace{0.17em}\mathrm{(}{\mathit{r̂}}^{\mathrm{\prime }}\mathrm{=}\mathit{r̂}\begin{array}{c}\mathrm{1}\mathit{/}\mathit{\delta }\\ \mathrm{0}\end{array}\mathit{,}{\mathit{v̂}\mathit{r}}^{\mathrm{\prime }}\mathrm{=}\mathrm{0}\mathrm{)}\mathrm{=}\mathrm{-}\mathrm{2}\hspace{0.17em}\frac{{\mathit{v̂}}_{\mathrm{crit}}^{\mathrm{2}}}{{\mathit{r̂}\mathrm{0}}^{\mathrm{1}\mathit{/}\mathit{\delta }}}\mathrm{+}\frac{{\mathit{v̂}\mathrm{rot}}^{\mathrm{2}}}{{\mathit{r̂}\mathrm{0}}^{\mathrm{2}\mathit{/}\mathit{\delta }}}\mathrm{·}$(61)Thus, from Eq. (60), the approximated wind solution reads $\begin{array}{ccc}\mathit{v̂}\mathit{r}\hspace{0.17em}\mathrm{\left(}\mathit{r̂}\mathit{,\theta }\mathrm{\right)}& \mathrm{=}& [\frac{\mathrm{2}}{\mathit{r̂}\mathrm{0}}\hspace{0.17em}\left({\mathit{v̂}}_{\mathrm{crit}}^{\mathrm{2}}\left(\frac{\mathit{r̂}\mathrm{0}}{\mathit{r̂}}\mathrm{-}{\mathit{r̂}\mathrm{0}}^{\mathrm{1}\mathrm{-}\mathrm{1}\mathit{/}\mathit{\delta }}\right)\mathrm{+}\frac{\mathit{ĝ}\mathrm{0}}{\mathit{\delta }\left(\mathrm{1}\mathrm{+}\mathit{\gamma }\right)}{\left(\mathrm{1}\mathrm{-}\frac{\mathit{r̂}\mathrm{0}}{{\mathit{r̂}}^{\mathit{\delta }}}\right)}^{\mathrm{1}\mathrm{+}\mathit{\gamma }}\right)\\ & & {\mathrm{+}{\mathit{v̂}\mathrm{rot}}^{\mathrm{2}}\hspace{0.17em}\left(\frac{\mathrm{1}}{{\mathit{r̂}\mathrm{0}}^{\mathrm{2}\mathit{/}\mathit{\delta }}}\mathrm{-}\frac{\mathrm{1}}{{\mathit{r̂}}^{\mathrm{2}}}\right)]}^{\mathrm{1}\mathit{/}\mathrm{2}}\mathit{,}\end{array}$(62)which can be expressed without the W function.

2.6.5. Comparison with the β velocity law

Eq. (62) yields a terminal velocity (as ) of ${\mathit{v̂}}_{\mathrm{\infty }}\hspace{0.17em}\mathrm{\left(}\mathit{\theta }\mathrm{\right)}\mathrm{=}\sqrt{\frac{\mathrm{2}}{\mathit{r̂}\mathrm{0}}\hspace{0.17em}\left(\frac{\mathit{ĝ}\mathrm{0}}{\mathit{\delta }\left(\mathrm{1}\mathrm{+}\mathit{\gamma }\right)}\mathrm{-}{\mathit{v̂}}_{\mathrm{crit}}^{\mathrm{2}}\hspace{0.17em}{\mathit{r̂}\mathrm{0}}^{\mathrm{1}\mathrm{-}\mathrm{1}\mathit{/}\mathit{\delta }}\right)\mathrm{+}\frac{{\mathit{v̂}\mathrm{rot}}^{\mathrm{2}}}{{\mathit{r̂}\mathrm{0}}^{\mathrm{2}\mathit{/}\mathit{\delta }}}}\mathit{,}$(63)dependent on angle θ, which is now comparable to the parameter in the (so-called) β velocity law (cf. Castor & Lamers 1979; CAK) $\mathit{v̂}\mathit{r}\hspace{0.17em}\mathrm{\left(}\mathit{r̂}\mathrm{\right)}\mathrm{=}{\mathit{v̂}}_{\mathrm{\infty }}\hspace{0.17em}\left(\mathrm{1}\mathrm{-}\frac{\mathit{r̂}\begin{array}{c}\mathrm{\prime }\\ \mathrm{0}\end{array}}{\mathit{r̂}}\right)\beta$(64)for a given latitude of interest. To be able to compare the γ (and δ) parameter in our wind law with the exponent in the β law (as we use β as an input parameter in our model atmosphere calculations), we express our line acceleration parameter $\mathit{ĝ}{\mathrm{0}}_{}$ in terms of by means of Eq. (63), i.e. $\mathit{ĝ}\mathrm{0}\mathrm{=}\left(\left(\frac{{\mathit{v̂}}_{\mathrm{\infty }}^{\mathrm{2}}}{\mathrm{2}}\mathrm{-}\frac{{\mathit{v̂}\mathrm{rot}}^{\mathrm{2}}}{{\mathit{r̂}\mathrm{0}}^{\mathrm{2}\mathit{/}\mathit{\delta }\mathrm{+}\mathrm{1}}}\right)\mathrm{+}\frac{{\mathit{v̂}}_{\mathrm{crit}}^{\mathrm{2}}}{{\mathit{r̂}\mathrm{0}}^{\mathrm{1}\mathit{/}\mathit{\delta }}}\right)\hspace{0.17em}\mathit{r̂}\mathrm{0}\hspace{0.17em}\mathit{\delta }\left(\mathrm{1}\mathrm{+}\mathit{\gamma }\right)\mathit{,}$(65)and insert it into Eq. (62), as to obtain $\begin{array}{ccc}\mathit{v̂}\mathit{r}\hspace{0.17em}\mathrm{\left(}\mathit{r̂}\mathit{,\theta }\mathrm{\right)}& \mathrm{=}& [\frac{\mathrm{2}}{\mathit{r̂}\mathrm{0}}\hspace{0.17em}({\mathit{v̂}}_{\mathrm{crit}}^{\mathrm{2}}\left(\frac{\mathit{r̂}\mathrm{0}}{\mathit{r̂}}\mathrm{-}{\mathit{r̂}\mathrm{0}}^{\mathrm{1}\mathrm{-}\mathrm{1}\mathit{/}\mathit{\delta }}\right)\\ & & \mathrm{+}\left(\left(\frac{\mathit{r̂}\mathrm{0}}{\mathrm{2}}{\mathit{v̂}}_{\mathrm{\infty }}^{\mathrm{2}}\mathrm{-}\frac{{\mathit{v̂}\mathrm{rot}}^{\mathrm{2}}}{{\mathit{r̂}\mathrm{0}}^{\mathrm{2}\mathit{/}\mathit{\delta }}}\right)\mathrm{+}{\mathit{v̂}}_{\mathrm{crit}}^{\mathrm{2}}\hspace{0.17em}{\mathit{r̂}\mathrm{0}}^{\mathrm{1}\mathrm{-}\mathrm{1}\mathit{/}\mathit{\delta }}\right){\left(\mathrm{1}\mathrm{-}\frac{\mathit{r̂}\mathrm{0}}{{\mathit{r̂}}^{\mathit{\delta }}}\right)}^{\mathrm{1}\mathrm{+}\mathit{\gamma }}{)}^{}\\ & & {\mathrm{+}{\mathit{v̂}\mathrm{rot}}^{\mathrm{2}}\hspace{0.17em}\left(\frac{\mathrm{1}}{{\mathit{r̂}\mathrm{0}}^{\mathrm{2}\mathit{/}\mathit{\delta }}}\mathrm{-}\frac{\mathrm{1}}{{\mathit{r̂}}^{\mathrm{2}}}\right)]}^{\mathrm{1}\mathit{/}\mathrm{2}}\mathit{,}\end{array}$(66)which now depends on , , γ   (θ), δ   (θ) and .

Since δ is of the order of 1, as is $\mathit{r̂}{\mathrm{0}}_{}$, one can approximate the expression ${\mathit{r̂}{\mathrm{0}}_{}}^{\mathrm{1}\mathrm{-}\mathrm{1}\mathit{/}\mathit{\delta }}$, in Eq. (66), as 1. Furthermore, as our line parameter $\mathit{r̂}{\mathrm{0}}_{}$ is defined as the parameter for which the line acceleration becomes zero at radius (cf. Eq. (39)), this radius is very close to the radius ${\mathit{r̂}{\mathrm{0}}_{}}^{\mathrm{\prime }}$ in the β–law (in Eq. (64)), where the wind velocity is assumed to be zero, i.e. .

Then, we can set the velocity in Eq. (64) equal to our velocity law in Eq. (66), to search for a relationship between the parameters β, γ and δ, which yields (for ) ${\left(\mathrm{1}\mathrm{-}\frac{\mathit{r̂}\mathrm{0}}{\mathit{r̂}}\right)}^{\mathrm{2}\mathit{\beta }\mathrm{-}\mathrm{1}}\stackrel{\mathrm{!}}{\mathrm{\approx }}\mathrm{-}\mathit{b}\mathrm{+}\mathrm{(}\mathrm{1}\mathrm{+}{\mathit{b}}^{\mathrm{)}}\hspace{0.17em}\frac{{\left(\mathrm{1}\mathrm{-}\frac{\mathit{r̂}\mathrm{0}}{{\mathit{r̂}}^{\mathit{\delta }}}\right)}^{\mathrm{1}\hspace{0.17em}\mathrm{+}\hspace{0.17em}\mathit{\gamma }}}{\left(\mathrm{1}\mathrm{-}\frac{\mathit{r̂}\mathrm{0}}{\mathit{r̂}}\right)}$(67)with $\begin{array}{ccc}\mathit{b}\mathrm{:}\mathrm{=}\frac{\mathrm{2}}{\mathit{r̂}\mathrm{0}}{\left(\frac{{\mathit{v̂}}_{\mathrm{crit}}}{{\mathit{v̂}}_{\mathrm{\infty }}}\right)}^{\mathrm{2}}\mathit{,}& & \end{array}$analogous to the one-dimensional case of a spherical wind (cf. Paper I).

Herein, for large radii (and especially for small values of b, e.g. b ≲ 0.1 for an O-V-star), the right hand side can be approximated by the last fraction only, which leads to ${\left(\mathrm{1}\mathrm{-}\frac{\mathit{r̂}\mathrm{0}}{\mathit{r̂}}\right)}^{\mathrm{2}\mathit{\beta }}\stackrel{\mathrm{!}}{\mathrm{\approx }}{\left(\mathrm{1}\mathrm{-}\frac{\mathit{r̂}\mathrm{0}}{{\mathit{r̂}}^{\mathit{\delta }}}\right)}^{\mathrm{1}\hspace{0.17em}\mathrm{+}\hspace{0.17em}\mathit{\gamma }}\mathit{,}$(68)or equivalently $\frac{\mathrm{2}\mathit{\beta }}{\mathrm{1}\mathrm{+}\mathit{\gamma }}\mathrm{\approx }\frac{\log \left(\mathrm{1}\mathrm{-}\frac{\mathit{r̂}\mathrm{0}}{{\mathit{r̂}}^{\mathit{\delta }}}\right)}{\log \left(\mathrm{1}\mathrm{-}\frac{\mathit{r̂}\mathrm{0}}{\mathit{r̂}}\right)}\mathrm{·}$(69)Next, the right hand side in Eq. (69) can be approximately set to 1, for values of δ −  → 1 or generally for smaller distances from the central star in the supersonic region as . This results the same relationship between β and γ as for a non-rotating spherical wind (see Paper I) $\mathrm{2}\hspace{0.17em}\mathit{\beta }\hspace{0.17em}\mathrm{\left(}\mathit{\theta }\mathrm{\right)}\mathrm{\approx }\mathrm{1}\mathrm{+}\mathit{\gamma }\hspace{0.17em}\mathrm{\left(}\mathit{\theta }\mathrm{\right)}\mathit{,}$(70)independent of δ, that is valid for the previously mentioned values for δ at smaller radii . It also applies at very large distances , since then, the numerical values inside the brackets of Eq. (68) are close to 1 and this equation is fulfilled for any value of the exponents β and γ in all cases. Only for intermediate distances from the star at lower values of δ (not close to 1), the relationship between β and γ is possibly not well approximated by Eq. (70).

2.6.6. Fitting formula for the line acceleration

Thus, we can provide another expression for the line acceleration (in Eq. (39)), now dependent on and on the parameters , γ (or β equivalently), δ and $\mathit{r̂}{\mathrm{0}}_{}$ (by eliminating $\mathit{ĝ}{\mathrm{0}}_{}$ using Eq. (65)): $\begin{array}{ccc}\mathit{ĝ}\begin{array}{c}\mathrm{line}\\ \mathrm{rad}\end{array}\mathrm{\left(}\mathit{r̂}\mathit{,\theta }\mathrm{\right)}\mathrm{=}\left(\left(\frac{{\mathit{v̂}}_{\mathrm{\infty }}^{\mathrm{2}}}{\mathrm{2}}\mathrm{-}\frac{{\mathit{v̂}\mathrm{rot}}^{\mathrm{2}}}{{\mathit{r̂}\mathrm{0}}^{\mathrm{2}\mathit{/}\mathit{\delta }\mathrm{+}\mathrm{1}}}\right)\mathrm{+}\frac{{\mathit{v̂}}_{\mathrm{crit}}^{\mathrm{2}}}{{\mathit{r̂}\mathrm{0}}^{\mathrm{1}\mathit{/}\mathit{\delta }}}\right)\mathit{r̂}\mathrm{0}\hspace{0.17em}\frac{\mathit{\delta }\left(\mathrm{1}\mathrm{+}\mathit{\gamma }\right)}{\mathit{r̂}\mathrm{1}\hspace{0.17em}\mathrm{+}\hspace{0.17em}\mathit{\delta }\hspace{0.17em}\left(\mathrm{1}\hspace{0.17em}\mathrm{+}\hspace{0.17em}\mathit{\gamma }\right)}\hspace{0.17em}{\mathrm{(}{\mathit{r̂}}^{\mathit{\delta }}\mathrm{-}{\mathit{r̂}\mathrm{0}}^{\mathrm{)}}}^{\mathit{\gamma }}\mathit{.}& & \end{array}$(71)This non-linear expression can then be used as fitting formula and applied to the results from a numerical calculation of $\mathit{ĝ}\begin{array}{c}\mathrm{line}\\ \mathrm{rad}\end{array}\hspace{0.17em}\mathrm{\left(}{\mathit{r̂}}_{\mathit{i}}\mathit{,\theta }\mathrm{\right)}$ for discrete radial grid points at a given latitude, in order to determine the line acceleration parameters γ   (θ) (or equivalently β   (θ) by means of Eq. (70)), δ   (θ), and the terminal velocity , cf. Fig. 1.

2.6.7. Physical interpretation of the equation of critical radius

Through the use of the exact wind solution, by using Eq. (42), valid for the critical radius at latitude with polar angle θ, we can solve for the line acceleration parameter ĝ₀ and insert it into Eq. (63) from the approximated wind solution, to provide another expression for the terminal velocity at θ, depending on the location of the critical (i.e. sonic) point $\begin{array}{ccc}{\mathit{v̂}}_{\mathrm{\infty }}\hspace{0.17em}\mathrm{\left(}\mathit{\theta }\mathrm{\right)}& \mathrm{=}& [\frac{\mathrm{2}}{\mathit{r̂}\mathrm{0}}\{{\left(\frac{{\mathit{r̂}}_{\mathrm{c}}^{\mathit{\delta }}}{{\mathit{r̂}}_{\mathrm{c}}^{\mathit{\delta }}\mathrm{-}\mathit{r̂}\mathrm{0}}\right)}^{\mathit{\gamma }}\frac{{\mathit{r̂}}_{\mathrm{c}}^{\mathit{\delta }\hspace{0.17em}\mathrm{-}\hspace{0.17em}\mathrm{2}}}{\mathit{\delta }\left(\mathrm{1}\mathrm{+}\mathit{\gamma }\right)}\mathrm{(}{\mathit{v̂}}_{\mathrm{crit}}^{\mathrm{2}}\hspace{0.17em}{\mathit{r̂}}_{\mathrm{c}}\mathrm{-}\mathrm{2}\hspace{0.17em}{\mathit{r̂}}_{\mathrm{c}}^{\mathrm{2}}\mathrm{-}{\mathit{v̂}\mathrm{rot}}^{\mathrm{2}}\mathrm{)}\\ & & { \mathrm{-}{\mathit{v̂}}_{\mathrm{crit}}^{\mathrm{2}}{\mathit{r̂}\mathrm{0}}^{\mathrm{1}\hspace{0.17em}\mathrm{-}\hspace{0.17em}\mathrm{1}\mathit{/}\mathit{\delta }}\mathrm{\}}\mathrm{+}\frac{{\mathit{v̂}\mathrm{rot}}^{\mathrm{2}}}{{\mathit{r̂}\mathrm{0}}^{\mathrm{2}\mathit{/}\mathit{\delta }}}]}^{\mathrm{1}\mathit{/}\mathrm{2}}\mathrm{·}\end{array}$(72)Or vice versa, by Eq. (42), the location of the critical point (through which the exact analytical wind solution of our Eq. of motion EOM (40) passes) is mainly determined, on the one hand, by the given terminal velocity v_∞, via the line acceleration parameter $\mathit{ĝ}{\mathrm{0}}_{}$. On the other hand, the position of must also be dependent on the given minimum velocity v_in   (θ) at the inner boundary radius R_in   (θ), where the velocity solution passes. This inner velocity v_in follows indirectly from the other remaining line acceleration or wind parameters γ, δ (which make up the shape of the velocity curve) and especially $\mathit{r̂}{\mathrm{0}}_{}$ (where the value of the latter parameter is determined by the radius ${\mathit{r̂}{\mathrm{0}}_{}}^{\mathrm{1}\mathit{/}\mathit{\delta }}$ at which ${\mathit{g}}_{\mathrm{rad}}^{\mathrm{line}}$ is zero).

Since the inner boundary condition of the velocity v_in   (θ) is connected to the mass-loss rate Ṁ   (θ), through the equation of mass continuity by Eq. (13) and the given density at the inner boundary, the position of the critical radius is uniquely specified by the values of v_∞   (θ) and Ṁ   (θ).

3. Numerical methods

3.1. Computing the radiative acceleration

As in Paper I, we first calculate the thermal, density and ionisation structure of the wind model by means of the non-LTE expanding atmosphere (improved Sobolev approximation) code ISA-Wind (de Koter et al. 1993, 1997). As a next step, we calculate the radiative acceleration as a function of distance by means of a Monte Carlo (MC) code MC-Wind (de Koter et al. 1997; Vink et al. 1999), accounting for the possibility that the photons can be scattered or eliminated (if they are scattered back into the star). The radiative transfer in MC-Wind involves multiple continuum and line processes using the Sobolev approximation (cf. Mazzali & Lucy 1993).

The radiative acceleration of the wind at a given constant co-latitude θ is calculated by following the fate of the large number of photons where the atmosphere is divided into a large number of concentric thin shells with radius r and thickness dr, and the loss of photon energy, due to all scatterings that occur within each shell, is determined. The total line acceleration per shell summed over all line scatterings in that shell equals (Abbott & Lucy 1985) ${\mathit{g}}_{\mathrm{rad}}^{\mathrm{line}}\hspace{0.17em}\mathrm{\left(}\mathit{r,\theta }\mathrm{\right)}\mathrm{=}\mathrm{-}\frac{\mathrm{1}}{\mathit{Ṁ}\hspace{0.17em}\mathrm{\left(}\mathit{\theta }\mathrm{\right)}}\frac{\mathrm{d}\mathit{L}}{\mathrm{d}\mathit{r}}\hspace{0.17em}\mathrm{\left(}\mathit{r,\theta }\mathrm{\right)}\mathit{,}$(73)here, (in contrast to the one-dimensional case) referred to a constant polar angle θ, where −dL   (r,θ) is the rate at which the radiation field loses energy by the transfer of momentum of the photons to the ions of the wind per time interval.

The line list that is used for the MC calculations consists of over 10⁵ of the strongest lines of the elements from H to Zn from a line list constructed by Kurucz & Bell (1995). Lines in the wavelength region between 50 and 10 000 Å are included with ionisation stages up to stage VII. The number of photon packets distributed over the spectrum in our wind model, followed from the lower boundary of the atmosphere, is 2–3.5 × 10⁷. The wind is divided into 90 spherical shells with a large number of narrow shells in the subsonic region and wider shells in the supersonic range.

3.2. Computing the mass-loss rate for known stellar and wind parameters

By neglecting the pressure term and using the expression for the line acceleration per shell (Eq. (73)), an integration of the Eq. of motion (35) at a given constant latitude, from stellar radius to infinity, yields (cf. Abbott & Lucy 1985) $\begin{array}{ccc}\frac{\mathrm{1}}{\mathrm{2}}\hspace{0.17em}\mathit{Ṁ}\hspace{0.17em}\mathrm{\left(}\mathit{\theta }\mathrm{\right)}\hspace{0.17em}\mathrm{(}{\mathit{v}}_{\mathrm{\infty }}^{\mathrm{2}}\hspace{0.17em}\mathrm{\left(}\mathit{\theta }\mathrm{\right)}\mathrm{+}{\mathit{v}}_{\mathrm{esc}}^{\mathrm{2}}\hspace{0.17em}\mathrm{(}\mathit{\theta }{\mathrm{)}}^{\mathrm{)}}\mathrm{=}\mathrm{\Delta }\mathit{L}\hspace{0.17em}\mathrm{\left(}\mathit{\theta }\mathrm{\right)}\mathit{,}& & \end{array}$or equivalently, $\mathit{Ṁ}\hspace{0.17em}\mathrm{\left(}\mathit{\theta }\mathrm{\right)}\mathrm{=}\frac{\mathrm{2}\hspace{0.17em}\mathrm{\Delta }\hspace{0.17em}\mathit{L}\hspace{0.17em}\mathrm{\left(}\mathit{\theta }\mathrm{\right)}}{{\mathit{v}}_{\mathrm{\infty }}^{\mathrm{2}}\hspace{0.17em}\mathrm{\left(}\mathit{\theta }\mathrm{\right)}\mathrm{+}{\mathit{v}}_{\mathrm{esc}}^{\mathrm{2}}\hspace{0.17em}\mathrm{\left(}\mathit{\theta }\mathrm{\right)}}$(74)(as in Paper I, but here depending on θ), where ΔL   (θ) is the total amount of radiative energy extracted per second, summed over all the shells (at co-latitude θ). This equation is now fundamental for determining mass-loss rates numerically from the total removed radiative luminosity, for the prespecified stellar and wind parameters v_esc   (θ) and v_∞   (θ), respectively.

3.3. The iteration method: determination of $\dot{\mathit{M}}\mathrm{\left(}\mathit{\theta }\mathrm{\right)}$ and the wind parameters

To compute the mass-loss rate and the wind model parameters from a given rotating central star with the fixed stellar parameters⁶L   (θ), T_eff   (θ), R   (θ), M, Γ   (θ), V_rot, the analogous iterative procedure as in the one-dimensional non-rotating case can be applied (cf. Paper I). However, the whole iteration cycle here has to be performed separately for each given co-latitude θ of interest:

• 1.

  By keeping the stellar and wind parametersṀ_n(θ), v_∞n(θ), β_n(θ) variable throughout our iteration process, we use arbitrary (but reasonable) starting values Ṁ_-1, v_∞-1, β_-1 in iteration step n =  −1 (cf. Tables A.1–A.4).

• 2.

  For these input parameters, a model atmosphere is calculated with ISA-Wind for constant co-latitude θ. The code yields the thermal structure, the ionisation and excitation structure, and the population of energy levels of all relevant ions. Then, the radiative acceleration ${\mathit{g}}_{\mathrm{rad}}^{\mathrm{line}}\hspace{0.17em}\mathrm{\left(}{\mathit{r}}_{\mathit{i}}\mathit{,\theta }\mathrm{\right)}$ is calculated for discrete radial grid points r_i at polar angle θ with MC-Wind and Eq. (73). In addition, an improved estimate for the mass-loss rate $\mathit{Ṁ}\begin{array}{c}\mathrm{out}\\ \mathit{n}\end{array}\mathrm{\left(}\mathit{\theta }\mathrm{\right)}$ is obtained by Eq. (74), which can be used as a new input value for the next iteration step. Moreover, one obtains a new output value for the sonic radius (which has to be equal to the critical radius of our wind theory).

• 3.

  To determine the improved line acceleration parameters γ_n(θ) (or equivalently β_n(θ)), δ_n(θ) and for the considered latitude, Eq. (71) (together with Eq. (70)) is used as the fitting formula to apply to the numerical results for ${\mathit{g}}_{\mathrm{rad}}^{\mathrm{line}}\hspace{0.17em}\mathrm{\left(}{\mathit{r}}_{\mathit{i}}\mathit{,\theta }\mathrm{\right)}$, cf. Fig. 1.

• 4.

  By applying Eq. (72) and inserting the current values of parameters γ_n, δ_n and $\mathit{r̂}{{\mathrm{0}}_{\mathit{n}}}_{}$, as well as the current sonic radius for , we obtain a new approximation of the terminal velocity v_∞n(θ), i.e. $\begin{array}{ccc}{\mathit{v}}_{{\mathrm{\infty }}_{\mathit{n}}}\mathrm{\left(}\mathit{\theta }\mathrm{\right)}& \mathrm{=}& \mathit{a}[\frac{\mathrm{2}}{\mathit{r̂}{\mathrm{0}}_{\mathit{n}}}({\left(\frac{{\mathit{r̂}}_{{\mathrm{s}}_{\mathit{n}}}^{{\mathit{\delta }}_{\mathit{n}}}}{{\mathit{r̂}}_{{\mathrm{s}}_{\mathit{n}}}^{{\mathit{\delta }}_{\mathit{n}}}\mathrm{-}\mathit{r̂}{\mathrm{0}}_{\mathit{n}}}\right)}^{{\mathit{\gamma }}_{\mathit{n}}}\frac{{\mathit{r̂}}_{{\mathrm{s}}_{\mathit{n}}}^{{\mathit{\delta }}_{\mathit{n}}\mathrm{-}\mathrm{2}}}{{\mathit{\delta }}_{\mathit{n}}\left(\mathrm{1}\mathrm{+}{\mathit{\gamma }}_{\mathit{n}}\right)}\mathrm{(}{\mathit{v̂}}_{\mathrm{crit}}^{\mathrm{2}}\hspace{0.17em}{\mathit{r̂}}_{{\mathrm{s}}_{\mathit{n}}}\mathrm{-}\mathrm{2}\hspace{0.17em}{\mathit{r̂}}_{{\mathrm{s}}_{\mathit{n}}}^{\mathrm{2}}\mathrm{-}{{\mathit{v̂}}_{\mathrm{rot}}^{\mathrm{2}}}^{\mathrm{)}}\\ & & {\mathrm{-}{\mathit{v̂}}_{\mathrm{crit}}^{\mathrm{2}}{\mathit{r̂}{\mathrm{0}}_{\mathit{n}}}^{\mathrm{1}\mathrm{-}\mathrm{1}\mathit{/}{\mathit{\delta }}_{\mathit{n}}})\mathrm{+}\frac{{\mathit{v̂}\mathrm{rot}}^{\mathrm{2}}}{{\mathit{r̂}{\mathrm{0}}_{\mathit{n}}}^{\mathrm{2}\mathit{/}{\mathit{\delta }}_{\mathit{n}}}}]}^{\mathrm{1}\mathit{/}\mathrm{2}}\mathrm{·}\end{array}$(75)

• 5.

  With these improved estimates of Ṁ_n(θ), v_∞n(θ), β_n(θ) as new input parameters, the whole iteration step, defined by items 2–4, is repeated until convergence (at given latitude) is achieved.

3.4. The adjustment of the wind formalism to ISA-Wind

The ISA-Wind code has already been described in detail by de Koter et al. (1993, 1996). Those model assumptions within ISA-Wind which affect our wind formalism, by using the code in our iteration process, have also been described in Paper I.

To be able to apply the analytical wind solution of our EOM (40) to model a stellar wind from a given central star (with fixed stellar parameters) by using ISA-Wind to find numerically the unique solution, we had to adjust our wind formalism, i.e. our more accurate EOM, to the assumed EOM and wind velocity structure in ISA-Wind. The EOM (40) (and therefore our exact analytical wind solution) considers (allows) a line acceleration throughout the whole wind regime, starting above the radius ${\mathit{r̂}}_{\mathrm{0}}^{\mathrm{1}\mathit{/}\mathit{\delta }}$, whereas the different EOM in ISA-Wind is only solved in the subsonic region by neglecting the line force. However then, ISA-Wind “switches on” the line force somewhere below a connection radius in the subsonic region by assuming a β velocity law above in the supersonic region.

This inconsistency through the use of ISA-Wind (compared to our model assumptions and solutions) has been eliminated by introducing the parameter into our formalism of Sect. 2 for which the line radiation force is zero at radius ${\mathit{r̂}}_{\mathrm{0}}^{\mathrm{1}\mathit{/}\mathit{\delta }}$ at a given co-latitude θ. Then, the final value of , together with the other remaining line acceleration or wind parameters, can be determined by fitting Eq. (71) and the iteration procedure.

Moreover, a further inconsistency would occur if we applied ISA-Wind, developed for non-rotating spherical winds, in our iteration method to compute the wind parameters from a rotating star: the EOM in ISA-Wind for the subsonic region and the assumed β velocity law for the supersonic wind region, therein, do not consider the additional centrifugal term in contrast to our more accurate EOM (40), which we solve for the radial wind velocity from rotating stars. Therefore, to compensate for this discrepancy, one has to neglect the centrifugal term in our simplified EOM (59) for the supersonic wind region, which results a simpler approximated wind solution and terminal velocity v_∞   (θ) than that in Eq. (62) and Eq. (63), respectively, where $\mathit{v̂}{\mathrm{rot}}_{}$ can be set to zero. Then, this in turn, leads also to a simpler expression for the fitting formula for the line acceleration where the $\mathit{v̂}{\mathrm{rot}}_{}$ term in Eq. (71) vanishes. However, the derivation of the expression for the terminal velocity v_∞   (θ), as a function of critical radius (or sonic radius ) in Sect. 2.6.7, yields then almost the same Eq. (72), where only the last -term (from the simplified approximated wind solution) vanishes but not the ${\mathit{v̂}{\mathrm{rot}}_{}}^{\mathrm{2}}$-term that originates from our equation of critical radius, Eq. (42). The latter takes the stellar rotational speed at co-latitude θ and its influence on the location of the critical point into consideration.

To avoid this inconsistency in our subsequent wind models for rotating O-stars in Sect. 4, we have thus applied the following simplified fitting formula $\mathit{ĝ}\begin{array}{c}\mathrm{line}\\ \mathrm{rad}\end{array}\mathrm{\left(}\mathit{r̂}\mathit{,\theta }\mathrm{\right)}\mathrm{=}\left(\frac{{\mathit{v̂}}_{\mathrm{\infty }}^{\mathrm{2}}}{\mathrm{2}}\mathrm{+}\frac{{\mathit{v̂}}_{\mathrm{crit}}^{\mathrm{2}}}{{\mathit{r̂}\mathrm{0}}^{\mathrm{1}\mathit{/}\mathit{\delta }}}\right)\mathit{r̂}\mathrm{0}\hspace{0.17em}\frac{\mathit{\delta }\left(\mathrm{1}\mathrm{+}\mathit{\gamma }\right)}{\mathit{r̂}\mathrm{1}\mathrm{+}\mathit{\delta }\hspace{0.17em}\left(\mathrm{1}\mathrm{+}\mathit{\gamma }\right)}\hspace{0.17em}\mathrm{(}{\mathit{r̂}}^{\mathit{\delta }}\mathrm{-}{\mathit{r̂}\mathrm{0}}^{\mathrm{)}}\gamma$(76)(instead of Eq. (71)) to determine gradually the improved line acceleration parameters by our iterative procedure for a given constant co-latitude θ, whereas the new estimate of the terminal velocity for the next step has been determined by the simplified iteration formula $\begin{array}{ccc}{\mathit{v}}_{{\mathrm{\infty }}_{\mathit{n}}}\mathrm{\left(}\mathit{\theta }\mathrm{\right)}& \mathrm{=}& \mathit{a}[\frac{\mathrm{2}}{\mathit{r̂}{\mathrm{0}}_{\mathit{n}}}({\left(\frac{{\mathit{r̂}}_{{\mathrm{s}}_{\mathit{n}}}^{{\mathit{\delta }}_{\mathit{n}}}}{{\mathit{r̂}}_{{\mathrm{s}}_{\mathit{n}}}^{{\mathit{\delta }}_{\mathit{n}}}\mathrm{-}\mathit{r̂}{\mathrm{0}}_{\mathit{n}}}\right)}^{{\mathit{\gamma }}_{\mathit{n}}}\frac{{\mathit{r̂}}_{{\mathrm{s}}_{\mathit{n}}}^{{\mathit{\delta }}_{\mathit{n}}\mathrm{-}\mathrm{2}}}{{\mathit{\delta }}_{\mathit{n}}\left(\mathrm{1}\mathrm{+}{\mathit{\gamma }}_{\mathit{n}}\right)}\mathrm{(}{\mathit{v̂}}_{\mathrm{crit}}^{\mathrm{2}}\hspace{0.17em}{\mathit{r̂}}_{{\mathrm{s}}_{\mathit{n}}}\mathrm{-}\hspace{0.17em}\mathrm{2}\hspace{0.17em}{\mathit{r̂}}_{{\mathrm{s}}_{\mathit{n}}}^{\mathrm{2}}\mathrm{-}{{\mathit{v̂}}_{\mathrm{rot}}^{\mathrm{2}}}^{\mathrm{)}}\\ & & {\mathrm{-}{\mathit{v̂}}_{\mathrm{crit}}^{\mathrm{2}}{\mathit{r̂}{\mathrm{0}}_{\mathit{n}}}^{\mathrm{1}\mathrm{-}\mathrm{1}\mathit{/}{\mathit{\delta }}_{\mathit{n}}})]}^{\mathrm{1}\mathit{/}\mathrm{2}}\mathit{,}\end{array}$(77)instead of the generally more accurate Eq. (75).

Further, since ISA-Wind begins its computations already below (however close to) the stellar (i.e. photospheric) radius, all formulae derived in Sect. 2 have been applied with reference to the inner boundary (core) radius R_in   (θ) from where the numerical calculations of the wind model start. Therefore, the dimensionless variable of distance at polar angle θ (in Sect. 4) refers to R = R_in   (θ).

3.5. The chosen boundary values in the wind models

In our following numerical wind models the inner boundary radius has been chosen constant throughout the whole iteration process for each chosen latitude and star. E.g., for the particular case of the wind from the undistorted rotating O-main-sequence star (cf. Table 1), it is even constantly R_in =11.757   R_⊙ at each co-latitude θ (e.g. at 0 and π/2), situated at a prescribed fixed Rosseland optical depth of about τ_R = 23. This corresponds then to a photospheric radius of R_phot = 11.828   R_⊙ (at each latitude), defined as where the thermal optical depth is ${\mathit{\tau }}_{\mathrm{e}}\mathrm{=}\mathrm{1}\mathit{/}\sqrt{\mathrm{3}}$, and an inner boundary density of ρ_in = 1.398 × 10^-8 g/cm³ at R_in.

Small changes of this particular chosen fixed value for τ_R, or corresponding ρ_in   (θ), at each step of the iteration cycle (generally dependent on the given co-latitude θ and stellar rotation V_rot in the case of a distorted star, cf. Table 3) would have no effect on the final wind parameters, i.e. in particular on the converged values of Ṁ   (θ) and v_∞   (θ), as already explained in our 1D Paper I.

4. Application: results for rotating O-stars

In this section we apply our theoretical results from Sect. 2 and the iterative procedure described in Sect. 3.3 to first compute the stellar wind parameters for a differentially rotating 40 M_⊙ O5–V main-sequence star with an equatorial rotation speed of V_rot = 300 km s^-1 (Ω = 0.42) and 500 km s^-1 (Ω = 0.70), respectively. Secondly, we determine the wind parameters for a 60 M_⊙ O-giant star rotating with V_rot = 300 km s^-1 (Ω = 0.55). In the first instance, we ignore the effects of stellar distortion and gravity darkening in order to be able to compare our procedures against previous models, such as those from FA, before we include the more realistic aspects involving stellar distortion in Sect. 4.2.

4.1. Rotating O-stars without distortion

The fixed parameters (see e.g. Martins et al. 2005 and references therein) of the O-main-sequence star are given in the upper part of Table 1, whereas the fixed stellar parameters of the 60 M_⊙ O-giant are displayed in the upper part of Table 2. The iteration steps are provided in the Appendix. The final results are listed in the lower part of Tables 1 and 2 respectively, and shown graphically for the OV star in Fig. 3. It can be noted that the predicted mass-loss has barely changed from the spherical case for the model rotating with 300 km s^-1. For the more rapidly rotating 500 km s^-1 model, the equatorial mass-loss rate increases slightly, by ~0.1 dex, in rough agreement with previous studies, such as FA. Furthermore, the terminal wind velocities also remain relatively constant, although they drop by up to 15% for the most rapid models.

The ratio of equatorial-to-polar density (ρ_eq/ρ_p) reaches a constant value at a sufficient distance from the stellar surface, which follows from Eq. (15) ${\left(\frac{{\mathit{\rho }}_{\mathrm{eq}}}{{\mathit{\rho }}_{\mathrm{p}}}\right)}_{\mathrm{\infty }}\mathrm{=}\frac{\mathit{Ṁ}\mathrm{eq}}{\mathit{Ṁ}\mathrm{p}}\hspace{0.17em}\frac{{\mathit{v}}_{{\mathrm{\infty }}_{\mathrm{p}}}}{{\mathit{v}}_{{\mathrm{\infty }}_{\mathrm{eq}}}}$(78)as illustrated, e.g., in the right diagram of Fig. 3 (lower panel). We therefore provide this ratio at infinity in all Tables 1–3.

For the rapid rotator (500 km s^-1) with the slightly enhanced equatorial mass-loss rate, we find a density contrast ratio of equator to the pole that varies with stellar distance, has a maximum of about 1.7 close to the star, and then approaches a constant value of ~1.5.

Next we turn our attention to the 60 M_⊙ O-giant rotating with V_rot = 300 km s^-1. Here we find the effects of rotation to be more pronounced than for the main sequence dwarf. Already at velocities as low as 300 km s^-1 the mass-loss rate at the equator is more than a factor two larger than for the spherical non-rotating case, and the equator-to-pole density contrast ratio is (ρ_eq/ρ_p)_∞ ≃ 2.4.

4.2. Considering oblateness and gravity darkening of the rotating O-main sequence star

We now turn our attention to the physically more realistic cases, including the effects of stellar oblateness and gravity darkening. The model parameter and results are listed in Table 3 and the iteration cycles are shown in the Appendix. Again the results are also shown graphically (in Fig. 4). The situation is notably different from the non-distorted case. Whilst the results for the 300 km s^-1 (Ω = 0.70) main-sequence dwarf are inconspicuous as they are not all that different from the spherical case, the behaviour of (ρ_eq/ρ_p) again varies somewhat with stellar distance, having its maximum of about 1.25 at about r/R =3.0, and then very slowly decreasing to approach a value of 0.97 at infinity. The more rapid rotator (at 500 km s^-1 and Ω = 0.92) shows, seemingly surprisingly, a lower mass-loss rate at the equator than at the pole, by almost an order of magnitude.

The right bottom panel of Fig. 4 also shows the opposite effect: the ratio of equatorial-to-polar density (ρ_eq/ρ_p) is significantly below 1, by about a factor of 5. In other words, the predictions including gravity darkening suggest a polarward rather than equatorial stellar wind. How can we understand this result? As can be noted from Eq. (33) the Eddington parameter has become modified to correct for the latitudinal surface distortion. In comparison to the 1D spherical case the luminosity L and effective temperature T_eff are lower and using the results of Vink et al. (2000) one can already expect the lowered equatorial mass-loss rate in comparison to the spherical case primarily due to the lower L. These results suggest that for rapid rotators, the big difference between the pole and equator would lead to significantly different diagnostics for a star that is observed from the pole or equator (see e.g. Hillier et al. 2003; Herrero et al. 2012). On a positive note, the pole-to-equator ratio of order 5 is sufficiently large to be measurable with linear polarisation data.

5. Discussion and conclusions

5.1. Comparison with previous models

The earliest works on CAK-type mass-loss predictions including the effects of stellar rotation via the effective mass in the equation of motion (but without considering gravity darkening) such as those by FA and Pauldrach et al. (1986), Petrenz & Puls (2000) resulted in moderate mass-loss rate enhancements, which were subsequently included in evolutionary calculations (Langer 2012). Our results, with mass-loss rates increases by at most 0.1 dex, are in good agreement with these earlier works based on the CAK theory.

When gravity darkening is included, previous studies, such as those of Owocki et al. (1996), Maeder (1999) and Pelupessy et al. (2000) all found relatively modest changes for low rotational velocities, but they all favour polar ejection for high rotation rates in close proximity to the break-up limit (Ω = 1). The most-oft used formulae for mass-loss enhancement on the basis of the CAK theory is the one given by Maeder & Meynet (2000). Their famous equation (4.29) reads: $\frac{\mathit{Ṁ}\mathrm{\left(}\mathrm{\Omega }\mathrm{\right)}}{\mathit{Ṁ}\mathrm{(}\mathrm{\Omega }\mathrm{=}\mathrm{0}\mathrm{)}}\mathrm{=}\frac{\mathrm{(}\mathrm{1}\mathrm{-}\mathrm{\Gamma }{\mathrm{)}}^{\frac{\mathrm{1}}{\mathit{\alpha }}\mathrm{-}\mathrm{1}}}{{\left[\mathrm{1}\mathrm{-}\frac{{\mathrm{\Omega }}^{\mathrm{2}}}{\mathrm{2}\mathit{\pi G}{\mathit{\rho }}_{\mathrm{m}}}\mathrm{-}\mathrm{\Gamma }\right]}^{\frac{\mathrm{1}}{\mathit{\alpha }}\mathrm{-}\mathrm{1}}}\mathit{,}$(79)where Γ = L/L_Edd = κL/(4πcGM) is the Eddington factor (with κ the total opacity), and α the T_eff −dependent CAK force multiplier parameter. Whilst it is claimed that the mass-loss increase is due to rotation (Ω) the more relevant reason for the high mass-loss rate is in fact the high Γ factor (see more recent computations by Vink 2006; Gräfener & Hamann 2008; Vink et al. 2011).

Surprisingly, our 2D results that take gravity darkening into account predict a equatorial decrease in the mass-loss rate. This would also imply that the total (cf. Eq. (14)) mass-loss rate is now lower than for the spherical 1D case. This is in contradiction to the above Maeder & Meynet equation that is oftentimes used in stellar evolution calculations for rotating massive stars.

5.2. Comparison with observations

There are two complementary types of astronomical techniques that can provide meaningful geometric constraints on large-scale 2D wind structures. These are interferometry and linear spectropolarimetry. The current status is that whilst interferometry has been very successful in obtaining results on objects such as classical Be stars, B[e] supergiants, and LBVs, results for O stars are – to the best of our knowledge – not yet available due to their smaller radii.

Linear spectropolarimetry can be employed with at least two levels of complexity. The first involves detailed line-profile predictions that can be used to determine the geometric and kinematic environments around young and massive stars (Harries 2000; Vink et al. 2005), as well as the polarimetric agents (Kuhn et al. 2007). The lower-level complexity application of linear spectropolarimetry is adopted here, with the emission line reasonably assumed to be unpolarised, but with the continuum responsible for the observed polarisation level: depolarisation (Poeckert & Marlborough 1976; Brown & McLean 1977). This latter approach is sufficient for our current comparisons. In this case a pole-to-equator density contrast of 5 might lead to a continuum polarisation of ≃1% according to the expectations of Harries et al. (1998).

If we assume that our most realistic density-contrast predictions concern those models that include gravity darkening, we would anticipate the vast majority of mostly moderately rotating O-dwarfs to be unpolarised, as the density contrast between the pole and equator is close to unity for rotation velocities up to 300 km s^-1. These overall findings are in good agreement with the linear Hα polarimetry survey of Vink et al. (2009).

An exception might involve the sub-group of Oe stars, for which Vink et al. (2009) list the highest rotation velocities, with vsini values approaching ≃400 km s^-1. As this concerns a lower limit (due to the inclination effects), the true rotational velocities of Oe stars might go all the way up to the break-up limit (Ω = 1). Interestingly, the Oe star HD 45314 was indeed shown to have a line effect, but the overall incidence of line effects in Oe stars was low: 1/6. On the basis of the general spectroscopic similarity between Oe and Be stars, and the presence of double-peaked line profiles many astronomers would expect Oe stars to be embedded in circumstellar disks, although it should be noted that spectroscopy alone cannot provide geometric constraints, without complementary imaging, interferometry, or polarimetry, because geometric information from spectral line profiles is not unique, due to the convolution that underlies spectral line formation. In particular, double-peaked profiles can even be noted in spherical shells (e.g. PP).

In the current paper the surprising result is that of polar outflows, whilst it should be noted that linear spectropolarimetry cannot yet distinguish between a prolate (polar) or oblate (equatorial) wind structure (e.g. Wood et al. 1997). It will be challenging to reconcile our axisymmetric density-contrast predictions with the general expectation that Oe/Be stars are embedded within equatorial structures. Clearly there are many puzzles remaining in the geometry and formation of circumstellar media around rotating OB stars, and quite possibly additional physical effects, such as magnetic fields, pulsations, and companions may need to be considered.

5.3. Summary and future work

In Paper I (Müller & Vink 2008), we proposed a new parametrisation of the radiative line acceleration, expressing it as a function of radius rather than of the velocity gradient (as in CAK theory) which generalised the classical thermal Parker (1958) wind and allowed solving it analytically. This has the advantage of higher accuracy and enables to check an influence of numerical viscosity if used as a test example for numerical simulations. In addition, the implementation of this formalism allows for local dynamical consistency as we are able to determine the momentum transfer at each location in the wind through the use of Monte Carlo simulations.

In this paper, we present a generalization of our model of spherical line-driven winds (in Paper I) to the case of axially symmetric rotating stars. We extend on the results of Paper I, deriving analytical solutions for the 2D case in an axisymmetric mass outflow or inflow scenario. The generalization approximates the solution of the two-dimensional problem of a rotating wind as a one-parametric set (parameterised by the latitude) of one-dimensional solutions. Here, the separation of the radial motion for individual latitudes is achieved by the basic assumption that the meridional flow velocity is negligible. Assuming furthermore only central external forces, then imply the conservation of angular momentum and leads to an additional centrifugal term in the equation of radial motion⁷. We also extend our iterative method for the simultaneous solution of the mass-loss rate and velocity field to the case of a rotating non-spherical stellar wind, using the parameterised description of the line acceleration that only depends on radius – at any given latitude.

Furthermore, an approximate analytical solution for the supersonic flow of a rotating wind is derived, which closely resembles the exact solution, and we apply the new expressions with our iterative method to (i) the stellar wind from a differentially rotating 40 M_⊙ O5–V main sequence star; as well as to (ii) a 60 M_⊙ O-giant star. We subsequently account for the effects of stellar oblateness and gravity darkening for the O5–V main sequence star. The results show an equatorial decrease in the mass-loss rate, which would imply a total mass-loss rate that is lower than for the spherical 1D case, in contradiction to the oft-used Maeder & Meynet (2000) formalism used in most current stellar evolution calculations for rotating massive stars.

In the future we plan to extend our method to B supergiants, in particular to tackle the problem of disk-formation in B[e] supergiants (Lamers & Pauldrach 1991; Pelupessy et al. 2000; Cure et al. 2005; Madura et al. 2007).

Online material

Appendix A: Tables of iteration cycles

---

Acknowledgments

We acknowledge financial report from the STFC and DCAL. We thank our anonymous referee for helpful comments.

References

All Tables

All Figures