A&A 417, 1003-1016 (2004)
DOI: 10.1051/0004-6361:20034030

NLTE models of line-driven stellar winds

I. Method of calculation and first results for O stars

J. Krticka^1,2,3 - J. Kubát³

1 - Ústav teoretické fyziky a astrofyziky PrF MU, 611 37 Brno, Czech Republic
2 - Department of Physics and Astronomy, University of Glasgow, Glasgow G12 8QQ, UK
3 - Astronomický ústav, Akademie ved Ceské republiky, 251 65 Ondrejov, Czech Republic

Received 1 July 2003 / Accepted 24 November 2003

Abstract
New numerical models of line-driven stellar winds of late O stars are presented. Statistical equilibrium (NLTE) equations of the most abundant elements are solved. Properly obtained occupation numbers are used to calculate consistent radiative force and radiative heating terms. Wind density, velocity and temperature are calculated as a solution of model hydrodynamical equations. Contrary to other published models we account for a multicomponent wind nature and do not simplify the calculation of the radiative force (e.g. using force multipliers). We discuss the convergence behaviour of our models. The ability of our models to predict correct values of mass-loss rates and terminal velocities of selected late O stars (mainly giants and supergiants) is demonstrated. The systematic difference between predicted and observed terminal velocities reported in the literature has been removed. Moreover, we found good agreement between the theoretical wind momentum-luminosity relationship and the observed one for Cyg OB2 supergiants.

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

1 Introduction

NLTE line-driven stellar wind models have been calculated by various authors (cf. de Koter et al. 1993; Hillier & Miller 1998; Vink et al. 2000; Pauldrach et al. 2001; Gräfener et al. 2002). These sophisticated models cover a wide range of OB and WR star parameters, and it may seem that there is no need for any additional NLTE wind models. However, we believe this is not the case. Firstly, these models employ different simplifying assumptions which cannot be necessarily fulfilled in all model stars. For example, wind velocity is often calculated using either a simplified $\beta$-velocity law (Vink et al. 2001; Gräfener et al. 2002) or using simplified line-force multipliers (Pauldrach et al. 2001; Kudritzki 2002) of Castor et al. (1975) and Puls et al. (2000). Some of these models (e.g., Vink et al. 2001) do not allow independent calculation of observed terminal velocities. Moreover, temperature stratification is often neglected (cf. Kudritzki 2002) assuming that the radiative force is not significantly influenced by temperature (see Pauldrach 1987) or calculated using a simplified expression (Pauldrach et al. 1994; Vink et al. 1999). Although the calculation of correct temperature stratification is included in modern wind codes (cf. Pauldrach et al. 2001), the latest published grid of NLTE calculations of wind temperature dates back to the paper of Drew (1989), who used a relatively constrained set of statistical equilibrium equations. Despite the abovementioned enormous effort in wind model calculations, many unsolved problems remain which may be connected either with the simplifying assumptions mentioned above or with some others, e.g. neglect of wind non-stationarity (cf. Runacres & Owocki 2002), wind 3D structure (Dessart & Owocki 2003), or wind magnetic fields (ud-Doula & Owocki 2002). Among these unsolved problems is the long-standing momentum problem of WR stars (for the latest progress see Gräfener et al. 2002) or problems with fitting of X-ray line profiles (Kramer et al. 2003).

Many of the simplifying assumptions employed in modern wind codes are probably acceptable for high density radiatively-driven stellar winds of O stars, but they may break down for low density stellar winds. Whereas the bulk wind material consists of hydrogen and helium, line-driven stellar winds are accelerated mainly by the absorption of radiation in the resonance lines of trace elements like C, N, O, or Fe. For low-density stellar winds this multicomponent nature may be important for the wind structure itself, influencing both temperature stratification (Springmann & Pauldrach 1992; Krticka & Kubát 2001, hereafter KKII) and wind dynamics (Babel 1995; Porter & Skouza 1999; KKII; Owocki & Puls 2002; Krticka & Kubát 2002). Clearly, such a complicated situation cannot be adequately described by models with a temperature-independent radiative force, or with a radiative force artificially depending on wind temperature via the thermal velocity as has been done, e.g., by KKII. Note that the artificial dependence on the thermal velocity may be removed using an alternative description of the radiative force by Gayley (1995).

The situation is especially appealing now, when wind models suitably describing the atmospheres of first generation stars have been calculated (Kudritzki 2002), because some of these models are also influenced by multicomponent effects (Krticka et al. 2003; see also Kudritzki 2002). Moreover, predicted mass-loss rates of low-metallicity winds of some SMC stars are more than one magnitude higher than the observed ones (Bouret et al. 2003). These discrepancies may be connected with problems with source function gradients (Owocki & Puls 1999).

Although satisfactory determination of mass-loss rates and terminal velocities of O stars can be performed from observations in different wavelength regions, the situation in the B star domain is not so clear. Even for massive B stars there is a significant discrepancy between observations in different wavelength regions (Vink et al. 2000). Even worse, for most main-sequence B stars there is no trustworthy determination of either mass-loss rates or terminal velocities, and, moreover, it is not clear which B stars posses a stellar wind (see Babel 1996 for the discussion of this problem). This is a challenge for sophisticated wind models to help supplement missing observational data. There exist several subclasses of B stars whose peculiar behaviour may be connected with their stellar wind. This is epecially the case of Be stars, because the plausible explanation of the Be phenomenon is still missing (excluding some binaries, for which a binary hypothesis of Kríz & Harmanec (1975) offers a trustworthy explanation).

Last but not least, there is an ambitious idea of the wind momentum-luminosity relationship (see Kudritzki & Puls 2000, and references therein) which may be used as an independent reliable method for determination of stellar distances. Thus, several types of detailed wind models are necessary to verify all of the assumptions of line-driven wind models.

Motivated by these considerations we offer a new independent type of NLTE line-driven wind model to obtain correct values of radiative force, mass-loss rates and terminal velocities in the case of low-density stellar winds.

2 Model assumptions

Hydrodynamic model equations solved in this paper are essentially the same as in KKII. An interested reader will find more details about the solution of multicomponent hydrodynamical equations in KKII and Krticka (2003). Models presented here differ in the method of calculation of the radiative force and the radiative heating/cooling term, and these new models assume slightly different boundary velocities. In the following sections we summarize our model equations. The basic model assumptions are the following:

• we assume a stationary spherically-symmetric flow;
• we solve the continuity equation, the momentum equation and the energy equation for each component of the flow, namely for absorbing ions, nonabsorbing ions (hydrogen and helium), and electrons. Nonabsorbing ions are ions for which the effect of a radiative force may be neglected (e.g., hydrogen and helium if Population I/II stars are considered);
• level occupation numbers of model ions are obtained from NLTE rate equations with an inclusion of a superlevel concept (Anderson 1989);
• radiation transfer is split into two parts. Line radiative transfer is solved in a Sobolev approximation (Rybicki & Hummer 1978) neglecting continuum opacity sources, line overlaps, and multiple scattering. Continuum radiative transfer is formally solved using the Feautrier method with an inclusion of all free-free and bound-free transitions of model ions, however neglecting line transitions. Thus, we neglect line-blocking;
• for the solution of NLTE rate equations together with the radiative transfer equation we apply the Newton-Raphson method for the radiative transfer in lines and lambda iterations in continuum. The simplification of a lambda iteration may, of course, fail for optically thick continua, like the He II continuum (however, see Sects. 2.3 and 3 for a discussion of this issue);
• derived occupation numbers are used to calculate radiative force in the hydrodynamic equations;
• the radiative cooling/heating term is obtained by adding the contributions from all free-free, bound-free and bound-bound transitions of model ions. For the calculation of this term we use the thermal balance of electrons method (Kubát et al. 1999);
• the description of a transition region between a quasi-static stellar atmosphere and supersonic stellar wind is simplified. This region is important for the continuum formation (Pauldrach 1987). This subsonic region was not included in our previous models since we did not solve the radiative transfer equation (see KKII). To obtain correct ionizing flux we included this subsonic region in our models, i.e. the boundary velocity at the stellar surface is set to some low value well below the sonic speed. The particular value of the base velocity is selected in such a way that lowering of this value by a factor of 10 does not change the model structure (i.e. mass-loss rate and terminal velocity) by more than 10%. For most of the wind models this means that the boundary velocity is of the order of 10^-3 of the hydrogen thermal speed;
• the calculation of models below and above the critical point is split in two parts, thus we neglect the influence of the radiation emitted in the region above the critical point on the region below the critical point. However, the influence of the subcritical parts on the supercritical ones is taken into account correctly;
• we neglect Gayley-Owocki heating (Gayley & Owocki 1994) which is not important for stars considered in this paper (see KKII);
• we also neglect the X-ray radiation which may influence the ionization balance.

These assumptions are discussed in more detail in the Sect. 5.

2.1 Equations of statistical equilibrium

Our treatment of statistical equilibrium equations is very similar to that of Pauldrach (1987). We solve rate equations in the form of (Mihalas 1978)

$$% \sum_{j\neq i}N_j P_{ji}-N_i \sum_{j\neq i}P_{ij}=0,$$ (1)

where N_i is the number density of ions in state i relative to the total number density of given atom, and P_ij, P_ji are rates of all processes that transfer an atom from a given state i to state j and back, respectively. They are sums of radiative and collisional rates (both excitation/de-excitation and ionization/recombination). The system of rate Eq. (1) is closed by the equation of normalisation of relative number densities for each model atom (i.e. abundance equation)

$$% \sum_{i} N_i=1.$$ (2)

This equation can be used instead of any equation from the set (1). However, it is computationally convenient to replace the rate equation for the level with the highest relative number density.

We use the relative number density N_i instead of the commonly used number density to be consistent with the hydrodynamical part of the code, because relative quantities are advantageous from the numerical point of view. Note that relative number densities were used for the solution of the equations of statistical equilibrium by Lucy (2001).

   
2.1.1 Model atoms and rates

Radiative excitation rates R_ij are given by

$$% N_i R_{ij} = N_i B_{ij} \bar J_{ij},$$ (3)

where $\bar J_{ij} = \int_0^\infty \phi_{ij}(\nu) J_\nu {\rm d}\nu$ is profile-weighed line intensity and

$$% B_{ij} = \frac{4\pi^2 {e}^2}{h\nu_{ij}m_{\rm e}c}f_{ij},$$ (4)

where $\nu_{ij}$ is line frequency, f_ij is corresponding oscillator strength, e is elementary charge and $m_{\rm e}$ is electron mass.

Similarly, the radiative de-excitation rates R_ji are

$$% N_j R_{ji} = N_j \left(A_{ji}+B_{ji}\bar J_{ij}\right),$$ (5)

where the remaining Einstein coefficients are given by:

$$% A_{ji}=\frac{2h\nu_{ij}^3}{c^2}B_{ji}, \qquad B_{ji}=\frac{g_i}{g_j} B_{ij},$$ (6)

and g_i, g_j are statistical weights of involved levels.

Radiative ionization and recombination rates are given by:

 

$$4\pi N_i \int_{\nu_i}^{\infty} \frac{\alpha_{i,\nu}}{h\nu}J_\nu~{\rm d}\nu,$$                                     N_i R_ik = (7a)

