3 Model equations

3.1 Basic assumptions

We assume that the stationary, spherically symmetric stellar wind consists of three components, namely absorbing ions, nonabsorbing hydrogen atoms and ions, and electrons, denoted by subscripts ${\rm i}$, ${\rm p}$, ${\rm e}$, respectively. Each of them is described by a density $\rho_a$, radial velocity ${ {v_{\rm r}}}_a$, temperature T_a, electrical charge q_a=e z_a(where e is an elementary charge and z_a denotes the ionization degree - may have a non-integer value), and particle mass m_a. Subscript a stands for $a={\rm i,p,e}$. Contrary to our previous models (KKI), we allow for different temperature of each component and for radial changes of electrical charge. We assume that chemical composition is given by the factor z_*, which is a stellar metallicity relative to the solar value.

3.2 Continuity equations

In the case of a stationary spherically symmetric stellar wind each component is described by the continuity Eq. (4) in the form of

$$\frac{1}{r^2}\frac{{\rm d}}{{\rm d} r}\left(r^2\rho_a{ {v_{\rm r}}}_a\right) = S_a,$$ (35a)

where the term S_a accounts for radial change of mass-loss rate of individual components due to the ionization. Whereas for all types of ions the mass-loss rate is constant through the wind and thus number of these particles is conserved ( ${S}_{{\rm i}}={S}_{{\rm p}} =0$), we account for the possibility of variation of electron number via ionization and recombination. Because the total electric charge is conserved,

$$\sum_a \frac{{\rm d}}{{\rm d} r}\left(r^2q_a\frac{\rho_a}{m_a}{ {v_{\rm r}}}_a\right)=0$$

and continuity equation holds separately for all components except electrons, we obtain electron continuity equation in the form of

$$\frac{1}{r^2}\frac{{\rm d}}{{\rm d} r}\left(r^2{\rho}_{{\rm e}}{{... ...{{m}_{{\rm e}}}{m_a} \rho_a{ {v_{\rm r}}}_a \frac{{\rm d} z_a}{{\rm d}r}\cdot$$ (35b)

Although inclusion of a term ${S}_{{\rm e}}$ into the electron continuity equation does not significantly alter the model, it is important to obtain well converged model.

3.3 Momentum equations

In the case of stationary spherically symmetric stellar wind the momentum equation Eq. (6) has the form of

 

$${ {v_{\rm r}}}_a\frac{{\rm d} { {v_{\rm r}}}_a}{{\rm d}r}= {g}_{a... ... {v_{\rm r}}}_b-{ {v_{\rm r}}}_a}{\vert{ {v_{\rm r}}}_b-{ {v_{\rm r}}}_a\vert},$$ (36)

where square of isothermal sound speed is a_a²=kT_a/m_a, E is a charge separation electric field. Gravitational acceleration has the form $g={G{\mbox{\eufont M}}}/{r^2}$, where $\mbox{\eufont M}$ is the stellar mass and G is the gravitational constant. The radiative acceleration acting on free electrons can be expressed as

$${g}_{{\rm e}}^{{\rm rad}}=\frac{{m}_{{\rm p}}}{{m}_{{\rm e}}}\frac{G\Gamma{\mbox{\eufont M}}}{r^2},$$ (37)

where the ratio of the radiative force caused by absorption of radiation by free electrons and gravitational force is

$$\Gamma=\frac{\sigma_{\rm e} L}{4\pi c G{\mbox{\eufont M}}}$$

L is the stellar luminosity and $\sigma_{\rm e}$ is the mass scattering coefficient of the free electrons (do not confuse it with Thomson scattering cross section).

The radiative acceleration acting on absorbing ions is taken in the form of Castor et al. (1975)

 

$${g}_{{\rm i}}^{{\rm rad}}= \frac{1}{{\mbox{\eufont Y}}_{{\rm i}}}... ...}_{{\rm i}}} \frac{{\rm d}{ {v_{\rm r}}}_{{\rm i}}}{{\rm d}r} \right)^{\alpha},$$ (38)

with force multipliers k, $\alpha$, $\delta$ after Abbott (1982). The finite disk correction factor (Friend & Abbott 1986; Pauldrach et al. 1986) is