$$4\pi N_k \left(\frac{N_i}{N_k}\right)^\ast \int_{\nu_i}^{\infty} ... ...frac{2h\nu^3}{c^2}+J_\nu\right] {\rm e}^{-\frac{h\nu}{kT_{\rm e}}} ~{\rm d}\nu,$$ N_k R_ki = (7b)

where $\alpha_{i,\nu}$ is the corresponding photoionization cross-section, $J_\nu$ is the mean continuum intensity (line absorption spectrum is neglected in calculations of bound-free rates), $\nu_i$ is the ionization frequency of level i, $T_{\rm e}$ is the electron temperature, and an asterisk denotes LTE values.

Upward (C_ij) and downward (C_ji) collisional rates (for both excitation and ionization) are calculated using the expressions

 

$$N_i C_{ij} = N_i n_{\rm e} \Omega_{ij} (T_{\rm e}),$$ (8a)

$$N_j C_{ji} = N_j n_{\rm e} \left(\frac{N_i}{N_j}\right)^* \Omega_{ij} (T_{\rm e}),$$ (8b)

where $\Omega_{ij} (T_{\rm e})$ is the integrated collisional cross-section (collision strength) and $n_{\rm e}$ is the electron number density. Finally, rate coefficients P_ij for both excitation and ionization are sums of radiative and collisional terms

 

P_ij=R_ij+C_ij. (9)

Table 1: Atoms and their ionization stages included in the NLTE calculations. In this table, level means either an individual level or a set of levels merged into a superlevel.

We solve NLTE rate Eq. (1) for each of the elements listed in Table 1 separately. This allows us (in addition to a fixed temperature) to fix the number of free electrons during the solution of the equations of statistical equilibrium. This step greatly reduces the computational needs, since instead of one large rate matrix we have to deal with a moderate number of smaller matrices for the elements considered.

The model atoms are predominantly taken from TLUSTY models (Hubeny 1988; Hubeny & Lanz 1992, 1995; Lanz & Hubeny 2003). The set of model atoms is slightly modified for the case of hot star wind modelling. For this purpose we used data from the Opacity Project (Seaton 1987; Luo & Pradhan 1989; Sawey & Berrington 1992; Seaton et al. 1992; Butler et al. 1993; Nahar & Pradhan 1993) and from the Iron Project (Hummer et al. 1993; Bautista 1996; Nahar & Pradhan 1996; Zhang 1996; Bautista & Pradhan 1997; Zhang & Pradhan 1997; Chen & Pradhan 1999).

We use a similar strategy for the preparation of model atoms that enter the calculations as was already applied in TLUSTY. With the exception of  Fe and  Ni, levels with low energies (relative to the ground level of a given ionization stage) are treated individually (neglecting the hyperfine splitting), whereas higher levels with similar properties are merged into superlevels. For  Fe and  Ni all levels (except the ground ones) are merged into superlevels. More details on the construction of model atoms and atomic data used are given in Appendix A.

Procedures for the calculation of $\Omega_{ij} (T_{\rm e})$ are taken from the TLUSTY code. The partition function approximations are mostly from Smith & Dworetsky (1988) and Smith (1992). Moreover, the partition function approximations taken from the TLUSTY code are applied for  C, N, O, Fe VI, and Fe VII.

2.1.2 Superlevels

Compared to the pioneering work of Anderson (1989), we use a simplified merging of individual levels into superlevels. Given the sublevels l with energies E_il, relative number densities N_il, and statistical weights g_il, the total population of superlevel is given by

$$% N_i=\sum_l N_{il},$$ (10)

and the energy of superlevel is

$$% E_i=\frac{\sum_l g_{il} E_{il}}{\sum_l g_{il}}\cdot$$ (11)

Thus, we effectively neglect the $\exp(-E_{il}/kT_{\rm e})$ term. Similarly, the photoionization cross section of superlevel  $\alpha_{i,\nu}$ is calculated by

$$% \alpha_{i,\nu}=\frac{\sum_l g_{il} \alpha_{il,\nu}} {\sum_l g_{il}},$$ (12)

and the relative number density of a sublevel of a given superlevel is approximated as

$$% N_{il}=\frac{g_{il}N_i}{\sum_l g_{il}}\cdot$$ (13)

Numerical tests showed that the neglect of the $\exp(-E_{il}/kT_{\rm e})$ term does not have a significant effect on the calculated radiative force.

2.2 Radiative transfer

Radiative flux at the base of the stellar wind $H_{\rm c}$ is taken from H- He spherically symmetric NLTE model atmospheres of Kubát (2003). The line transfer equation and continuum transfer equation are solved separately.

2.2.1 Line radiative transfer

The averaged line intensity is calculated using the Sobolev approximation (Rybicki & Hummer 1978)

$$% \bar J_{ij}=(1-\beta)S_{ij} + \beta_{\rm c} I_{\rm c},$$ (14)

where escape probability functions $\beta$ and $\beta_{\rm c}$ are given by integrals

  

$$\beta_{\rm c} =\frac{1}{2}\int_{\mu_*}^{1}{\rm d}\mu\frac{1-{\rm e}^{-\tau_\mu}} {\tau_\mu},$$ (15)

$$\beta =\frac{1}{2}\int_{-1}^{1}{\rm d}\mu\frac{1-{\rm e}^{-\tau_\mu}}{\tau_\mu},$$ (16)

$\mu_*=\left(1-R_*^2/r^2\right)^{1/2}$, and $I_{\rm c}=4H_{\rm c}$. The Sobolev optical depth $\tau_\mu$ is given by (Castor 1974; Rybicki & Hummer 1978)

$$% \tau_\mu=\frac{\pi e^2}{m_{\rm e}\nu_{ij}} \left(\frac{n_i}... ...f_{ij} \frac{r}{{{v_r}}_{{\rm i}} \left(1+\sigma\mu^2\right)},$$ (17)

where ${{v_r}}_{{\rm i}}$ is the radial velocity of absorbing ions, n_i, n_j are the number densities of individual states,

$$% n_i=N_iz_{\rm atom}n_{\rm H}$$ (18)

(similarly for n_j), where $z_{\rm atom}=n_{\rm atom}/n_{\rm H}$ is the number density of a given atom relative to the hydrogen number density $n_{\rm H}$,

$$% n_{\rm H}= \frac{\rho_{\rm p}}{m_{\rm H}+z_{{\rm He}}m_{{\rm He}}},$$ (19)

where $\rho_{\rm p}$ is the density of a passive component (hydrogen and helium). The variable $\sigma$ was introduced by Castor (1974),

$$% \sigma=\frac{{\rm d}\ln {{v_r}}_{{\rm i}} }{{\rm d}\ln r} -1.$$ (20)

A suitable method for the calculation of escape probability functions $\beta$ and  $\beta_{\rm c}$ was described by KKII. Finally, the line source function is

$$% S_{ij}=\frac{N_j A_{ji}}{N_i B_{ij}-N_j B_{ji}}\cdot$$ (21)

2.2.2 Continuum radiative transfer

For the calculation of the continuum intensity we use a procedure for the solution of the transfer equation based on the Feautrier method in the spherical coordinates (for a description, see Mihalas & Hummer 1974 or Kubát 1993). However, for the solution of the obtained tridiagonal system of equations we use the numerical package LAPACK (http://www.cs.colorado.edu/~lapack, Anderson et al. 1999) instead of the classical Gaussian scheme. All bound-free transitions of explicit ions, free-free transitions of hydrogen and helium and light scattering due to free electrons contribute to the source function. The radiative transfer equation is solved only above the hydrogen ionization edge. For frequencies lower than the hydrogen ionization limit we use an optically thin approximation $J_\nu(r) = 4 W(r) H_{\rm c}$, where $W(r)=\frac{1}{2}\left(1-\mu_*\right)$ is the dilution factor.

   
2.3 Solution of NLTE equations

For the solution of the NLTE rate Eqs. (1) and (2) together with the radiative transfer equation we use a combination of lambda iterations and of a complete linearization. For the continuum transfer we use simple lambda iterations. Although the stellar wind is optically thick near to the star below the He II ionization limit, these iterations converge relatively well (see also Sect. 3). In lines, we use complete linearization of the radiative transfer equation in a simple Sobolev approximation (neglecting continuum radiation, see Eq. (14)). The line intensity  $\bar J_{ij}$ depends only on the level populations of the given atom for this case. Thus, the statistical equilibrium equations for different atoms are independent (for a fixed hydrodynamical structure). Hence, the following scheme may be obtained from statistical equilibrium Eqs. (1) and (2) for the calculation of the correction  $\delta N_i$ of relative number densities in the given iteration step

 

$$\sum_{j\neq i}\left[N_j\frac{\partial P_{ji}}{\partial N_i}\delta... ...\right)\delta N_i+ N_i\frac{\partial P_{ij}}{\partial N_j} \delta N_j\right]=0,$$ (22)

where the terms ${\partial P_{ij}}/{\partial N_j}$ are calculated using Eqs. (14)-(21). For details see Appendix B. NLTE Eq. (22) are solved for each element separately. Note that Eq. (22) is an application of a Newton-Raphson method to statistical equilibrium Eq. (1). For the solution of the system of NLTE Eq. (22) we use the numerical package LAPACK. Note that we can again take full advantage of LAPACK subroutines designed for general band matrices because our NLTE rate equation matrices have a band structure.

2.4 Derivatives of occupation numbers

The hydrodynamic equations are solved using a Newton-Raphson method. To obtain a well converged model we have to know the explicit values of derivatives of the radiative force and of the radiative heating/cooling term with respect to fluid variables (i.e. densities, temperatures, etc.). However, because the radiative force and the radiative heating/cooling term depend on level populations, we have to calculate derivatives of level populations with respect to fluid variables. These can be calculated in the following way. For clarity, we rewrite NLTE rate Eq. (1) in a symbolic form

 

$$% \sum_{j\neq i}P_{ji} (\vec{N},J_\nu(\vec{N},{\vec{h}}),\vec{h}) N_j -N_i\sum_{j\neq i} P_{ij} (\vec{N},J_\nu(\vec{N},\vec{h}),\vec{h}) =0,$$ (23)

where we denoted all explicit dependencies of rate equations. The elements of vector $\vec{N}$ are level populations N_i, the elements of vector $\vec{h}$ are fluid variables $\vec{h}=(n_{\rm e}, \rho_{\rm p},T_{\rm e},{{v_r}}_{{\rm i}} )$, and the mean continuum intensity $J_\nu$ depend on $\vec{h}$ and $\vec{N}$ via the radiative transfer equation. Clearly, it is not feasible to obtain correct explicit values of $\partial N_i / \partial h$ (where h denotes elements of vector $\vec{h}$), since the level populations depend on the solution of radiative transfer equation in continuum. However, it is possible to obtain approximate values of these derivatives neglecting the explicit dependence of rates P_ij on the continuum intensity. Deriving Eq. (23) with respect to h and neglecting the dependence on $J_\nu(\vec{N},\vec{h})$ we obtain a system of equations for $\partial N_i / \partial h$

 

$$% \sum_{j\neq i}\left[\frac{\partial P_{ji}}{\partial h} N_j+ \fr... ...j}}{\partial N_i} N_i+P_{ij}\right) \frac{\partial N_i}{\partial h} \right] =0.$$ (24)

All calculated derivatives $\partial N_i / \partial h$ are inserted into the Newton-Raphson matrix of hydrodynamic equations (corresponding to the radiative force, radiative heating/cooling term and ionization charge, see KKII and Krticka 2001).

Note that the system of equations for derivatives (24) has the same form as the system of NLTE rate Eq. (22), however with a different right-hand side. Thus, because the LAPACK package uses an LU decomposition for the solution of the system of linear equations, the values of $\partial N_i / \partial h$ can be calculated using data from the previous solution of NLTE Eq. (22).

This method is especially useful in the case when the intensity of radiation in the continuum does not depend on local properties (i.e., for an optically thin continuum, when it is given by the intensity emerging from the underlying photosphere). From the mathematical point of view, this method splits the Jacobi matrix into several parts and it is equivalent to a linearization of the whole system of model equations (i.e. hydrodynamics, NLTE equations and radiative transfer). However, from the numerical point of view, it saves a huge amount of computer time because it is not necessary to invert large matrices. A similar method has been used, e.g., by Anderson (1985), and Kubát (1994,1997) referred to it as to the implicit linearization of the b-factors.

2.5 Calculation of the radiative force

The radiative force is calculated using the Sobolev approximation directly as the sum of contributions of individual lines after Castor (1974),

$$% {g}_{{\rm i}}^{{\rm rad}}=\frac{8\pi}{{\rho}_{{\rm i}} c^2}... ...gma\mu^2\right) \left(1-{\rm e}^{-\tau_\mu}\right){\rm d}\mu,$$ (25)

where $\tau_\mu$ is given by Eq. (17). Occupation numbers N_i calculated during the NLTE iteration step are used for a calculation of the radiative force. Note that in Eq. (17) number densities of individual levels are used instead of the total number density of a superlevel. These number densities may be calculated from Eq. (13) for those levels which are explicitly included in the NLTE rate equations. However, non-explicit levels, i.e., some of the higher excited levels or excited levels of the highest ion of an atom, are also included in the calculation of the radiative force. For these non-explicit levels we use either the Boltzmann equilibrium at the wind temperature in the case of metastable levels or the radiative excitation equilibrium of Abbott & Lucy (1985) for other levels. The influence of these non-explicit levels on the radiative force is very small, however.

The term $\partial {g}_{{\rm i}}^{{\rm rad}} / \partial\left({\rm d}{{v_r}}_{{\rm i}} /{\rm d}r\right)$ is calculated explicitly after (25) where the variable $\sigma$ depends on a velocity gradient via (20). This term is included in the critical point condition (generalised CAK condition, KKII, Eq. (57) therein). Derivatives of number densities calculated by means of Eq. (24) are used to obtain derivatives of the radiative force with respect to fluid variables (densities, velocities, and temperatures) for the Newton-Raphson iterations.

Oscillator strengths necessary for the calculation of the radiative force are extracted from the VALD database (Piskunov et al. 1995; Kupka et al. 1999). This set of data consists of about 180 000 line transitions in the wavelength interval 250-10 000 Å.

2.6 Radiative cooling and heating

The radiative cooling and/or heating term is calculated using the thermal balance of electrons (Kubát et al. 1999). For the calculation of this term we include all the explicit radiative bound-free and collisional (both bound-bound and bound-free) transitions and free-free transitions of hydrogen and helium. Again, in order to obtain a well-converged model, derivatives of number densities calculated by means of Eq. (24) are used to obtain derivatives of the radiative cooling and/or heating terms with respect to fluid variables. The calculated derivatives of the radiative cooling/heating term are inserted into the Jacobi matrix of the Newton-Raphson iteration of the model hydrodynamic structure. Detailed expressions for the derivatives of heating and cooling rates are presented in Appendix C.

2.7 Electron density

The charge of the passive component is calculated from ionization fractions of both hydrogen and helium. This enables us to calculate the correct value of the electron density, which is important for the proper solution of NLTE rate equations.

   
2.8 Iteration scheme for the model calculation

	Figure 1: Iteration scheme for the model calculation.
Open with DEXTER

The initial estimate of fluid variables for the Newton-Raphson iterations (namely the so-called beta-law for velocities and density obtained from the continuity equation) was described by KKII. For the initial value of the intensity we use an optically thin approximation, i.e., the radiation field is defined by the lower boundary condition at the beginning of our calculation.

The computational procedure involves two similar basic steps, namely solution below and above the wind critical point. This splitting into two steps is necessary and is related to the critical point condition. In the first step the value of the base density, which allows the model to smoothly pass the critical point, is obtained. Afterwards, the model above this critical point is calculated.

The calculation of a given model below the critical point consists of three nested iteration cycles (see Fig. 1). In the innermost iteration cycle the NLTE rate Eq. (1) together with the radiative transfer equation both in lines (14) and in continuum are solved. For this innermost iteration cycle the fluid variables (i.e. densities, velocities and temperatures of individual wind components) are kept fixed.

In the middle cycle we iterate the hydrodynamical structure of the model for fixed level populations. This iteration cycle is essentially the same as the iteration cycle of three-component nonisothermal wind models described by KKII. To save computing time, for large relative changes of fluid variables (greater than 10^-2) we perform only 3 iterations of NLTE rate equations, but when the relative changes of fluid variables are lower, and we may expect that we are closer to the correct solution, we perform more iterations of NLTE rate equations (up to 50). We assume that model is well converged if the relative changes of both fluid variables and occupation numbers are lower than 10^-5.

In the outermost iteration cycle we solve for the value of the base density for which a model smoothly passes through the critical point. Again, these iterations were already described by KKII. For the base density (which determines the value of mass-loss rate ${\rm d}{\mbox{\eufont M}}/{\rm d}t$) an accuracy of 1% is sufficient.

After the base wind density is known, we calculate wind model above the critical point by a sequence of iterations of fluid variables. Each of these iterations is preceded by the set of NLTE iterations. Boundary conditions of the model above the critical point are specified at the critical point and are taken from the model below the critical point. The model below the critical point is fixed during the solution above the critical point.

Splitting of a wind solution into two parts may cause problems since the radiation emitted from the outer region (i.e. above the CAK point) may influence the inner region. However, test calculations where we changed the outer boundary of inner solution showed that this is not the case, i.e. that the inner wind solution depends only marginally on the outer solution. The reason for such a behaviour is probably the strong dependence of the wind continuum optical depth on frequency. Calculations show that the stellar wind above the CAK point is essentially optically thin for continuum radiation with frequencies below the He II ionization limit, thus radiation is not absorbed in the continuum above the CAK point and does not influence the stellar wind below the CAK point. On the other hand, stellar wind above the CAK point is optically very thick for radiation with frequencies above the He II ionization limit. Again, due to this large optical thickness the radiation emitted in the regions above the CAK point does not influence the inner stellar wind. The influence of lines above the critical point on the subcritical region has been also neglected.

   
3 Convergence behaviour of models

Figure 2: The convergence of the models below ( upper panels) and above ( lower panels) the critical point. Upper left panel: the convergence of level populations below the critical point. Note that NLTE iterations in consecutive hydrodynamical iterations are combined in this plot. This explains the relatively large number of iterations and peaks that occur after the hydrodynamical iteration step. The oscillating behaviour for the relative change of order 10-5is due to the weakly populated levels, because there are no oscillations in the plot of convergence of relative change of major levels (i.e. levels with Ni>0.01). Upper right panel: the convergence of hydrodynamical structure (maximum relative change of all hydrodynamical variables, i.e. densities, velocities, temperatures and charges of individual wind components) of the model below the critical point. The slightly oscillating behaviour around the relative change 10-3is caused by the change of maximum allowed number of NLTE iteration steps (see the text). Lower left panel: the convergence of level populations of the model above the critical point. Again, NLTE iterations in different hydrodynamical iterations are put together and peaks occur after the hydrodynamical iteration step. Note the fast convergence of NLTE population between successive hydrodynamical iterations. Lower right panel: the convergence of the hydrodynamical structure of the model above the critical point.

Open with DEXTER

As an example we study the convergence behaviour of a model corresponding to $\alpha$ Cam (HD 30614) with parameters given in Table 2. In Fig. 2 we plotted relative change of population numbers, and hydrodynamical variables (i.e. densities, velocities, temperatures and charges of individual wind components) for the model below the critical point and for the model above the critical point.

First we discuss the convergence behaviour of the model below the critical point. For the plot of relative change of level populations (upper left panel of Fig. 2) we combine NLTE iterations performed during individual hydrodynamical iteration steps for fixed base wind density (corresponding to the final CAK solution). Thus, in this plot the peaks corresponding to the hydrodynamical iteration steps occur. They are caused by a reaction of NLTE populations to a significant change of the fluid variables after the calculation of the new estimate of the hydrodynamic structure. For the relative change of about 10^-5 relative changes oscillate. However, if we plot only relative changes of populations of levels with largest occupation numbers (levels with N_i>0.01), no oscillations appear. Thus, we conclude that the oscillatory behaviour is caused by weakly populated levels and that their influence on the convergence behaviour of the hydrodynamical structure is negligible (upper right panel of Fig. 2). This oscillatory behaviour is not a typical feature of convergence of our models.

The plot of hydrodynamical structure convergence (upper right panel of Fig. 2) shows very fast convergence which is slightly slowed down by the oscillatory behaviour around the relative change 10^-3. This is caused by the change of the maximum allowed number of NLTE iteration steps (see Sect. 2.8). Although due to this change the convergence is slightly slowed down, the total computing time is less because we do not lose time with calculation of precise occupation numbers for hydrodynamical structure which is far from the converged one. Finally, note that the relatively small initial relative change of fluid variables (of order 10^-2) of the hydrodynamical structure is caused by the way the critical point is calculated, i.e. that models with similar base densities are calculated. Thus, the initial change of fluid variables of the final model below the critical point is very small. The convergence behaviour of previous models (i.e. models which do not correspond to the critical solution) is similar, however the initial change of fluid variables is higher.