$$f=\frac{\left(1+\sigma\right)^{\alpha+1}- \left(1+\sigma\mu_{\rm ... ...lpha+1\right)\left(1-\mu_{\rm c}^2\right)\sigma\left(1+\sigma\right)^{\alpha}},$$ (39)

where W is a dilution factor. The thermal speed $v_{{\rm th}}=\sqrt{2k{T}_{{\rm i}}/{m}_{{\rm p}}}$ is, owing to the Abbott's normalization, hydrogen thermal speed. However, because $v_{{\rm th}}$ physically describes ionic thermal speed, it depends on ionic temperature. Finally, ${\mbox{\eufont Y}}_{{\rm i}}$ is the photospheric ratio of the metallic ion density to the passive plasma density. We selected the same value as in KKI, namely ${\mbox{\eufont Y}}_{{\rm i}}=0.0127$ (this value corresponds to the solar ratio of sum of densities of C, N, O, Fe to the density of bulk plasma).

Constant of friction evaluated using Fokker-Planck approximation (cf. Burgers 1969) has the following form:

$${K}_{ab}={n}_a{n}_b\frac{4\pi {q}_a^2{q}_b^2}{k T_{ab}}\ln\Lambda,$$ (40)

where n_a and n_b are number densities of individual components and mean temperature of both components

$$T_{ab}=\frac{m_aT_b+m_bT_a}{m_a+m_b}\cdot$$ (41)

The Coulomb logarithm is of the form

$$\ln\Lambda= \ln\left[\frac{3k{T}_{{\rm e}}}{e^2}\left(\frac{k{T}_{{\rm e}}}{4\pi n e^2}\right)^{1/2}\right],$$ (42)

where n is the particle density ( $n={n}_{{\rm p}}+{n}_{{\rm e}}+{n}_{{\rm i}}$). Finally, the Chandrasekhar function G(x), defined in terms of the error function ${\rm erf}(x)$ (Dreicer 1959) is

$$G(x)=\frac{1}{2 x^2}\left({\rm erf}(x)-\frac{2x}{\sqrt\pi}\exp\left(-x^2\right)\right).$$ (43)

The argument of the Chandrasekhar function is

$$x_{ab}=\frac{\vert{ {v_{\rm r}}}_b-{ {v_{\rm r}}}_a\vert}{\alpha_{ab}},$$ (44)

where

$$\alpha_{ab}^2=\frac{2k\left(m_aT_b+m_bT_b\right)} {m_am_b}\cdot$$ (45)

3.4 Energy equation

Energy Eq. (9) in the case of a stationary, spherically symmetric multicomponent flow has the form of (cf. Burgers 1969)

 

$$\frac{3}{2}k{ {v_{\rm r}}}_a\frac{\rho_a}{m_a}\frac{{\rm d}T_a}{{\rm d}r} a_a^2\rho_a \frac{1}{r^2}\frac{{\rm d}}{{\rm d}r}\left(r^2 { {v_{\rm r}}}_a\right)= Q_a^{{\rm rad}}$$ +

$$\frac{1}{\sqrt\pi}\sum_{b\neq a}K_{ab} \frac{2k\left(T_b-T_a\right)}{m_a+m_b} \frac{\exp\left(-x_{ab}^2\right)}{\alpha_{ab}}$$ +

$$\sum_{b\neq a}\frac{m_b}{m_a+m_b}K_{ab}G(x_{ab})\vert{ {v_{\rm r}}}_b-{ {v_{\rm r}}}_a\vert.$$ + (46)

Two terms on the left-hand side stand for advection and adiabatic cooling, respectively. The right-hand side terms describe radiative heating/cooling, heat exchange by encounters of different particles caused by unequal temperatures of the components, and frictional heating.

There are two sources of radiative heating/cooling. First source are bound-free and free-free transitions and the second is Gayley-Owocki heating/cooling.

3.4.1 Radiative energy term of bound-free and free-free transitions

Bound-free and free-free transitions (which will be called "classical'' radiative transitions) deposit energy directly on electrons. Therefore, this classical radiative energy term should be considered in an electron energy equation. We decided to estimate the radiative heating/cooling term $Q^{{\rm rad}}$ using two mechanisms only, hydrogen Lyman bound-free and free-free transitions. The detailed form of heating and cooling in the above mentioned transitions is nearly the same as in KKI and will not be repeated here (see also Kubát et al. 1999). The only difference is that the temperature in these equations is now the electron temperature. $J_{\nu}$ at the base of the wind is taken as an emergent radiation from a spherically symmetric static hydrogen model atmosphere for a corresponding stellar type (Kubát 2001).

3.4.2 Gayley-Owocki heating/cooling

Contrary to bound-free and free-free transitions Gayley-Owocki heating/cooling deposits energy directly to absorbing ions. GO heating/cooling term has the following form (GO, Eq. (32))

 

$${Q}_{{\rm i}}^{{\rm rad}}\equiv{Q}_{{\rm i}}^{{\rm GO}} ={\rho}_{... ...}{\sigma_{\rm e}{\rho}_{{\rm i}}v_{{\rm th}}r}\right)^{\alpha} G(\sigma,\mu_*),$$ (47)

where the function $G(\sigma,\mu_*)$ is given by Eq. (33). Here $v_{{\rm th,i}}$ is a thermal speed of driving ions, thus, $v_{{\rm th,i}}=\sqrt{2k{T}_{{\rm i}}/{m}_{{\rm i}}}$. Inspecting the GO heating term (47) we conclude that wind is heated via the Doppler effect by direct radiation whereas the Gayley-Owocki cooling is caused by the diffuse radiation. The sign of GO heating depends on the value of a variable $\sigma$. In the case where $\sigma=0$, when the expansion is locally isotropic, the GO heating and cooling vanishes. When the $\sigma$ is positive, the cooling term dominates and the wind is cooled by the GO heating. In the opposite case, when the $\sigma$ is negative, the heating term dominates and thus, the wind is heated.

3.5 Charge separation electric field

The equation for charge separation electric field can be obtained directly from the third Maxwell equation, which in the case of spherical symmetry can be written as

 

$$\frac{1}{r^2}\frac{{\rm d}}{{\rm d}r}\left(r^2 E\right) =4\pi\sum_aq_an_a.$$ (48)

3.6 Ionic charge

The ionization structure of stellar wind should be derived using time consuming NLTE calculations (e.g. Pauldrach et al. 1994). Because we want to determine only a mean charge of selected elements, we use a simpler approximate method. As described by Mihalas (1978, Eq. (5.46) therein), the ionization equilibrium in stellar winds can be approximated by

 

$$\frac{n_{a,j}}{n_{a,j+1}}\approx \frac{1}{2} \frac{{n}_{{\rm e}}}... ...R}}{{T}_{{\rm e}}}\right)^{1/2}\exp\left(\frac{\chi_{a,j}}{k T_{\rm R}}\right),$$ (49)