The convergence behaviour of the model above the critical point (lower panels of Fig. 2) is similar to that below the critical point. The convergence of occupation numbers (for a given hydrodynamical structure) is faster because the stellar wind is optically thinner (note that this fast convergence of NLTE iterations combined with larger number of iterations of hydrodynamic structure yields something like a "forest'' of iterations in the plot, which is caused by a large change of occupation numbers after a new hydrodynamic structure was calculated. The convergence of a hydrodynamical structure (lower right panel of Fig. 2) is however slower than below the critical point. This is partially caused by lower maximum change of wind parameters allowed during Newton-Raphson iteration step. This constraint of the maximum change of wind parameters is necessary for the stability of the solution (see Krticka 2003).

The convergence of the model temperature is very fast in the case of a model below the critical point. The relative change of temperature is typically ten times lower than the relative change of the rest of the hydrodynamical variables in this case. On the other hand, the convergence behaviour of a model above the critical point is dominated by the temperatures since the relative changes of temperature are the highest.

4 Comparison of calculated and observed quantities

Table 2: Stellar and wind parameters of selected O stars from our test sample. Stellar parameters (columns  ${\mbox{\eufont M}}$, R_*, $\mbox{$T_{{\rm eff}}$ }$) were adopted from Lamers et al. (1995). Observed mass loss rates were taken from [1] - Leitherer (1988); [2] - Scuderi et al. (1992); [3] - Lamers & Leitherer (1993); [4] - Puls et al. (1996); [5] - Crowther & Bohannan (1997); [6] - Scuderi et al. (1998); [7] - Lamers et al. (1999); [8] - BG (see the text), predicted values of  $\dot{\mbox{\eufont M}}$ and  $v_{\infty }$ are determined using our code. Note that terminal velocities calculated by KKII were obtained using simplified force multipliers after Abbott (1982).

Because in future we intend to calculate mainly wind models of B stars, to test the reliability of our assumptions we calculated NLTE wind models of stars with spectral type relatively close to B, i.e. later O stars. For this purpose we selected O stars with a spectral type O6 or later from the set of Lamers et al. (1995, hereafter LSL). We included only stars for which reliabile determination of their mass-loss rate exists in the literature. Stellar parameters (masses, radii and effective temperatures) taken from LSL are listed in Table 2.

4.1 Terminal velocities

Wind terminal velocities of galactic O stars can be measured with a typical uncertainty of $100~{\rm km}~{\rm s}^{-1}$. Thus, in principle, they may provide a relatively strict test for any wind model. Unfortunately, because terminal velocities depend on surface escape velocities (see Castor et al. 1975) the precise prediction of wind terminal velocity relies on the precise knowledge of stellar mass and radius. Consequently, the long-standing problem of "mass discrepancy'' in O stars (cf. Herrero et al. 1992; Lanz et al. 1996) complicates a straightforward comparison of calculated and observed terminal velocities. This is especially the case of stars HD 30614, HD 188209, and HD 218915, for which we used the same stellar parameters, but which have different observed terminal velocities.

The comparison of observed (taken from LSL and Bianchi & Garcia 2002, hereafter BG) and calculated terminal velocities for model stars shown in Fig. 3 implies that our NLTE wind models are able to predict relatively correct values of wind terminal velocities. The fit is improved compared to our previous calculations presented in KKII (see Table 2). Note that our calculations consistently solved the problem of high theoretical terminal velocities obtained by LSL using the so called "cooking formula'' of Kudritzki et al. (1989), see also Table 2 and Fig. 12 in KKII. However, since our models include several simplifying assumptions, it is necessary to test this conclusion against more advanced models.

It remains unclear why LSL obtained too high values of terminal velocities when theoretical calculations (KKII) with the same force multipliers yielded terminal velocities consistent with observations. Note that Pauldrach et al. (2001) also predicted correct terminal velocities, however for a different set of stars.

	Figure 3: Comparison of calculated and observed (taken from LSL, BG) terminal velocities.
Open with DEXTER

	Figure 4: Comparison of calculated (crosses) and observed (taken from LSL, BG, errorbars) ratio of the terminal velocity to the escape velocity $v_\infty /v_{{\rm esc}}$.
Open with DEXTER

The calculated terminal velocities are, in fact, not too different from those calculated by KKII, although the terminal velocities presented by KKII were obtained simply using the force multipliers of Abbott (1982). This relatively surprising conclusion can be attributed to the fact that the $\alpha$ parameter calculated properly using line statistics (Puls et al. 2000) is very close to the values calculated by Abbott (1982). And because the terminal velocities depend on the $\alpha$ value, terminal velocities calculated in this paper and those using Abbott's multipliers (KKII) do not differ significantly. However, properly calculated k parameters (Puls et al. 2000) are significantly lower than those calculated by Abbott (1982). Indeed, our mass-loss rates (which are, in fact, influenced by the k parameter) are systematically lower than those of KKII and are close to the observed values (see also the next section for the discussion of mass-loss rates).

The stellar wind theory predicts a strong correlation between wind terminal velocity and stellar escape velocity. Thus, we have also compared the ratio of the observed and calculated terminal velocity to the escape velocity $v_\infty /v_{{\rm esc}}$ (see Fig. 4, $v_{{\rm esc}}$ is calculated using stellar parameters from Table 2). Our derived mean ratio $v_\infty/v_{{\rm esc}} \sim 2.6$ is in a good agreement with the mean observed $v_\infty /v_{{\rm esc}}$ ratio 2.5 derived by LSL for stars with an effective temperature higher than 25 000 K.

4.2 Mass loss rates

The credibility of observed mass-loss rates in most stars is much worse than their wind terminal velocities. Although in general the wind models are able to predict correct mass-loss rates (cf. Vink et al. 2000), there is a large scatter for individual stars. This situation is demonstrated in Fig. 5 where we compare calculated and observed mass-loss rates. We compiled observed mass-loss rates based on the H$\alpha$ measurement (Leitherer 1988; Scuderi et al. 1992; Puls et al. 1996; Crowther & Bohannan 1997; Scuderi et al. 1998; Lamers et al. 1999) and on the radio emission (Lamers & Leitherer 1993; Scuderi et al. 1998) together with mass-loss determinations obtained using UV lines by BG. Generally, UV line-diagnostic is believed to be badly affected by the unknown wind ionization fractions (Lamers et al. 1999). This is especially the case of the O star winds for which unsaturated wind lines originate only from trace elements. However, precise modelling can overcome these difficulties (BG).

We conclude that there is generally a good agreement between predicted and observed mass-loss rates. However, significant differencies exist for individual stars. Similarly, comparison of mass-loss rates calculated in this paper and by Vink et al. (2000) shows relative good agreement between these predicted mass-loss rates. This possibly supports the conclusion of Puls (1987) that multiple scattering (included into Vink et al. (2000) models) does not significantly influence mass-loss rates of weaker stellar winds.

An important caveat has to be mentioned. The observed mass-loss rates may be influenced by wind clumping (e.g. Crowther et al. 2002; Massa et al. 2003; Bouret et al. 2003) and thus, the observed mass-loss rates shall be lower in such a case (see also Kramer et al. 2003, for the discussion of effect of clumping on the theoretical profiles of X-ray lines). However, in such a case clumping would influence also the theoretical mass-loss rates (cf. Gräfener 2003).

	Figure 5: Comparison of calculated and observed mass-loss rates.
Open with DEXTER

	Figure 6: Comparison of theoretical mass-loss rates.
Open with DEXTER

4.3 Wind momentum luminosity relationship

	Figure 7: Comparison of calculated (crosses) and observed (errorbars) wind momentum (determined by Herrero et al. (2002) for Cyg OB2 association). Linear regression of calculated data is also plotted.
Open with DEXTER

The wind momentum-luminosity relationship (see Kudritzki & Puls 2000) may provide a possibility to obtain precise measurements of stellar and galaxy distances. Thus, the exact calibration of this relationship is necessary. Moreover, the wind momentum depends on the stellar mass only marginally, and thus it provides a more reliable tool for the investigation of the consistency between theoretical models and reality (Puls et al. 1996).

We have compared a theoretical modified wind momentum-luminosity relationship $\dot{\mbox{\eufont M}} v_\infty \left(R/R_{\hbox{$\odot$ }}\right)^{1/2}$ obtained for considered stars with our NLTE wind models with a relationship derived from observations of supergiants by Herrero et al. (2002) for the Cyg OB2 association (see Fig. 7). Again, there is a relatively good agreement between calculated and observed data.

However, recent observational studies (Repolust et al. 2004; Markova et al. 2004) indicate that the situation is probably more complicated. They concluded that there are two different wind momentum-luminosity relationships for supergiants (or, probably more correctly, for stars with the H$\alpha$ line in emission), and for other stars (mainly for giants and main-sequence stars). However, theoretical studies (cf. Vink et al. 2000) do not predict such a difference. Also results from our sample do not exhibit any significant difference between wind momentum-luminosity relationships for supergiant and for non-supergiant stars. Repolust et al. (2004) conclude that the observed difference may be a consequence of clumping and that observed mass-loss rates of high-density winds (mostly supergiant winds) may be artificially enhanced by clumping. To test this idea we calculated our predicted wind momentum-luminosity relationship and compared it with both the theoretical one (Vink et al. 2000) and with the observed wind momentum-luminosity relationship for supergiants and non-supergiant stars (see Table 3). Apparently, our predicted regression to the linear wind momentum-luminosity relationship

$$% \log\left[\dot{\mbox{\eufont M}}v_\infty\left(R/R_{\hbox{$\... ...ght)^{1/2}\right]=x\log(L/L_\odot)+\log D_0, \quad {\rm (CGS)}$$ (26)

where x and D₀ are fit parameters, agrees well with the theoretical results of Vink et al. (2000) and with the observed relationship for non-supergiant stars. However, the observed wind momentum-luminosity relationship for supergiants is different. This confirms the finding of Repolust et al. (2004) that the theoretical wind momentum-luminosity relationship corresponds to the observed one for non-supergiant stars. Note, however, that most supergiants observed by Herrero et al. (2002) represent exceptions to this rule since most of them do not have an H$\alpha$ line in emission.

Table 3: Comparison of wind momentum-luminosity relationship from different sources, both theoretical and observational (see Eq. (26)). Parameters of the sample denoted as "supergiants'' were obtained by Repolust et al. (2004) combining their own results and results by Herrero et al. (2002) for supergiants. Similarly values of the sample "giants and dwarfs'' were obtained.

   
5 Discussion

The main improvement of the presented models is that the radiative force is calculated without simplifying assumptions about the line-distribution function (Puls et al. 2000), i.e., the parameters k, $\alpha$, and $\delta$ (Castor et al. 1975; Abbott 1982) are no longer required. A correctly-calculated wind temperature structure influences occupation numbers and thus, also the radiative force. Moreover, these models are, in principle, able to consistently calculate a multicomponent wind structure, which is not necessary for the case of O star winds, but that will be invaluable for calculations of low-density B star winds.

On the other hand, there are several assumptions involved in the calculation of these models which can affect their validity.

First, we assume a stationary and spherically symmetric stellar wind. Both of these assumptions put severe limitations on our models. Variable phenomena, which are often observed in spectra of O stars, cannot be properly described with a stationary spherically symmetric model of a wind. Spherical symmetry also prevents a proper description of rotating winds. On the other hand, a detailed numerical study of wind stability (cf. Feldmeier et al. 1997; Runacres & Owocki 2002), and of effects of stellar rotation on the stellar wind (cf. Petrenz & Puls 2000) showed that basic wind parameters (i.e. terminal velocity and mass-loss rate) can be well reproduced by stationary and spherically symmetric wind models (for rotational velocities lower than the critical rotational velocity). However, wind clumping may affect the values of theoretical mass-loss rates (Gräfener 2003). Finally, shocks caused by wind instabilities or by the presence of strong magnetic fields (ud-Doula & Owocki 2002) are sources of strong X-ray and EUV emission. This emission influences the ionization balance due to the Auger ionization (Cassinelli & Olson 1979) and direct photoionization (Pauldrach et al. 1994). On the other hand, MacFarlane et al. (1994) showed, that for O stars the X-ray radiation influences the ionization balance of trace elements only. Thus, we conclude that at least wind models of hotter star from our sample are not significantly influenced by X-ray shock emission.

Another sources of possible errors are included in the approximate form of model equations itself. First, multiple scattering and line overlap may influence the mass-loss rate, especially for high-density winds (see Abbott & Lucy 1985; Puls 1987; Vink et al. 1999). However, we plan to apply our code first to models of low-density winds, for which the effect of multiple scattering is not important.

Second, radiative transfer in lines may be influenced by the continuum effects. To overcome this limitation, we plan to include continuum radiation in the formulation of the radiative transfer equation in the Sobolev approximation into our models after Hummer & Rybicki (1985). For some lines, specially for the strong ones, the Sobolev approximation may not be credible. Thus, we plan to test the reliability of the application of the Sobolev approximation for the calculation of the radiative force by using a correct solution of the radiative transfer equation in moving media (Korcáková & Kubát 2003). This will also allow as to account for the self-shadowing by photospheric lines which is important for thin winds (Babel 1996, together with calculation of model atmospheres with inclusion of more elements). Moreover, it will enable a study of instabilities connected with the source-function gradients (Owocki & Puls 1999).

The effect of wind blanketing (Abbott & Hummer 1985) was neglected in our calculations to make the treatment of the wind easier. However, it is clear that this effect causes backwarming of stellar surface and may lead to different values of effective temperatures (e.g. Herrero et al. 2002).

The lambda-iterations employed for the solution of rate equations together with radiative transfer in continuum may be source of another potentially very serious problem. However, the stellar wind is optically thin in most continuum frequencies and the optically thick region occurs mostly relatively near to the star. On the other hand, the radiative transfer in lines in a Sobolev approximation is solved consistently with rate equations. Thus, due to these reasons and due to a relatively good convergence of our models we conclude that we are able to calculate consistent occupation numbers and wind structure even with this approximation.

The processes of ionization and excitation can be very slow compared to the typical flow-time. Thus, the ionization state of the wind can be "frozen'' in a gas parcel and, consequently, can be advected out. Such a case in the stationary spherically symetric stellar wind is described by the rate equation in the form of (Mihalas 1978)

$$% \frac{1}{r^2}\frac{{\rm d}\left(r^2 n_i {{v_r}}_{{\rm i}} \... ...rm d}r} = \sum_{j\neq i}n_j P_{ji}-n_i \sum_{j\neq i}P_{ij}=0,$$ (27)

instead of Eq. (1). However, numerical tests showed that this phenomenon can be neglected for most considered levels (see also Lamers & Morton 1976). Moreover, there would not be any radiatively driven wind if the opposite is true since radiatively driven stellar wind is enabled by high transition rates P_ij (Gayley 1995). On the other hand, at first sight it is not clear whether this assumption is also fulfilled in the case of strong wind-shocks.

Different approximations of NLTE model atoms may influence the wind parameters obtained. Besides the well-known errors in the atomic databases, inclusion of superlevels may badly deteriorate calculated models. Whereas the effect of our simplified inclusion of superlevels is found to be small, neglected collisional transitions between individual levels forming the superlevel may affect wind temperature structure (cf. Hillier & Miller 1998). Moreover, another subtle basic effect may influence the calculated occupation numbers. The calculated ionization fractions depend on the partition function approximation. Especially, in the relatively thin wind environment (where a large number of higher excited levels have high occupation probabilities) different partition function approximations result in different final ionization fractions. Thus, it would be convenient to include more sophisticated methods for the calculation of partition functions (cf. Hummer & Mihalas 1988) in the model equations than those mentioned in Sect. 2.1.1.

6 Conclusions

We have presented new models of line-driven winds of hot stars. The radiative force and the radiative heating term are consistently calculated using statistical equilibrium equations. The statistical equilibrium equations are solved for 15 most abundant elements. Obtained occupation numbers are used for the calculation of the radiative force (in total about 180 000 lines are accounted for in the calculation of the radiative force) and for the calculation of radiative cooling and heating term. All bound-bound collisional transitions and bound-free collisional and radiative transitions explicitly included in the models are allowed to contribute to the cooling and heating terms. Moreover, free-free transitions of H and He are also included.

We test the ability of our new models to predict correct wind parameters of late O stars. In particular, we compare calculated terminal velocities, mass-loss rates and wind momenta with the observed ones. We show that wind parameters calculated using our models are consistent with observation. Moreover, our new wind models have solved the problem introduced by LSL of high theoretical wind terminal velocities. However, this conclusion has to be tested again against more advanced models, because our models involve several simplifying assumptions. Moreover, the relatively good correspondence between observed and theoretical wind parameters can be reduced by our insufficient knowledge of stellar and wind parameters.

There is a good agreement between our calculated wind momentum-luminosity relationship and that relationship observed by Herrero et al. (2002), and Repolust et al. (2004) for non-supergiant stars. However, contrary to the theoretical prediction of Vink et al. (2000) and ours, there is a difference between the wind momentum-luminosity relationship observed for non-supergiant stars and for supergiants (or, more precisely, between stars with a different type of H$\alpha$ profile). Repolust et al. (2004) attribute this difference to wind clumping.

On the other hand, many steps remains to be done. Besides the overall improvements of our code we have to compare observed and calculated wind parameters for B stars before we start to study multicomponent models.

We postpone detailed discussion of wind temperature structure for a separate paper.

Acknowledgements

The authors would like to thank Dr. J. Puls for his valuable comments on the manuscript. This research has made use of NASA's Astrophysics Data System and the SIMBAD database, operated at CDS, Strasbourg, France. This work was supported by a PPARC Rolling Grant and by grants GA CR 205/01/0656 and 205/02/0445. The Astronomical Institute Ondrejov is supported by projects K2043105 and Z1003909.

References

• Abbott, D. C. 1982, ApJ, 259, 282 [NASA ADS] [CrossRef] (In the text)
• Abbott, D. C., & Lucy, L. B. 1985, ApJ, 288, 679 [NASA ADS] [CrossRef] (In the text)
• Abbott, D. C., & Hummer, D. G. 1985, ApJ, 294, 286 [NASA ADS] [CrossRef] (In the text)
• Anderson, E., Bai, Z., Bischof, C., et al. 1999, LAPACK Users' Guide, 3rd ed., SIAM, Philadelphia (In the text)
• Anderson, L. S. 1985, ApJ, 298, 848 [NASA ADS] [CrossRef] (In the text)
• Anderson, L. S. 1989, ApJ, 339, 558 [NASA ADS] [CrossRef] (In the text)
• Babel, J. 1995, A&A, 301, 823 [NASA ADS] (In the text)
• Babel, J. 1996, A&A, 309, 867 [NASA ADS] (In the text)
• Bianchi, L., & Garcia, M. 2002, ApJ, 581, 610 [NASA ADS] [CrossRef] (BG) (In the text)
• Bautista, M. A. 1996, A&AS, 119, 105 [NASA ADS] (In the text)
• Bautista, M. A., & Pradhan, A. K. 1997, A&AS, 126, 365 [NASA ADS] (In the text)
• Bouret, J.-C., Lanz, T., Hillier, D. J., et al. 2003, ApJ, 595, 1182 [NASA ADS] [CrossRef] (In the text)
• Butler, K., Mendoza, C., & Zeippen, C. J. 1993, J. Phys. B, 26, 4409 [NASA ADS] (In the text)
• Cassinelli, J. P., & Olson, G. L. 1979, ApJ, 229, 304 [NASA ADS] [CrossRef] (In the text)
• Castor, J. I. 1974, MNRAS, 169, 279 [NASA ADS] (In the text)
• Castor, J. I., Abbott, D. C., & Klein, R. I. 1975, ApJ, 195, 157 [NASA ADS] [CrossRef] (CAK) (In the text)
• Chen, G. X., & Pradhan, A. K. 1999, A&AS, 136, 395 [NASA ADS] (In the text)
• Crowther, P. A., & Bohannan, B. 1997, A&A, 317, 532 [NASA ADS] (In the text)
• Crowther, P. A., Hillier, D. J., Evans, C. J., et al. 2002, ApJ, 579, 774 [NASA ADS] [CrossRef] (In the text)
• de Koter, A., Schmutz, W., & Lamers, H. J. G. L. M. 1993, A&A, 277, 561 [NASA ADS] (In the text)
• Dessart, L., & Owocki, S. P. 2003, A&A, 406, L1 [EDP Sciences] [NASA ADS] [CrossRef] (In the text)
• Drew, J. E. 1989, ApJS, 71, 267 [NASA ADS] [CrossRef] (In the text)
• Feldmeier, A., Puls, J., & Pauldrach, A. W. A. 1997, A&A, 322, 878 [NASA ADS] (In the text)
• Gayley, K. G. 1995, ApJ, 454, 410 [NASA ADS] [CrossRef] (In the text)
• Gayley, K. G., & Owocki, S. P. 1994, ApJ, 434, 684 [NASA ADS] [CrossRef] (GO) (In the text)
• Gräfener, G. 2003, in Stellar Atmosphere Modelling, ed. I. Hubeny, D. Mihalas, & K. Werner, ASP Conf. Ser., 288, 533 (In the text)
• Gräfener, G., Koesterke, L., & Hamann, W.-R. 2002, A&A, 387, 244 [EDP Sciences] [NASA ADS] [CrossRef] (In the text)
• Herrero, A., Kudritzki, R. P., Vilchez, J. M., et al. 1992, A&A, 261, 209 [NASA ADS] (In the text)
• Herrero, A., Puls, J., & Najarro, F. 2002, A&A, 396, 949 [EDP Sciences] [NASA ADS] [CrossRef] (In the text)
• Hillier, D. J., & Miller, D. L. 1998, ApJ, 496, 407 [NASA ADS] [CrossRef] (In the text)
• Hubeny, I. 1988, Comput. Phys. Commun., 52, 103 [NASA ADS] [CrossRef] (In the text)
• Hubeny, I., & Lanz, T. 1992, A&A, 262, 501 [NASA ADS] (In the text)
• Hubeny, I., & Lanz, T. 1995, ApJ, 439, 875 [NASA ADS] [CrossRef] (In the text)
• Hummer, D. G., Berrington, K. A., Eissner, W., et al. 1993, A&A, 279, 298 [NASA ADS] (In the text)
• Hummer, D. G., & Mihalas, D. 1988, ApJ, 331, 794 [NASA ADS] [CrossRef] (In the text)
• Hummer, D. G., & Rybicki, G. B. 1985, ApJ, 293, 258 [NASA ADS] [CrossRef] (In the text)
• Korcáková, D., & Kubát, J. 2003, A&A, 401, 419 [EDP Sciences] [NASA ADS] [CrossRef] (In the text)
• Kramer, R. H., Cohen, D. H., & Owocki, S. P. 2003, ApJ, 592, 532 [NASA ADS] [CrossRef] (In the text)
• Krticka, J. 2001, Ph.D. Thesis, Masaryk University Brno (In the text)
• Krticka, J. 2003, in Stellar Atmosphere Modelling, ed. I. Hubeny, D. Mihalas, & K. Werner, ASP Conf. Ser., 288, 259 (In the text)
• Krticka, J., & Kubát, J. 2001, A&A, 377, 175 [EDP Sciences] [NASA ADS] [CrossRef] (KKII) (In the text)
• Krticka, J., & Kubát, J. 2002, A&A, 388, 531 [EDP Sciences] [NASA ADS] [CrossRef] (In the text)
• Krticka, J., Owocki, S. P., Kubát, J., Galloway, R. K., & Brown, J. C. 2003, A&A, 402, 713 [EDP Sciences] [NASA ADS] [CrossRef] (In the text)
• Kríz, S., & Harmanec, P. 1975, Bull. Astron. Inst. Czechosl., 26, 65 [NASA ADS] (In the text)
• Kubát, J. 1993, Ph.D. Thesis, Astronomický ústav AV CR, Ondrejov (In the text)
• Kubát, J. 1994, A&A, 287, 179 [NASA ADS] (In the text)
• Kubát, J. 1997, A&A, 326, 277 [NASA ADS]
• Kubát, J. 2003, in Modelling of Stellar Atmospheres, ed. N. E. Piskunov, W. W. Weiss, & D. F. Gray, IAU Symp. 210, ASP, A8 (In the text)
• Kubát, J., Puls, J., & Pauldrach, A. W. A. 1999, A&A, 341, 587 [NASA ADS] (In the text)
• Kudritzki, R. P. 2002, ApJ, 577, 389 [NASA ADS] [CrossRef] (In the text)
• Kudritzki, R. P., Pauldrach, A. W. A., Puls, J., & Abbott, D. C. 1989, A&A, 219, 205 [NASA ADS] (In the text)
• Kudritzki, R. P., & Puls, J. 2000, ARA&A, 38, 613 [NASA ADS] (In the text)
• Kupka, F., Piskunov, N. E., Ryabchikova, T. A., Stempels, H. C., & Weiss, W. W. 1999, A&AS, 138, 119 [NASA ADS] (In the text)
• Lamers, H. J. G. L. M., Haser, S., de Koter, A., & Leitherer, C. 1999, ApJ, 516, 874 [NASA ADS] (In the text)
• Lamers, H. J. G. L. M., & Leitherer, C. 1993, ApJ, 412, 771 [NASA ADS] [CrossRef] (In the text)
• Lamers, H. J. G. L. M., & Morton, D. C. 1976, ApJS, 32, 715 [NASA ADS] [CrossRef] (In the text)
• Lamers, H. J. G. L. M., Snow, T. P., & Lindholm, D. M. 1995, ApJ, 455, 269 [NASA ADS] [CrossRef] (LSL) (In the text)
• Lanz, T., de Koter, A., Hubeny, I., & Heap, S. R. 1996, ApJ, 465, 359 [NASA ADS] [CrossRef] (In the text)
• Lanz, T., & Hubeny, I. 2003, ApJS, 146, 417 [NASA ADS] [CrossRef] (In the text)
• Leitherer, C. 1988, ApJ, 326, 356 [NASA ADS] [CrossRef] (In the text)
• Lucy, L. B. 2001, MNRAS, 326, 95 [NASA ADS] [CrossRef] (In the text)
• Luo, D., & Pradhan, A. K. 1989, J. Phys. B, 22, 3377 [NASA ADS] (In the text)
• MacFarlane, J. J., Cohen, D. H., & Wang, P. 1994, ApJ, 437, 351 [NASA ADS] [CrossRef] (In the text)
• Markova, N., Puls, J., Repolust, T., & Markov, H. 2004, A&A, 413, 693 [EDP Sciences] [NASA ADS] [CrossRef] (In the text)
• Massa, D., Fullerton, A. W., Sonneborn, G., & Hutchings, J. B. 2003, ApJ, 586, 996 [NASA ADS] [CrossRef] (In the text)
• Mihalas, D. 1978, Stellar Atmospheres (San Francisco: Freeman & Co.) (In the text)
• Mihalas, D., & Hummer, D. G. 1974, ApJS, 28, 343 [NASA ADS] [CrossRef] (In the text)
• Nahar, S. N., & Pradhan, A. K. 1996, A&AS, 119, 509 [NASA ADS] (In the text)
• Nahar, S. N., & Pradhan, A. K. 1993, J. Phys. B, 26, 1109 [NASA ADS] (In the text)
• Owocki, S. P., & Puls, J. 1999, ApJ, 510, 355 [NASA ADS] [CrossRef] (In the text)
• Owocki, S. P., & Puls, J. 2002, ApJ, 568, 965 [NASA ADS] [CrossRef] (In the text)
• Pauldrach, A. W. A. 1987, A&A, 183, 295 [NASA ADS] (In the text)
• Pauldrach, A. W. A., Hoffmann, T. L., & Lennon, M. 2001, A&A, 375, 161 [EDP Sciences] [NASA ADS] [CrossRef] (In the text)
• Pauldrach, A. W. A., Kudritzki, R. P., Puls, J., Butler, K., & Hunsinger, J. 1994, A&A, 283, 525 [NASA ADS] (In the text)
• Petrenz, P., & Puls, J. 2000, A&A, 358, 956 [NASA ADS] (In the text)
• Piskunov, N. E., Kupka, F., Ryabchikova, T. A., Weiss, W. W., & Jeffery, C. S. 1995, A&AS, 112, 525 [NASA ADS] (In the text)
• Porter, J. M., & Skouza, B. A. 1999, A&A, 344, 205 [NASA ADS] (In the text)
• Puls, J. 1987, A&A, 184, 227 [NASA ADS] (In the text)
• Puls, J., Kudritzki, R.-P., Herrero, A., et al. 1996, A&A, 305, 171 [NASA ADS] (In the text)
• Puls, J., Springmann, U., & Lennon, M. 2000, A&AS, 141, 23 [NASA ADS] (In the text)
• Repolust, T., Puls, J., & Herrero, A. 2004, A&A, 415, 349 [EDP Sciences] [NASA ADS] [CrossRef] (In the text)
• Runacres, M. C., & Owocki, S. P. 2002, A&A, 381, 1015 [EDP Sciences] [NASA ADS] [CrossRef] (In the text)
• Rybicki, G. B., & Hummer, D. G. 1978, ApJ, 219, 645 [NASA ADS] [MathSciNet] (In the text)
• Sawey, P. M. J., & Berrington, K. A. 1992, J. Phys. B, 25, 1451 [NASA ADS] (In the text)
• Schmutz, W., Abbott, D. C., Russell, R. S., Hamann, W.-R., & Wessolowski, U. 1990, ApJ, 355, 255 [NASA ADS] [CrossRef]
• Scuderi, S., Bonanno, G., di Benedetto, R., Spadaro, D., & Panagia, N. 1992, ApJ, 392, 201 [NASA ADS] [CrossRef] (In the text)
• Scuderi, S., Panagia, N., Stanghellini, C., Trigilio, C., & Umana, G. 1998, A&A, 332, 251 [NASA ADS] (In the text)
• Seaton, M. J. 1987, J. Phys. B, 20, 6363 [NASA ADS] (In the text)
• Seaton, M. J., Zeippen, C. J., Tully, J. A., et al. 1992, Rev. Mex. Astron. Astrofis., 23, 19 [NASA ADS] (In the text)
• Smith, K. C. 1992, Ph.D. Thesis, Univ. London (In the text)
• Smith, K. C., & Dworetsky, M. M. 1988, in Elemental abundance analyses, ed. S. J. Adelman, & T. Lanz (Switzerland: Institut d'Astronomie de l'Univ. de Lausanne), 32 (In the text)
• Springmann, U. W. E., & Pauldrach, A. W. A. 1992, A&A, 262, 515 [NASA ADS] (In the text)
• ud-Doula, A., & Owocki, S. P. 2002, ApJ, 576, 413 [NASA ADS] [CrossRef] (In the text)
• Van Regemorter, H. 1962, ApJ, 136, 906 [NASA ADS] [CrossRef] (In the text)
• Vink, J. S., de Koter, A., & Lamers, H. J. G. L. M. 1999, A&A, 350, 181 [NASA ADS] (In the text)
• Vink, J. S., de Koter, A., & Lamers, H. J. G. L. M. 2000, A&A, 362, 295 [NASA ADS] (In the text)
• Vink, J. S., de Koter, A., & Lamers, H. J. G. L. M. 2001, A&A, 369, 574 [EDP Sciences] [NASA ADS] [CrossRef] (In the text)
• Zhang, H. L. 1996, A&AS, 119, 523 [NASA ADS] (In the text)
• Zhang, H. L., & Pradhan, A. K. 1997, A&AS, 126, 373 [NASA ADS] (In the text)

  
Online Material

   
Appendix A: Atomic data for NLTE calculations

The source of all data used for the calculation of model atoms is given in Table A.1. Most important models are taken from the TLUSTY web page http://tlusty.gsfc.nasa.gov.

Table A.1: Details of atomic models for NLTE calculation.

These atomic models were complemented by data constructed for some light elements using the Opacity Project (hereafter OP; Seaton 1987; Luo & Pradhan 1989; Seaton et al. 1992; Butler et al. 1993; Nahar & Pradhan 1993) database http://vizier.u-strasbg.fr/OP.html. Energy levels, photoionization cross sections and oscillator strengths were taken from OP data. To avoid complicated calculation of numerous resonances (which shall be, in fact, accounted for using solution of radiative transfer equation in moving medium), we applied linear regression of photoionization data. For the collisional ionization we applied the so-called Seaton's formula (Mihalas 1978, Eq. (5.78) therein) with the photoionization cross section at the ionization edge taken from OP data. Finally, for the collisional excitation we used the Van Regemorter formula (Mihalas 1978, Eq. (5.75) therein, Van Regemorter 1960).

Another important atom is iron. For Fe III we used energy levels and oscillator strengths from Nahar & Pradhan (1996), photoionization cross sections from OP (Sawey & Berrington 1992), collision rate coefficients for low-lying levels from Zhang (1996), and the Van Regemorter formula for levels not included in Zhang (1996). Radiative rates for Fe IV and Fe V were calculated using data from Bautista & Pradhan (1997) or Bautista (1996), and collision excitation rates are again calculated using either data from Zhang & Pradhan (1997) or the Van Regemorter formula. Finally, radiative cross sections (both bound-bound and bound-free) for Fe VI are taken from OP data and collisional rate coefficients are taken from Chen & Pradhan (1999).

   
Appendix B: Details of linearization of rates

Solution of NLTE rate equations is performed using Newton-Raphson iterations (the so called linearization). Linearization matrix follows from Eq. (22), where derivatives of P_ij corresponding to the line transitions are evaluated using Eqs. (3), (5), and (14) - (21) as:

$$\frac{\partial P_{ij}}{\partial N_i}=B_{ij}\left[ (1-\beta)\frac{... ...artial N_i}S_{ij}+ \frac{\partial\beta_{\rm c}}{\partial N_i} I_{\rm c}\right],$$ (B.1a)

$$\frac{\partial P_{ij}}{\partial N_j}=B_{ij}\left[ (1-\beta)\frac{... ...artial N_j}S_{ij}+ \frac{\partial\beta_{\rm c}}{\partial N_j} I_{\rm c}\right],$$ (B.1b)

$$\frac{\partial P_{ji}}{\partial N_i}=B_{ji}\left[ (1-\beta)\frac{... ...artial N_i}S_{ij}+ \frac{\partial\beta_{\rm c}}{\partial N_i} I_{\rm c}\right],$$ (B.1c)

$$\frac{\partial P_{ji}}{\partial N_j}=B_{ji}\left[ (1-\beta)\frac{... ...artial N_j}S_{ij}+ \frac{\partial\beta_{\rm c}}{\partial N_j} I_{\rm c}\right],$$ (B.1d)

where

$$\frac{\partial S_{ij}}{\partial N_i}= -\frac{S_{ij}B_{ij}}{N_i B_{ij}-N_j B_{ji}},$$ (B.2a)

$$\frac{\partial S_{ij}}{\partial N_j}= \frac{S_{ij}}{N_j }+\frac{S_{ij}B_{ji}}{N_i B_{ij}-N_j B_{ji}},$$ (B.2b)

$$\frac{\partial\beta}{\partial N_k}= \frac{1}{2}\int_{\mu_*}^{1}{\... ...au_\mu}}{\tau_\mu^2}\right] \frac{\partial\tau_\mu}{\partial N_k}, \quad k=i,j,$$ (B.2c)

$$\frac{\partial\beta_{\rm c}}{\partial N_k}= \frac{1}{2}\int_{-1}^... ...au_\mu}}{\tau_\mu^2}\right] \frac{\partial\tau_\mu}{\partial N_k}, \quad k=i,j,$$ (B.2d)

$$\frac{\partial\tau_\mu}{\partial N_i}=\frac{\tilde\tau_\mu}{g_i},$$ (B.2e)

$$\frac{\partial\tau_\mu}{\partial N_j}=-\frac{\tilde\tau_\mu}{g_j},$$ (B.2f)

$$\tilde\tau_\mu=\frac{\pi e^2}{m_{\rm e}\nu_{ij}} z_{\rm atom}n_{\rm H} g_if_{ij} \frac{r}{{{v_r}}_{{\rm i}} \left(1+\sigma\mu^2\right)}\cdot$$ (B.2g)

Note that the derivatives of collisional rates with respect to N_i are zero. Finally, the derivatives of normalization condition (2) with respect of all relative populations N_i are:

$$% \frac{\partial P_{ki}}{\partial N_i}=1,$$ (B.3)

where k is index of level with the highest occupation number.

Similarly, for the calculation of derivatives of relative number densities N_i with respect to fluid variables $n_{\rm e}$, $\rho_{\rm p}$, $T_{\rm e}$, $\sigma$ and  ${{v_r}}_{{\rm i}}$ in Eq. (24) it is necessary to evaluate quantities  ${\partial P_{ij}}/{\partial h}$.

The derivatives of excitation and deexcitation rates with respect to electron density n_e follow from Eqs. (3), (5), and (8), and are relatively simple,