where $\chi_{a,j}$ is the ionization potential, U_a,j is the partition function and $T_{\rm R}$ is the radiation temperature (we set $T_{\rm R}=\frac{3}{4}T_{\rm eff}$). Partition function approximations are taken from Smith & Dworetsky (1988). Electrical charge of absorbing and passive ions are then derived using a formula

$$z_a=\frac{\sum_j j n_{a,j}}{\sum_jn_{a,j}}\cdot$$ (50)

Clearly, electron electrical charge is ${z}_{{\rm e}}=-1$ everywhere.

3.7 Critical points

Critical points are points where derivatives of variables cannot be determined directly from differential equations. The derivation of critical point conditions for our set of equations is simpler than that of KKI because we use different temperatures of each component and correct momentum equation for electrons. Due to these generalizations, the set of critical point equations is not as coupled as it was in KKI and thus, the obtained critical point conditions are simpler.

We write model equations in a simplified form, where we explicitly write only terms containing derivatives of individual variables and other terms are included into the terms $F_{\rm i}$. Thus, the continuity Eqs. (35a), (35b) are

$$\rho_a\frac{{\rm d}{ {v_{\rm r}}}_a}{{\rm d} r}+{ {v_{\rm r}}}_a\frac{{\rm d}\rho_a}{{\rm d} r}= F_{a,1}(r,\rho_a,{ {v_{\rm r}}}_a).$$ (51a)

In the electron continuity Eq. (35b) we neglected the derivatives of ionic charge because their contribution to electron continuity equation is only marginal. However, inclusion of such term influences critical point and regularity conditions for electrons only, which will not be used (see bellow).

Similarly we can rewrite momentum Eq. (36). In the momentum equations of absorbing ions we shall linearize a term containing the velocity gradient. Note that because model equations are not quasi-linear (i.e. linear with respect to the derivatives of the independent variables), the mathematically more correct method would employ some form of transformation to the quasi-linear form (cf. Courant & Hilbert 1962). However, because the results are essentially the same in this case, we present analysis of critical points in a simplified form. Thus, momentum equations are

$${ {v_{\rm r}}}_a\frac{{\rm d} { {v_{\rm r}}}_a}{{\rm d}r}- \frac{... ...}{T_a}\frac{{\rm d}T_a}{{\rm d}r}= F_{a,2}(r,E,\rho_b, {v_{\rm r}}_b,T_b,z_b).$$ (51b)

Similarly, due to the dependence of the Doppler term on the velocity gradient (in the Sobolev approximation) we shall write energy equations in the form of

$$\frac{3}{2}{ {v_{\rm r}}}_a\rho_a\frac{{a}_a^2}{T_a}\frac{{\rm d}... ... { {v_{\rm r}}}_a}{{\rm d}r}= F_{a,3}(r,\rho_b, {v_{\rm r}}_b,T_b,z_{\rm b}).$$ (51c)

The system of equations is closed by the equation for charge separation electric field, which has a simple form,

$$\frac{{\rm d}E}{{\rm d}r}=F_{4}(r,\rho_b,z_b).$$(51d)

The system of Eq. (51) can be simplified by inserting the derivatives of density from the Eq. (51a) and derivatives of temperature (51c) into the momentum Eq. (51b). We obtain modified linearized momentum equations

 

$$\left[{ {v_{\rm r}}}_a- \frac{\partial {g}_{a}^{{\rm rad}}}{\part... ...{{\rm d} { {v_{\rm r}}}_a}{{\rm d}r}= F_{a}(r,E,\rho_b, {v_{\rm r}}_b,T_b,z_b),$$ (52)

where F_a is a combination of F_a,i from Eqs. (51a), (51b), and (51c). Because the charge separation electric field Eq. (51d) does not introduce any physically interesting critical point, the system of Eq. (52) consists of only three independent critical point equations. Each of them will be discussed separately in the following subsections.

3.7.1 Critical point of passive ions

For the passive plasma the critical point condition (52) has a simple form

 

$${{v_{\rm r}}^2}_{{\rm p}}=\frac{5}{3}{a}_{{\rm p}}^2.$$ (53)

In order to obtain continuous solution of model equations, at this point should the quantity $F_{a}(r,E,\rho_b, {v_{\rm r}}_b,T_b,z_b)$ vanish (regularity condition). This requirement is fulfilled if

 

$$\frac{10{a}_{{\rm p}}^2}{3r}-g+\frac{{q}_{{\rm p}}}{{m}_{{\rm p}}... ...b}-{ {v_{\rm r}}}_{\rm p}}{\vert{ {v_{\rm r}}}_{b}-{ {v_{\rm r}}}_{\rm p}\vert}$$

$$~~~~~-\frac{2}{3}\frac{1}{{\rho}_{{\rm p}}{{v_{\rm r}}}_{{\rm p}}... ...m_{{\rm p}}+m_b} \frac{\exp\left(-x_{{{\rm p}}b}^2\right)}{\alpha_{{{\rm p}}b}}$$

$$~~~~~- \frac{2}{3}\frac{1}{{\rho}_{{\rm p}}{{v_{\rm r}}}_{{\rm p}... ...}b}G(x_{{{\rm p}}b})\vert{ {v_{\rm r}}}_b-{ {v_{\rm r}}}_{{\rm p}}\vert=0. \,\,$$ (54)

This equation is a generalization of well-known regularity condition for the coronal wind.

3.7.2 Critical point of absorbing ions

Critical point condition Eq. (52) for absorbing ions has form

$${ {v_{\rm r}}}_{{\rm i}}- \frac{\partial {g}_{{\rm i}}^{{\rm rad}... ...^{{\rm rad}}}{\partial\left({\rm d}{{v_{\rm r}}}_{{\rm i}} /{\rm d}r\right)}=0.$$ (55)

However, as was discussed by KKI, this critical point condition is not met anywhere in the wind. Therefore, to obtain CAK type solution of the stellar wind we use another condition, which is connected to generalized Abbott waves (see KKI). Such a condition can be derived by multiplying the critical point conditions Eq. (52) by density of corresponding component and summing them

$$\sum_a \rho_a \frac{{\rm d} { {v_{\rm r}}}_a}{{\rm d}r} \left[{ {... ... i}} /{\rm d}r\right)}\right]= \sum_a F_{a}(r,E,\rho_b, {v_{\rm r}}_b,T_b,z_b).$$ (56)