$$\frac{\partial R_{ij}}{\partial {n}_{{\rm e}}}=\frac{\partial R_{ji}}{\partial {n}_{{\rm e}}} =0,$$ (B.4a)

$$\frac{\partial C_{ij}}{\partial {n}_{{\rm e}}}=\frac{C_{ij}}{{n}_{{\rm e}}},$$ (B.4b)

$$\frac{\partial C_{ji}}{\partial {n}_{{\rm e}}}=\frac{C_{ji}}{{n}_{{\rm e}}},$$ (B.4c)

whereas for the calculation of the derivatives of the ionization rates (Eqs. (7) and (8)) it is necessary to account for the dependency of Saha-Boltzmann ratio on the electron density

$$\frac{\partial R_{ij}}{\partial {n}_{{\rm e}}} =0,$$ (B.5a)

$$\frac{\partial R_{ji}}{\partial {n}_{{\rm e}}} =\frac{R_{ji}}{{n}_{{\rm e}}},$$ (B.5b)

$$\frac{\partial C_{ij}}{\partial {n}_{{\rm e}}}=\frac{C_{ij}}{{n}_{{\rm e}}},$$ (B.5c)

$$\frac{\partial C_{ji}}{\partial {n}_{{\rm e}}}=\frac{2C_{ji}}{{n}_{{\rm e}}}\cdot$$ (B.5d)