To assure that common equation of motion (obtained by summing of individual momentum Eq. (36)) does not depend on the derivatives of any variables, a condition

 

$$\sum_a \rho_a \frac{{\rm d} { {v_{\rm r}}}_a}{{\rm d}r} \left[{ {... ... rad}}}{\partial\left({\rm d}{{v_{\rm r}}}_{{\rm i}} /{\rm d}r\right)}\right]=0$$ (57)

shall be fulfilled. This condition we use to fix the mass-loss rate.

3.7.3 Critical point of electrons

The last critical point condition for electrons has again a very simple form

 

$${ {v_{\rm r}}}_{{\rm e}}^2=\frac{5}{3}{a}_{{\rm e}}^2.$$ (58)

This critical point condition has a similar form to the condition of nonabsorbing ions Eq. (53). Thus, regularity condition for electrons resembles the regularity condition for passive ions Eq. (54),

 

$$\frac{10{a}_{{\rm e}}^2}{3r} {g}_{{\rm e}}^{{\rm rad}}-g+\frac{{q}_{{\rm e}}}{{m}_{{\rm e}}}E+... ...}}_b-{ {v_{\rm r}}}_{\rm e}}{\vert{ {v_{\rm r}}}_b-{ {v_{\rm r}}}_{\rm e}\vert}$$ +

$$\frac{2}{3}\frac{1}{{\rho}_{{\rm e}}{{v_{\rm r}}}_{{\rm e}} } \fr... ...m_{{\rm e}}+m_b} \frac{\exp\left(-x_{{{\rm e}}b}^2\right)}{\alpha_{{{\rm e}}b}}$$ -

$$\frac{2}{3}\frac{1}{{\rho}_{{\rm e}}{{v_{\rm r}}}_{{\rm e}} } \su... ..._{{{\rm e}}b}G(x_{{{\rm e}}b})\vert{ {v_{\rm r}}}_b-{ {v_{\rm r}}}_{\rm e}\vert$$ -

$$\frac{2}{3}\frac{1}{{\rho}_{{\rm e}}{{v_{\rm r}}}_{{\rm e}} }{Q}_{{\rm e}}^{{\rm rad}}=0.$$ - (59)

However, numerical tests showed that electron regularity condition (59) is approximately fulfilled at the electron critical point if the zero current condition is used as a boundary condition. Thus, this condition should not necessarily be included into the set of model equations.

3.8 Boundary conditions

3.8.1 Boundary conditions for temperatures

We assume that the flow at the inner boundary is in radiative equilibrium and that the boundary temperature of all components is the same, thus, we write boundary condition for temperatures in the form of

$$T_{a}={T}_{{\rm e}},\quad a={\rm i, p},$$ (60a)

$$Q^{{\rm rad}}(R_{*})=0.$$ (60b)

Boundary values of ionic charges can be directly obtained from the condition of ionization equilibrium (49).

3.8.2 Boundary condition for velocity

Conditions (54), (57) can be generally used to fix the boundary values of model quantities. However, inclusion of two inner conditions directly into model equations sometimes leads to numerical problems. Therefore, we use a more secure method, which gives essentially the same results.

We start to calculate our models at the passive plasma critical point. Consequently, the boundary condition for the passive plasma velocity is the critical point condition Eq. (53). Boundary condition for the velocity of absorbing ions may be obtained from the passive plasma regularity condition Eq. (54). Because we suppose equal boundary temperatures of each component Eq. (60a), the regularity condition may be simplified

 

$$\frac{10{a}_{{\rm p}}^2}{3r}-g+\frac{{q}_{{\rm p}}}{{m}_{{\rm p}}... ...b}-{ {v_{\rm r}}}_{\rm p}}{\vert{ {v_{\rm r}}}_{b}-{ {v_{\rm r}}}_{\rm p}\vert} ~~~~- \frac{2}{3}\frac{1}{{\rho}_{{\rm p}}{{v_{\rm r}}}_{{\rm p}}... ...{\rm p}}b}G(x_{{{\rm p}}b})\vert{ {v_{\rm r}}}_b-{ {v_{\rm r}}}_{\rm p}\vert=0.$$ (61)

The boundary value of electron velocity is chosen to fulfil the electron regularity condition Eq. (59) at the electron critical point Eq. (58). As was already mentioned, this condition is approximately satisfied if the zero current condition

 

$${n}_{{\rm e}}{ {v_{\rm r}}}_{{\rm e}}={z}_{{\rm p}}{n}_{{\rm p}}{ {v_{\rm r}}}_{{\rm p}} + {z}_{{\rm i}} {n}_{{\rm i}}{ {v_{\rm r}}}_{{\rm i}}$$ (62)

is used as boundary condition for electron velocity. The latter condition (62) was applied in our models.

3.8.3 Boundary conditions for density

We write the boundary condition for the passive plasma density in the same form as in KKI,

 

$$z_*{\rho}_{{\rm p}}(R_{*})= \frac{1}{{\mbox{\eufont Y}}_{{\rm i}}... ...}) \frac{{{v_{\rm r}}}_{{\rm i}} (R_{*})}{{{v_{\rm r}}}_{{\rm p}} (R_{*})}\cdot$$ (63)

Here we only newly introduced the relative abundance z_* which accounts for different chemical composition in stellar atmospheres.

The boundary value of ionic density is determined numerically to obtain CAK-type solution (see Sect. 3.9). Boundary electron density is calculated from the condition of quasi-neutrality

 

$${n}_{{\rm e}}=\sum_{a\neq{\rm e}}z_a n_a.$$ (64)

3.8.4 Boundary conditions for electric field

Because we have no critical point condition to determine the intensity of the electric field at the stellar surface, we used the condition of neutrality, which simply sets the gradient of the electric field at the stellar surface to zero (cf. Eq. (48)).

   
3.9 Numerical method

We apply the Henyey method (Henyey et al. 1964), which is a modification of the well-known Newton-Raphson method to solve equations described here together with the appropriate boundary conditions. We use essentially the same method as KKI, except that the vector of variables at each grid point d has the form of

$$\vec{\psi}_d =\left(\rho_{a,d},v_{{\rm r}a,d},T_{a,d},z_{a,d}, E_{\rm d},\Delta v_{{\rm r},d}\right), \quad a={\rm e},{\rm i},{\rm p}$$ (65)

where the velocity difference

$$\Delta {v_{\rm r}}={{v_{\rm r}}}_{{\rm i}} -{{v_{\rm r}}}_{{\rm p}}$$ (66)

may be added to the set of variables to assure better convergence of the models.

First of all we search for the boundary density $\rho_{0}={\rho}_{{\rm i}}(R_{*})$. We compute several wind models (each of them is a result of several Newton-Raphson iterative steps) for the region near the star for different values of $\rho_{0}$ (for more details see KK0, KKI). We select such value of $\rho_{0}$ which allows wind model to pass smoothly through the point defined by the Eq. (57) and to obtain CAK-type solution. After the appropriate value of $\rho_{0}$ is chosen, we compute the wind model downstream of the point defined by the Eq. (57) again using several Newton-Raphson iterative steps.

The detailed method of calculation of Gayley-Owocki heating/cooling term is given in Appendix A.