The derivatives of the excitation rates with respect to electron temperature are also relatively simple,

  

$$\frac{\partial R_{ij}}{\partial {T}_{{\rm e}}}=\frac{\partial R_{ji}}{\partial {T}_{{\rm e}}} =0,$$ (B.6a)

$$\frac{\partial C_{ij}}{\partial {T}_{{\rm e}}}=n_{\rm e} \frac{\partial \Omega_{ij} (T_{\rm e})}{\partial{T}_{{\rm e}}},$$ (B.6b)

$$\frac{\partial C_{ji}}{\partial {T}_{{\rm e}}}=n_{\rm e} \frac{\p... ...{N_j}\right)^* \frac{\partial \Omega_{ij} (T_{\rm e})}{\partial {T}_{{\rm e}}},$$ (B.6c)

however, for the evaluation of derivatives of ionization rates it is necessary to perform an integration

$$\frac{\partial R_{ij}}{\partial {T}_{{\rm e}}}=0,$$ (B.7a)

$$\frac{\partial R_{ji}}{\partial {T}_{{\rm e}}}= 4\pi \frac{\parti... ...\frac{2h\nu^3}{c^2}+J_\nu\right] {\rm e}^{-\frac{h\nu}{kT_{\rm e}}} ~{\rm d}\nu$$

$$\qquad + 4\pi \left(\frac{N_i}{N_k}\right)^*\int_{\nu_i}^{\infty}... ...^2}+J_\nu\right]{\rm e}^{-h\nu/kT_{\rm e}} \frac{{\rm d}\nu}{kT_{\rm e}^2}\cdot$$ (B.7b)

Derivatives of collisional rates are given by Eqs. (B.6b) and (B.6c).

Using the dependence of the radiative transfer equation in Sobolev approximation on the velocity gradient $\sigma$ and on  ${{v_r}}_{{\rm i}}$ and  $\rho_{\rm p}$ via the fraction $y={\rho}_{{\rm p}}/{{v_r}}_{{\rm i}}$ we calculate derivatives of excitation rates with respect to these variables, namely

$$\frac{\partial R_{ij}}{\partial \sigma}=B_{ij}\left[ \frac{\parti... ...\partial \sigma} I_{\rm c}- \frac{\partial\beta}{\partial \sigma}S_{ij}\right],$$ (B.8a)

$$\frac{\partial R_{ji}}{\partial \sigma}=B_{ji}\left[ \frac{\parti... ...\partial \sigma} I_{\rm c}- \frac{\partial\beta}{\partial \sigma}S_{ij}\right],$$ (B.8b)

$$\frac{\partial R_{ij}}{\partial {\rho}_{{\rm p}}}=\frac{1}{{{v_r}}_{{\rm i}} }\frac{\partial R_{ij}}{\partial y},$$ (B.8c)

$$\frac{\partial R_{ji}}{\partial {\rho}_{{\rm p}}}=\frac{1}{{{v_r}}_{{\rm i}} }\frac{\partial R_{ji}}{\partial y},$$ (B.8d)

$$\frac{\partial R_{ij}}{\partial {{v_r}}_{{\rm i}} }=-\frac{y}{{{v_r}}_{{\rm i}} }\frac{\partial R_{ij}}{\partial y},$$ (B.8e)

$$\frac{\partial R_{ji}}{\partial {{v_r}}_{{\rm i}} }=-\frac{y}{{{v_r}}_{{\rm i}} }\frac{\partial R_{ji}}{\partial y},$$ (B.8f)

where

$$\frac{\partial\beta}{\partial \sigma}= \frac{1}{2}\int_{\mu_*}^{1... ...\mu^2} \left[\frac{1-{\rm e}^{-\tau_\mu}}{\tau_\mu}-{\rm e}^{-\tau_\mu}\right],$$ (B.9a)

$$\frac{\partial\beta_{\rm c}}{\partial \sigma}= \frac{1}{2}\int_{-... ...\mu^2} \left[\frac{1-{\rm e}^{-\tau_\mu}}{\tau_\mu}-{\rm e}^{-\tau_\mu}\right],$$ (B.9b)

$$\frac{\partial R_{ij}}{\partial y}=B_{ij}\left[ \frac{\partial\beta_{\rm c}}{\partial y} I_{\rm c}- \frac{\partial\beta}{\partial y}S_{ij}\right],$$ (B.9c)

$$\frac{\partial R_{ji}}{\partial y}=B_{ji}\left[ \frac{\partial\beta_{\rm c}}{\partial y} I_{\rm c}- \frac{\partial\beta}{\partial y}S_{ij}\right],$$ (B.9d)

$$\frac{\partial\beta}{\partial y}= \frac{1}{2}\int_{\mu_*}^{1}{\rm... ...[{\rm e}^{-\tau_\mu}- \frac{1-{\rm e}^{-\tau_\mu}}{\tau_\mu}\right]\frac{1}{y},$$ (B.9e)

$$\frac{\partial\beta_{\rm c}}{\partial y}= \frac{1}{2}\int_{-1}^{1... ...m e}^{-\tau_\mu}- \frac{1-{\rm e}^{-\tau_\mu}}{\tau_\mu}\right]\frac{1}{y}\cdot$$ (B.9f)

Note that the variable $\sigma$ is not an indepentend variable during linearization. Derivatives with respect to $\sigma$ are transformed to derivatives of ionic velocity using Eq. (20).

   
Appendix C: Details of linearization of raditive heating/cooling term

Radiative heating/cooling term in the case of thermal balance of electrons is given by (Kubát et al. 1999)

$$% Q^{{\rm rad}}=Q_{{\rm ff}}^{\rm H}+Q_{{\rm bf}}^{\rm H}+Q_{... ...Q_{{\rm ff}}^{\rm C}-Q_{{\rm bf}}^{\rm C}-Q_{{\rm c}}^{\rm C},$$ (C.1)

where individual terms are

 

$$Q_{{\rm ff}}^{\rm H}\left[n_i\right] = 4\pi {n}_{{\rm e}} \sum_{i... ..._{i}\int_{0}^{\infty} \alpha_{{\rm ff},i}(\nu,{T}_{{\rm e}}) J_{\nu}{\rm d}\nu,$$ (C.2a)

$$Q_{{\rm ff}}^{\rm C}\left[n_i\right] = 4\pi {n}_{{\rm e}} \sum_{i... ... ii}, He {\sc iii}} \int_{0}^{\infty}n_{i} \alpha_{{\rm ff}}(\nu,{T}_{{\rm e}})$$

$$\qquad\quad \times \left(J_{\nu} + \frac{2h\nu^3}{c^2}\right) {\rm e}^{-{h\nu}/{k{T}_{{\rm e}}}}{\rm d}\nu,$$ (C.2b)

$$Q_{{\rm bf}}^{\rm H}\left[n_i\right] = 4\pi \sum_{i} n_i \int_{\nu_i}^{\infty} \alpha_{i,\nu} J_{\nu}\left(1-\frac{\nu_{i}}{\nu}\right){\rm d}\nu,$$ (C.2c)

$$Q_{{\rm bf}}^{\rm C}\left[n_{k_i}\right] = 4\pi \sum_{i} \int_{\n... ..._i}{N_{k_i}}\right)^* \alpha_{i,\nu} \left(J_{\nu} + \frac{2h\nu^3}{c^2}\right)$$

$$\qquad\quad\; \times {\rm e}^{-h\nu/k{T}_{{\rm e}}} \left(1-\frac{\nu_{i}}{\nu}\right){\rm d}\nu,$$ (C.2d)

$$Q_{{\rm c}}^{\rm H}\left[n_j\right] = {n}_{{\rm e}} \sum_{ij} n_j \left(\frac{N_i}{N_j}\right)^* \Omega_{ij} (T_{\rm e}) h\nu_{ij},$$ (C.2e)

$$Q_{{\rm c}}^{\rm C}\left[n_i\right] = {n}_{{\rm e}} \sum_{ij} n_i \Omega_{ij} (T_{\rm e}) h\nu_{ij},$$ (C.2f)

where k_i denotes level into which is level i ionized. From these equations derivatives of heating/cooling term with respect to  ${{v_r}}_{{\rm i}}$, $\rho_{\rm p}$, ${\rho}_{{\rm e}}$, $T_{\rm e}$ for the Newton-Raphson iteration step of hydrodynamical structure can be calculated,

$$\frac{\partial Q_{{\rm ff}}^{\rm H}}{\partial {T}_{{\rm e}}} = Q_... ...\rm e}}}\right]- \frac{1}{2{T}_{{\rm e}}} Q_{{\rm ff}}^{\rm H}\left[n_i\right],$$ (C.3a)

$$\frac{\partial Q_{{\rm ff}}^{\rm H}}{\partial {n}_{{\rm e}}} = \f... ...Q_{{\rm ff}}^{\rm H}\left[\frac{\partial n_{i}}{\partial {n}_{{\rm e}}}\right],$$ (C.3b)

$$\frac{\partial Q_{{\rm ff}}^{\rm H}}{\partial {\rho}_{{\rm p}}} =... ...{\rm ff}}^{\rm H}\left[\frac{\partial n_{i}}{\partial {\rho}_{{\rm p}}}\right],$$ (C.3c)

$$\frac{\partial Q_{{\rm ff}}^{\rm H}}{\partial {{v_r}}_{{\rm i}} }... ...rm ff}}^{\rm H}\left[\frac{\partial n_{i}}{\partial {{v_r}}_{{\rm i}} }\right],$$ (C.3d)

$$\frac{\partial Q_{{\rm ff}}^{\rm H}}{\partial \sigma} = Q_{{\rm ff}}^{\rm H}\left[\frac{\partial n_{i}}{\partial\sigma}\right],$$ (C.3e)

$$\frac{\partial Q_{{\rm ff}}^{\rm C}}{\partial {T}_{{\rm e}}}= Q_{... ...u n_{i}\right]- \frac{1}{2{T}_{{\rm e}}}Q_{{\rm ff}}^{\rm C}\left[n_{i}\right],$$ (C.3f)

$$\frac{\partial Q_{{\rm ff}}^{\rm C}}{\partial {n}_{{\rm e}}}= \fr... ...Q_{{\rm ff}}^{\rm C}\left[\frac{\partial n_{i}}{\partial {n}_{{\rm e}}}\right],$$ (C.3g)

$$\frac{\partial Q_{{\rm ff}}^{\rm C}}{\partial{\rho}_{{\rm p}}}= Q_{{\rm ff}}^{\rm C}\left[\frac{\partial n_{i}}{\partial {\rho}_{{\rm p}}}\right],$$ (C.3h)

$$\frac{\partial Q_{{\rm ff}}^{\rm C}}{\partial {{v_r}}_{{\rm i}} }... ...rm ff}}^{\rm C}\left[\frac{\partial n_{i}}{\partial {{v_r}}_{{\rm i}} }\right],$$ (C.3i)

$$\frac{\partial Q_{{\rm ff}}^{\rm C}}{\partial \sigma} = Q_{{\rm ff}}^{\rm C}\left[\frac{\partial n_{i}}{\partial\sigma}\right],$$ (C.3j)

$$\frac{\partial Q_{{\rm bf}}^{\rm H}}{\partial {T}_{{\rm e}}} = Q_{{\rm bf}}^{\rm H}\left[\frac{\partial n_{i}}{\partial {T}_{{\rm e}}}\right],$$ (C.3k)

$$\frac{\partial Q_{{\rm bf}}^{\rm H}}{\partial {n}_{{\rm e}}} = Q_{{\rm bf}}^{\rm H}\left[\frac{\partial n_{i}}{\partial {n}_{{\rm e}}}\right],$$ (C.3l)

$$\frac{\partial Q_{{\rm bf}}^{\rm H}}{\partial {\rho}_{{\rm p}}} =... ...{\rm bf}}^{\rm H}\left[\frac{\partial n_{i}}{\partial {\rho}_{{\rm p}}}\right],$$ (C.3m)

$$\frac{\partial Q_{{\rm bf}}^{\rm H}}{\partial {{v_r}}_{{\rm i}} }... ...rm bf}}^{\rm H}\left[\frac{\partial n_{i}}{\partial {{v_r}}_{{\rm i}} }\right],$$ (C.3n)

$$\frac{\partial Q_{{\rm bf}}^{\rm H}}{\partial \sigma} = Q_{{\rm bf}}^{\rm H}\left[\frac{\partial n_{i}}{\partial\sigma}\right],$$ (C.3o)

$$\frac{\partial Q_{{\rm bf}}^{\rm C}}{\partial {T}_{{\rm e}}}= Q_{... ...}}\right]- \frac{h}{k{T}_{{\rm e}}}Q_{{\rm bf}}^{\rm C}\left[\nu n_{k_i}\right]$$

$$\qquad + Q_{{\rm bf}}^{\rm C}\left[n_{k_i}\left\{\left(\frac{N_i}... ...ac{\partial}{\partial {T}_{{\rm e}}} \left(\frac{N_i}{N_{k_i}}\right)^*\right],$$ (C.3p)

$$\frac{\partial Q_{{\rm bf}}^{\rm C}}{\partial {n}_{{\rm e}}}= Q_{{\rm bf}}^{\rm C}\left[\frac{\partial n_{k_i}}{\partial{n}_{{\rm e}}}\right]$$

$$\qquad + Q_{{\rm bf}}^{\rm C}\left[n_{k_i}\left\{\left(\frac{N_i}... ...ac{\partial}{\partial {n}_{{\rm e}}} \left(\frac{N_i}{N_{k_i}}\right)^*\right],$$ (C.3q)

$$\frac{\partial Q_{{\rm bf}}^{\rm C}}{\partial{\rho}_{{\rm p}}}= Q... ...rm bf}}^{\rm C}\left[\frac{\partial n_{k_i}}{\partial {\rho}_{{\rm p}}}\right],$$ (C.3r)

$$\frac{\partial Q_{{\rm bf}}^{\rm C}}{\partial {{v_r}}_{{\rm i}} }... ... bf}}^{\rm C}\left[\frac{\partial n_{k_i}}{\partial {{v_r}}_{{\rm i}} }\right],$$ (C.3s)

$$\frac{\partial Q_{{\rm bf}}^{\rm C}}{\partial \sigma} = Q_{{\rm bf}}^{\rm C}\left[\frac{\partial n_{k_i}}{\partial\sigma}\right],$$ (C.3t)

$$\frac{\partial Q_{{\rm c}}^{\rm H}}{\partial {T}_{{\rm e}}}= Q_{{... ...}_{{\rm e}}} \left(\frac{N_i}{N_j}\right)^* \Omega_{ij}({T}_{{\rm e}})h\nu_{ij}$$

$$\qquad + {n}_{{\rm e}} \sum_{ij} n_j\left(\frac{N_i}{N_j}\right)^* \frac{\partial\Omega_{ij}({T}_{{\rm e}})}{\partial {T}_{{\rm e}}}h\nu_{ij},$$ (C.3u)

$$\frac{\partial Q_{{\rm c}}^{\rm H}}{\partial {n}_{{\rm e}}}= \fra... ...+ Q_{{\rm c}}^{\rm H}\left[\frac{\partial n_{j}}{\partial {n}_{{\rm e}}}\right]$$

$$\qquad + {n}_{{\rm e}} \sum_{ij} n_j\frac{\partial}{\partial {n}_{{\rm e}}} \left(\frac{N_i}{N_j}\right)^* \Omega_{ij}({T}_{{\rm e}})h\nu_{ij},$$ (C.3v)

$$\frac{\partial Q_{{\rm c}}^{\rm H}}{\partial {\rho}_{{\rm p}}}= Q_{{\rm c}}^{\rm H}\left[\frac{\partial n_{j}}{\partial {\rho}_{{\rm p}}}\right],$$ (C.3w)

$$\frac{\partial Q_{{\rm c}}^{\rm H}}{\partial{{v_r}}_{{\rm i}} }= Q_{{\rm c}}^{\rm H}\left[\frac{\partial n_{j}}{\partial{{v_r}}_{{\rm i}} }\right],$$ (C.3x)

$$\frac{\partial Q_{{\rm c}}^{\rm H}}{\partial\sigma}= Q_{{\rm c}}^{\rm H}\left[\frac{\partial n_{j}}{\partial\sigma}\right],$$ (C.3y)

$$\frac{\partial Q_{{\rm c}}^{\rm C}}{\partial {T}_{{\rm e}}}= Q_{{... ...n_i \frac{\partial\Omega_{ij}({T}_{{\rm e}})}{\partial {T}_{{\rm e}}}h\nu_{ij},$$ (C.3z)

$$\frac{\partial Q_{{\rm c}}^{\rm C}}{\partial {n}_{{\rm e}}}= \fra... ... Q_{{\rm c}}^{\rm C}\left[\frac{\partial n_{i}}{\partial {n}_{{\rm e}}}\right],$$ (C.3A)

$$\frac{\partial Q_{{\rm c}}^{\rm C}}{\partial {\rho}_{{\rm p}}}= Q_{{\rm c}}^{\rm C}\left[\frac{\partial n_{i}}{\partial {\rho}_{{\rm p}}}\right],$$ (C.3B)

$$\frac{\partial Q_{{\rm c}}^{\rm C}}{\partial{{v_r}}_{{\rm i}} }= Q_{{\rm c}}^{\rm C}\left[\frac{\partial n_{i}}{\partial{{v_r}}_{{\rm i}} }\right],$$ (C.3C)

$$\frac{\partial Q_{{\rm c}}^{\rm C}}{\partial\sigma}= Q_{{\rm c}}^{\rm C}\left[\frac{\partial n_{i}}{\partial\sigma}\right],$$ (C.3D)

where

$$\frac{\partial n_{i}}{\partial x}=\frac{n_i}{N_i} \frac{\partial ... ...\qquad {\rm where}\;x={T}_{{\rm e}}, {n}_{{\rm e}}, {{v_r}}_{{\rm i}} , \sigma,$$ (C.4a)

$$\frac{\partial n_{i}}{\partial{\rho}_{{\rm p}}}=\frac{n_i}{{\rho}_{{\rm p}}}+ \frac{n_i}{N_i}\frac{\partial N_{i}}{\partial{\rho}_{{\rm p}}}\cdot$$ (C.4b)