A&A 458, 597-608 (2006)
DOI: 10.1051/0004-6361:20054666

2D non-LTE modeling for axisymmetric winds

Method and test cases

L. N. Georgiev^1,2 - D. J. Hillier² - J. Zsargó²

1 - Instituto de Astronomía, UNAM, CD. Universitaria, Apartado Postal 70-264, 04510 México DF, México
2 - Physics and Astronomy Department, University of Pittsburgh, 3941 O'Hara Street, Pittsburgh PA 15260, USA

Received 8 December 2005 / Accepted 28 June 2006

Abstract
Aims. This paper describes a code which is capable of modeling axisymmetric stellar winds.
Methods. At present the code uses a short characteristic technique for the continuum transfer, while bound-bound transitions are treated using the Sobolev approximation. The simultaneous solution of the transfer equations, and the equations of statistical equilibrium, is achieved by a combination of preconditioning and linearization.
Results. Extensive tests were performed which show that, in the extreme conditions of gray opacity or spherical symmetry, the code performs similar to other established codes. Simple 2D tests have been undertaken which confirm the performance of the code, and serve to illustrate simple effects that can occur in axisymmetric stellar winds. The code has been designed with flexibility in mind, with the radiation transport section treated, as far as possible, distinctly from the solution of the statistical equilibrium and energy balance equations. This will facilitate the inclusion of other transport solvers, and the treatment of more complex geometries.

Key words: stars: atmospheres - stars: early-type - stars: winds, outflows - line: formation - radiative transfer

1 Introduction

Massive stars are one of the principal sources of heavy elements in galaxies. These stars have stellar winds which expel the outer layers of the star, mixing their material with the interstellar gas. Recent developments in the theory of stellar evolution (Maeder & Meynet 2003) show that rotation has a major influence on the evolution of a star and the chemical composition of its atmosphere. Further, rotation can deform the star, deviating its atmosphere and wind from spherical symmetry (Bjorkman & Cassinelli 1993; Owocki et al. 1994). The influence of deviations from spherical geometry on observed spectra are complicated, and their effect on results obtained using quantitative 1D detailed modeling is unclear. Although difficult, multi-dimensional models are essential for reliable and quantitative analysis of real stars.

The need for multi-dimensional models also arises in interacting binary systems. In O star binary systems, wind-wind collisions occur and these can have a substantial effect on the observed wind spectra (Hill et al. 2002). The analysis of the wind eclipses in WR+O binaries (Auer & Koenigsberger 1994) is a promising tool for studying the wind structure. However the structure of the primary wind is influenced by the radiation field of the secondary star, and thus its correct application requires 2D modeling. The Hatchett-McCray effect (Hatchett & McCray 1977) is another example of an interaction between stars in a binary system, and for which correct interpretation is impossible with 1D codes. Even radial velocity curves of massive binaries are affected by profile changes due to the deviation from spherical symmetry, and these changes need to be taken into account if accurate mass ratios are to be obtained.

The importance of multi-dimensional modeling has been known since the beginning of computer-based stellar atmosphere modeling. But the straight forward application of 1D techniques to 2D problems was not possible and the development of codes was delayed. Several important developments have now made progress possible. These include the dramatic increase in computational power, the availability of the necessary atomic data, and new computational techniques. For example, Mihalas et al. (1978) and Auer & Paletou (1994) introduced the short characteristic method for solution of the equation of radiative transport (ERT). The method is flexible and can be generalized for different geometries and dimensions (e.g., van Noort et al. 2002). Olson et al. (1986) introduced the Approximate Lambda Operator (ALO) to overcome the slow convergence of the usual Lambda iteration. The method combines the simplicity of the Lambda iteration with good convergence properties. Variants of this technique, which use linearization, have also been shown to be extremely useful for construction of 1D non-LTE model atmospheres (Hillier 1990; Hubeny & Lanz 1995). Finally, Rybicki & Hummer (1991) combined the ALO with preconditioning of the statistical equilibrium equations (SEE) which replaces the linearization technique.

Other codes, and techniques, are being developed to perform multi-dimensional radiative transfer. For example, Monte-Carlo techniques have also been utilized to perform multidimensional radiative transfer calculations. Both Juvela (1997) and Park & Hong (1998) have used Monte-Carlo techniques to model CS molecule spectra produced by clumpy cloud models. Monte-Carlo codes have also been extensively applied to the study SEDs of protostars (Whitney et al. 2003b,a). Despite the use of the MC approach, techniques have been developed that allow the temperature structure of the envelope to be solved for (Lucy 1999a; Bjorkman & Wood 2001). Most, if not all, the Monte-Carlo codes use the Sobolev approximation for the line transfer. A major proponent of the Monte-Carlo technique is Lucy, who has devoted considerable effort to its development (e.g., Lucy 1999a,b, 2003).

Given the complexity of non-LTE multi-dimensional radiative transfer, and the wide variety of astrophysical situations where they can be applied, it is important to develop independent codes and techniques. With such independent codes, issues of accuracy and efficiency can be addressed. Importantly, the absolute accuracy of results can be be discerned by a detailed comparison between results obtained using independent codes.

In this paper we describe an early version of a non-LTE radiative transfer code that can be used to undertake detailed spectral analysis in 2D, and in particular spectral analysis of stars with non-spherical stellar winds. The code utilizes several of the numerical techniques mentioned previously. It solves self consistently the ERT simultaneously with the system of SEE and the energy conservation equation. At present we limit ourselves to models with axisymmetry although we allow for the influence of rotation. The solution of the continuum ERT is similar to the one described by Busche & Hillier (2000).

The paper is organized as follows. Section 2 describes the coordinate system and the method used to solve the ERT. In Sect. 3 we discuss the solution of the system of statistical equilibrium equations, the use of the Sobolev approximation for bound-bound transitions, and the electron thermal balance (ETB). The technique we use for the simultaneous solution of the ERT, SEE, and ETB is also discussed. Section 4 presents a comparison of 1D models with those of CMFGEN, and detailed results for simple 2D models.

2 Solution of the radiative transfer equation

The solution of the radiative transfer equation follows closely the method described by Busche & Hillier (2000). For a more consistent description of the code we briefly summarize the main ideas.

The parameters of the gas are described in two coordinate systems. One is a Cartesian system (x y z) with the z axis along the symmetry axis. The second system is spherical ( $r \beta \gamma$), defined so that the angle $\beta$ is measured from the z axis (Fig. 1). Due to the symmetry, variables describing the conditions of the gas and the radiation do not depend on $\gamma$. In the spherical system we define a grid $(r_i,\beta_j)$ (where $1 \le i \le N_r$ and $1 \le j \le N_\beta$). As usual, r_i > r_i+1 and $\beta_j < \beta_{j+1}$. The points with $\beta =0$ lie on the symmetry axis (the pole) and $\beta = \pi/2$ corresponds to the equator. At each point O we define a third local coordinate system $(\vec i,\vec j,\vec k)_{\rm O}$ so that $\vec k \vert\vert \vec r$ (Fig. 1). In this system, a unitary vector $\vec n~=~(\cos{\phi} \sin{\theta},\sin{\phi} \sin{\theta},\cos{\theta})$ describes the direction at which the radiation is propagated. At each point $r_i,\beta_j$ we define a second local grid of directions $\mu_k,\phi_l$ (where $1 \le l \le N_\mu$ and $1 \le k \le N_\phi$, $\mu = \cos \theta$). The angles $\phi_k$ are chosen for a Gaussian integration in the interval $[0,\pi]$. $N_\phi$ is odd, so the direction $\phi = \pi/2$ is always present. The $\mu$ grid is chosen as in the spherically symmetric codes. A set of impact parameters is calculated as follows

$$p_m = \left\{ \begin{array}{l@{\quad:\quad}l} r_{\min} \left... ... \\ r_{m-N_{\rm c}+1} & m \ge N_{\rm c} \end{array} \right.$$ (1)

where $\alpha$ is a parameter which defines how dense are the impact parameters near the edge of the core and $N_{\rm c}$ is the number of impact parameters in the core. Then for each r_i the grid points $\mu_{ik}$ are calculated by $\mu_{ik} = \sqrt{1-p_k^2/r_i^2}$. In that way, a line starting at r_i in the direction $\mu_{ik}$ crosses the shell with radius r_i-1 at angle $\mu_{i-1,k}$. The grid $\mu_k,\phi_l$is the same for all points with the same radius r_i.

	Figure 1: Definition of the coordinate systems.
Open with DEXTER

The specific intensity I_i,j,k,l at a point O with coordinates $(r_i,\beta_j)$ and at the direction $(\mu_k,\phi_l)$is obtained from the formal solution of the ERT along the segment OM.

$$I_{i,j,k,l} = I_{\rm M} \exp(-\Delta \tau) + \int_0^{\Delta \tau} S(t) {\rm e}^{-t} {\rm d}t$$ (2)

where S and $\chi$ are, respectively, the source function and the opacity along the line, and $\Delta \tau$ is the optical depth of the segment OM

$$\Delta \tau = \int_{\rm O}^{\rm M} \chi(l) ~ {\rm d}l .$$ (3)

For the radiation propagated from the outer boundary inward, the segment starts on the shell r_i at the point O and finishes at the point M on the shell with radius r_i-1. We define P to be a point on the same line but lying at the shell with radius r_i+1. If the length of the segment OM is small with respect to the characteristic length scale of the wind, then one can interpolate the source function S and the opacity $\chi$along the line MOP with second order polynomials and rewrite (2) in the form

$$I_{i,j,k,l} = I_{\rm M} \exp(-\Delta \tau) + S_{\rm P}~w_{\rm P} + S_{\rm O}~w_{\rm O} + S_{\rm M}~w_{\rm M}.$$ (4)

The algorithm for calculating the coordinates of the points M and P, the weight w and the interpolation of S, $\chi$ and the intensity I is given in Auer & Paletou (1994), Dullemond & Turolla (2000), Busche & Hillier (2000) and van Noort et al. (2002).

At the points M and P only r and $\mu$ coincide with the predefined grid. Consequently, the source function and the opacity need to be interpolated in the angle $\beta$. For this we use monotonic cubic polynomials as in Busche & Hillier (2000). The intensity depends on all four coordinates and we used bilinear interpolation on $\beta$ and $\phi$. Linear interpolation is known to introduce a numerical scatter (Kunasz & Auer 1988) and can introduce inaccuracies in the flux integrals. However, it is simple and guarantees the positivity of the calculated intensity. At this early stage of the development of the code, we prefer the simple techniques but we keep in mind the possible drawbacks when we analyze the accuracy of the results. More sophisticated techniques can be easily incorporated at a later stage.

If $I(r,\beta,\mu,\phi)$ is known, the mean intensity can be calculated by

 

$$J_\nu \frac{1}{4\pi} \int_0^\pi {\rm d}\phi \int_{-1}^{1} I_\nu(\mu,\phi) {\rm d}\mu$$ = (5)

$$\frac{1}{4\pi} \int_0^\pi {\rm d}\phi \int_{-1}^{1}{\rm d}\mu I_{\rm M}(\mu,\phi) ~ {\rm e}^{-\Delta \tau\left(\mu,\phi\right)}$$ =

$$+\frac{1}{4\pi} \int_0^\pi {\rm d}\phi \int_{-1}^{1}{\rm d}\mu S_{\rm M}(\mu,\phi) w_{\rm M}(\mu,\phi)$$

$$+\frac{1}{4\pi} S_{\rm O} \int_0^\pi {\rm d}\phi \int_{-1}^{1}{\rm d}\mu w_{\rm O}(\mu,\phi)$$

$$+\frac{1}{4\pi} \int_0^\pi {\rm d}\phi \int_{-1}^{1}{\rm d}\mu S_{\rm P}(\mu,\phi) w_{\rm P}(\mu,\phi)$$ (6)

$$\Lambda[S].$$ = (7)

The ALO method splits the $\Lambda$ operator in two: $\Lambda = \Lambda^* + (\Lambda-\Lambda^*)$where $\Lambda^*$ is chosen such that it is easily computed, and is easily inverted. Following Olson et al. (1986), $\Lambda^*$ can be calculated assuming $S_{\rm O} = 1$ and $S_{\rm M} = S_{\rm P} =0$ for all values of $\mu$ and $\phi$. As this operator is local, it is both easy to compute and invert. Applying this method to (7) one obtains

$$\Lambda^*(r_i,\beta_j)= \frac{1}{4\pi} \int_0^\pi {\rm d}\phi \int_{-1}^{1} w_{\rm O}(r_i,\beta_j,\mu,\phi) {\rm d}\mu .$$ (8)

The algorithm to compute $J_\nu$ is as follows. The grid $r_i,\beta_j$ is constructed in the way it samples the density structure of the wind. Then at each point $(r_i,\beta_j)$, the local grid $\mu_k,\phi_l$ is constructed. The lengths of all characteristics on this grid are then calculated, together with the coordinates of their end points as well as the weights for the interpolation of S, $\chi$ and I. All these values depend only on the chosen grids and do not change from iteration to iteration.

Once the grids are ready, the total opacity $\chi_\nu$and the total source function $S_\nu$ are calculated from the current level populations for a given frequency $\nu$. The algorithm calculates the optical depths for each characteristic and then calculates the intensity $I_\nu(r_i,\beta_j,\mu_k,\phi_l)$ shell by shell, and angle by angle, starting from the outer boundary inward and then from the inner boundary outward. At the same time the mean intensity $J_\nu$ and the flux  $\vec{H_\nu}$ are calculated by

   

$$J_\nu \frac{1}{4\pi} \int_0^\pi {\rm d}\phi \int_{-1}^{1} I_\nu(\mu,\phi) {\rm d}\mu$$ = (9)

$$H^r_\nu \frac{1}{4\pi} \int_0^\pi {\rm d}\phi \int_{-1}^{1} I_\nu(\mu,\phi) \mu {\rm d}\mu$$ = (10)

$$H^\beta_\nu \frac{1}{4\pi} \int_0^\pi {\rm d}\phi \int_{-1}^{1} I_\nu(\mu,\phi) \cos{\phi} \sqrt{1-\mu^2}{\rm d}\mu$$ = (11)

where $H^r_\nu,H^\beta_\nu$ are the r and $\beta$ components of the flux-like variable $\vec{H_\nu}$. The integrals on $\mu$are calculated assuming linear dependence of I on $\mu$ and the integrals over $\phi$ are calculated with Gaussian weights.

The final step is the definition of the boundary conditions. The outer boundary is simple. For a single isolated star $I^-(r_1,\beta_j) = 0$. The axial symmetry of the wind allows us to model stars in binary systems. In that case, the intensity at the outer boundary is set by the illumination of the companion.

The inner boundary is more complicated. The diffusion approximation is usually assumed in spherically symmetric atmospheres but if the star is distorted by rotation, the flux at the surface of the stellar core will depend on latitude (Maeder 1999). There is no simple model that describes how the interface between the stellar core and the base of the wind change with rotation nor how the wind hydrodynamics reacts to this change. To keep the model simple, we assume that in the deep layers (high optical depth) the grid points with the same r_i interchange only a small amount of information in the tangential ($\beta$) direction. That is, we assume that the $H_\beta$ component vanishes, and the usual diffusion approximation is independently valid for each $\beta$:

$$I^+\left(r_{\min},\beta,\mu\right) = S\left(r_{\min},\beta\ri... ...\beta,\mu\right)}{\partial \tau}\right\vert _{r=r_{\min}}\cdot$$ (12)

If the opacity $\chi_\nu$ is pure gray scattering then $\left.\frac{\partial S}{\partial \tau}\right\vert _{r=r_{\min}} = 3H$where H is the given total flux emitted by the star. When $\chi_\nu$ depends on the frequency $\nu$, and if the radiation field tends to its local thermodynamic equilibrium (LTE) value, the usual diffusion approximation

$$S_\nu(r_{\min}) = B_\nu(T(r_{\min})) \\$$ (13)

and

$$\left.\frac{\partial S}{\partial \tau}\right\vert _{r=r_{\min... ...eft.\frac{\partial B_\nu}{\partial T}\right\vert _{r=r_{\min}}$$ (14)

is used. Here $B_\nu(T)$ is the Planck function, $\chi_{\rm ross}$ is the Rosseland mean opacity and H is the total flux defined by the stellar luminosity, ${\mathcal L}$, as $H = {\mathcal L}/16\pi^2 r_{\min}^2$. In the present models we have assumed, for simplicity, that at the inner boundary the flux H is independent of $\beta$ and that all the asymmetries in the radiation field arise from the asymmetries in the density structure of the wind. The luminosity is a parameter of the model.

At the end of this procedure, we have the radiation field $I_\nu$, its moments $J_\nu$ and $\vec{H_\nu}$ and the approximate lambda operator $\Lambda^*_\nu$at each point of the grid. One can look at this procedure as a $\Lambda$ operator applied on a known source function $S_\nu$.

3 Solution of the statistical equilibrium equation (SEE)

The radiation field and the thermodynamical conditions of the gas are strongly coupled in the stellar atmosphere. The transport of the radiation is determined by the opacity and the emissivity of the gas which are determined by statistical equilibrium and energy conservation which, in turn, depend on the radiation field. All of these processes work together so the system of the equations which describes them should be solved simultaneously. Unfortunately, the system is highly nonlinear, which makes its solution both difficult and unstable.

Before proceeding further we need to decide how to treat the energy conservation. There are at least two possibilities. The first one is the use of the equation of the radiative equilibrium (ERE). It balances the total absorbed energy with the total emitted energy at each point of the atmosphere. The second possibility is the balance of the electron kinetic energy. At each point the amount of energy transferred to the kinetic motion of the free electrons is balanced by the kinetic energy converted to radiation and therefore lost from the thermal pool. This process is known as electron thermal balance (ETB) (Kubát et al. 1999). The two methods are physically and mathematically equivalent (Hillier 2003), but the numerical application of each of them is different. For example, the collisional terms appear explicitly in the ETB method, but only appear in the ERE method through the coupling with the SE equations.

The ETB method has been used for years in nebula work. It has the advantage that it directly cancels some of the radiation scattering terms. In the present code, we treat the bound-bound transitions in the Sobolev approximation (see below). These terms do not appear explicitly in the ETB equation and we adopt this method of treating the energy balance. However, one has to keep in mind that the ETB method will cause problems at depth. At depth the collisional terms are actually irrelevant and this may cause cancellation problems.

The complete system of equations to be solved is

    

$$J_\nu \Lambda\left[S_\nu\right]$$ = (15)

$$S_\nu \frac{\sigma}{\chi_{\rm tot}} J_\nu + \frac{\eta_{\rm therm}}{\chi_{\rm tot}}$$ = (16)

$$-\sum_{j<i} \Gamma_{ji} + \sum_{j>i} \Gamma_{ij} - \sum_{j<i} \Xi_{ji} + \sum_{j>i} \Xi_{ij}$$ 0 = (17)

$$\Delta Q_{\rm bf} + \Delta Q_{\rm coll} + \Delta Q_{\rm ff}$$ 0 = (18)

where $\Xi_{ij} = n_i C_{ij} - n_j C_{ji}$ and $\Gamma_{ij} = n_i R_{ij} - n_j R_{ji}$ are the collisional and radiative net rates respectively, $\Delta Q_{\rm bf}$, $\Delta Q_{\rm ff}$ and $\Delta Q_{\rm coll}$ are the differences between heating and cooling due to bound-free, free-free and collisional transitions, $\chi_{\rm tot}$ is the total opacity at a given frequency $\nu$ and $\eta_{\rm therm}$ is the total thermal emissivity at the same frequency. As we mentioned above, this system is highly nonlinear and contains a huge number of equations. Below, we use the ALO to solve Eqs. (15) and (16). The result is substituted in (17) and (18) which converts the initial system to a set of equations in only the unknown populations and the temperature.

   
3.1 Bound-free transitions

The analytical reduction of the system starts with Eqs. (15) and (16). Equation (15) is written in the ALO form

$$J_\nu = \Lambda^* S_\nu + (\Lambda-\Lambda^*) S_\nu^\dagger$$ (19)

and then $S_\nu$ from Eq. (16) is substituted in it. As a result we obtain an explicit dependence of $J_\nu$on the emission and the opacity.

$$J_\nu = \left(\Lambda^* \frac{\eta_{\rm therm}}{\chi_{\rm tot... ...) \left(1 - \frac{\sigma}{\chi_{\rm tot}}\Lambda^*\right)^{-1}$$ (20)

where $\Delta J_\nu^\dagger = (\Lambda-\Lambda^*) S_\nu^\dagger$ is calculated using the "old'' source function $S_\nu^\dagger$, calculated in the previous iteration. Following Rybicki & Hummer (1991) we also use $\chi_{\rm tot}$ from the previous iteration. The source function $S_\nu$ depends on the scattering and therefore on the new $J_\nu$. We iterate the Eqs. (16) and (20), with fixed opacity $\chi_{\rm tot}$ and emissivity $\eta_{\rm therm}$ until $J_\nu$ and $S_\nu$ are consistent with the current scattering opacity $\sigma$ and then use that solution in the transition rates $\Gamma_{il}$(see below).

The thermal emissivity in Eq. (20) can be written as the sum of the bound-free part $\eta_{\rm bf}$ and free-free part $\eta_{\rm ff}$which are proportional to the level population as

 

$$\eta_{\rm bf}(\nu) \sum_k n_k A_k(\nu)$$ =

$$\eta_{\rm ff}(\nu) \sum_m n_m B_m(\nu)$$ = (21)

where n_k are the populations of the upper levels of all transitions which emit at frequency $\nu$, n_m are the populations of the levels included in the free-free transitions, A_k and B_m collect all the terms which depend of the frequency and temperature, but do not depend on the level populations (Appendix A). Substituting (21) in (20) we obtain for $J_\nu$

 

$$J_\nu \left(\frac{\Lambda^*}{\chi_{\rm tot}} \left(\sum_k n_k A_k + \su... ...ta J^\dagger\right) \left(1-\frac{\sigma}{\chi_{\rm tot}} \Lambda^*\right)^{-1}$$ =

$$\tilde{\Lambda}^* \left(\sum_k n_k A_k + \sum_m n_m B_m\right) + \Delta \tilde{J}^\dagger$$ = (22)

which explicitly depends on the level populations and the radiation field obtained at the previous iteration. Here we define $\tilde{\Lambda}^*$ and $\Delta \tilde{J}^\dagger$as

 

$$\tilde{\Lambda}^* \frac{\Lambda^*}{\chi_{\rm tot}} \left(1-\frac{\sigma}{\chi_{\rm tot}} \Lambda^*\right)^{-1}$$ =

$$\Delta \tilde{J}^\dagger \Delta J^\dagger \left(1-\frac{\sigma}{\chi_{\rm tot}} \Lambda^*\right)^{-1} .$$ = (23)

The bound-free transition net rate between the bound level i and the target level l is

 

$$\Gamma_{il} n_i R_{il} - n_l R_{li} = n_i 4\pi \int_{\nu_0}^\infty \frac{\alpha_{il}}{h\nu} J_{\nu} {\rm d}\nu$$ =

$$- n_l \left(\frac{n_i^*}{n_l^*}\right) 4\pi \int_{\nu_0}^\infty \... ...left(\frac{2 h \nu^3}{c^2}+ J_\nu\right) {\rm e}^{-\frac{h\nu}{kT}} {\rm d}\nu.$$ (24)

Substituting (22) into (24) and taking into account (23) we obtain a form of the bound-free transition rates which depend only on the unknown level populations and are coupled to the radiation field through $\Delta \tilde{J}^\dagger$

 

$$\Gamma_{il} n_i ~ 4\pi \int_{\nu_{0_{il}}}^\infty \frac{\alpha_{il}}{h\nu} ~ \Delta \tilde{J}^\dagger_\nu ~ {\rm d}\nu$$ =

$$- n_l \left(\frac{n_i^*}{n_l^*}\right) ~ 4\pi ~ \int_{\nu_{0_il}}... ...2} + \Delta \tilde{J}^\dagger_\nu\right)~ {\rm e}^{-\frac{h\nu}{kT}}~{\rm d}\nu$$

$$+ \sum_k n_k ~ 4\pi ~ \int_{\nu_{0_{il}}}^\infty \frac{\alpha_{il... ...left(\frac{n_i^*}{n_l^*}\right)~ {\rm e}^{-\frac{h\nu}{kT}}\right] ~ {\rm d}\nu$$

$$+ \sum_m n_m ~ 4\pi ~ \int_{\nu_{0_{il}}}^\infty \frac{\alpha_{il... ...ft(\frac{n_i^*}{n_l^*}\right)~ {\rm e}^{-\frac{h\nu}{kT}}\right] ~ {\rm d}\nu .$$ (25)

In this expression, the sum over k contains all levels of all ions included in bound-free transitions. Therefore the level l is included in this sum also. For the rest of the levels, the product $\alpha_{il} A_k \sim \alpha_{il}~ \alpha_{nk}$(see Eq. (A.2) for the definition of A_k) decreases with increasing difference of the threshold frequencies of the transitions $i \rightarrow l$ and $n \rightarrow k$. Thus, the product is a measure of the overlapping of the continua and based on it, one can choose which of the transitions to take into account in the final system.

3.2 Bound-bound transitions

Our solution of the ERT is good for a static atmosphere. This is an acceptable approximation for the continuum transitions in the wind, but not for the lines. The generalization of the short characteristics method and the construction of $\Lambda^*$ for lines in non-hydrostatic media is not trivial and has been treated in another paper of this series (Zsargó et al. 2006). An alternative is the generalization of the Sobolev approximation for 3D flows (Rybicki & Hummer 1978). Its disadvantage is that, with the present technique, it does not easily allow the influence of the lines on the continuous radiation field (i.e., line blanketing) to be taken into account.

Similar to the ALO method described above, in a fast expanding wind, the Sobolev approximation allows the equations of radiative transport and the one describing the source function to be solved analytically.

If the Doppler shift between any two points in the wind displaces the line profile more than its thermal width, the points become radiatively decoupled and the equation of radiative transport becomes

$$\bar{J_{ij}} = (1-\beta_{ij}) S_{ij} + \langle \beta_c I_\nu \rangle$$ (26)

where $\bar{J_{ij}} = \int_0^{2\pi} {\rm d}\varphi \int_{-1}^1 {\rm d}\mu \int I_\nu \phi(\nu) {\rm d}\nu$ is the intensity averaged over the angles and the line profile $\phi(\nu)$. The Sobolev escape probability $\beta_{ij}$ is

$$\beta_{ij} = \frac{1}{4\pi} \int_0^{2\pi} \int_{-1}^1 \frac{1... ...\tau_{ij}/Q\left(\mu,\varphi\right)} {\rm d}\mu {\rm d}\varphi$$ (27)

and $\langle \beta_c I_\nu \rangle$ accounts for the interaction of the continuum radiation with the transition at a given point:

$$\langle \beta_c I_\nu \rangle = \frac{1}{4\pi} \int_0^{2\pi} ... ...tau_{ij}/Q\left(\mu,\varphi\right)} {\rm d}\mu {\rm d}\varphi.$$ (28)

The optical dept of the transition is

$$\tau_{ij} = \frac{\chi_{ij} c}{ \nu_{ij} % }$$ (29)

and $Q\left(\mu,\varphi\right)$ is

$$Q\left(\mu,\varphi\right) =\frac{1}{2} \sum n_k n_l \left(\fr... ...}{\partial r_l}+\frac{\partial v_l}{\partial r_k} \right)\cdot$$ (30)

Here n_k is the kth component of the vector $\vec{n}$ and r_l is the lth component of the radius vector $\vec{r}$ of the point of interest (Rybicki & Hummer 1978). The main difference between 1D and the 3D case is the way Q is calculated. In our code we calculate Q using all three components $v_r,v_\beta, v_\gamma$ of the velocity field and their nine derivatives.

The source function S_ij of the transition $i \leftrightarrow j$ can be written as

$$S_{ij} = \frac{n_j A_{ji}} { n_i B_{ij} - n_j B_{ji}}$$ (31)

where A_ji, B_ij are the Einstein coefficients. Substituting Eq. (26), we can write the bound-bound radiation rates as

 

$$\Gamma_{ij} n_i R_{ij} - n_j R_{ji} = \left( n_i B_{ij} - n_j B_{ji} \right) \bar{J} - n_j A_{ji}$$ = (32)

$$\left( n_i B_{ij} - n_j B_{ji} \right) \frac{n_j A_{ji}}{\left( n_i B_{ij} - n_j B_{ji} \right) } (1-\beta_{ij})$$ =

$$+ \left( n_i B_{ij} - n_j B_{ji} \right) \langle \beta_c I_\nu \rangle - n_j A_{ji}$$ (33)

$$\left( \left( n_i B_{ij} - n_j B_{ji} \right) \frac{\langle \beta_c I_\nu\rangle }{\beta_{ij}} - n_j A_{ji} \right) \beta_{ij},$$ = (34)

when $\tau_{ij} \rightarrow \infty$, $\beta \rightarrow 0$ the hence the net rate $\rightarrow 0.0$ assuring detailed balance in the lines.

The idea of the Sobolev approximation is to substitute the solution of the ERT into the rate equation. This is the same idea as in the ALO method. So, one can look at the Sobolev approximation as a ALO with $\Lambda^* = (1-\beta_{ij})$.

3.3 Equation of electron thermal balance

The next step is the formulation of the equation of electron thermal balance (18) as a function of the populations. We start with the bound-free transitions. Above, we defined the radiative net rates $\Gamma_{il}$ as integrals over frequency. Let us now define $\Gamma_{il}(\nu)$to be the frequency dependent part of Eq. (25). Then $\Gamma_{il} = \int_{\nu_{il}}^\infty \Gamma_{il}(\nu) {\rm d}\nu$. In this notation, $\Gamma_{il}(\nu)~(h\nu-h\nu_{il})$ is the amount of energy added or subtracted to the electron kinetic energy due to the transition $i \leftrightarrow l$. Then ETB due to all bound-free transitions is

$$\Delta Q_{\rm bf} = \sum_{i,l} \int_{\nu_{il}}^\infty \Gamma_{il}(\nu)~\left(h\nu-h\nu_{il}\right) {\rm d}\nu .$$ (35)

The substitution of the net rates from (25) convert $\Delta Q_{\rm bf}$ into a function of the unknown level populations only.

The electron thermal balance due to free-free transitions is

 

$$\Delta Q_{\rm ff} = 4 \pi N_{\rm e} \sum_m n_m \int_0^\infty \alp... ...kT}}\right)- \frac{2 h \nu^3}{c^2}{\rm e}^{-\frac{h\nu}{kT}} \right]{\rm d}\nu,$$ (36)

and the contribution of the collisional transitions is

$$\Delta Q_{\rm coll} = \sum \Xi_{ij}~\varepsilon_{ij}$$ (37)

where both bound-bound and bound-free transitions are taken into account and $\varepsilon_{ij}$ is the transition energy.

Following the idea described in Sect. 3.1, we substitute the radiation field from (22) in (35) and (36) (see Appendix A for details). The final equation is again a function only of the unknown level populations and the temperature.

3.4 Removing the nonlinearity from the final system

Substituting (25), (35), (36) and (37) in (17) and (18) at each point of the wind we obtain a system of N_l+1 equations for the N_l unknown level populations and the unknown temperature. The system is nonlinear on the level populations due to the products n_k n_i, n_k n_l in Eqs. (25), (35) and (36). Following Rybicki & Hummer (1991) and de Koter et al. (1993) we remove the nonlinearity by preconditioning. We assume that the level populations and the opacity $\chi_{\rm tot}$ do not vary rapidly from iteration to iteration. The quantities in the square bracket in (25), (36) and the opacity $\chi_{\rm tot}$ can be evaluated using the populations from the previous iteration. The resulting system becomes linear in the level populations but it is still nonlinear in the temperature. Unfortunately, the nonlinearity of the temperature is mainly in the form $\exp(-\frac{h\nu}{kT})$and it is not as easy to precondition as the level populations. So we linearize the preconditioned system in the temperature. As a result, we obtain the system

 

$$-\sum_{j<i} \Gamma_{ji} + \sum_{j>i} \Gamma_{ij} - \sum_{j<i} \Xi_{ji} + \sum_{j>i} \Xi_{ij}$$

$$+\left[-\sum_{j<i} \frac{\partial }{\partial T}\Gamma_{ji} + \sum... ...}\Xi_{ji} + \sum_{j>i} \frac{\partial }{\partial T}\Xi_{ij}\right] \Delta T = 0$$

$$\Delta Q_{\rm bf} + \Delta Q_{\rm ff} + \Delta Q_{\rm coll}$$

$$+\left[\frac{\partial }{\partial T}\Delta Q_{\rm bf} + \frac{\par... ...{\rm ff} + \frac{\partial }{\partial T}\Delta Q_{\rm coll}\right] \Delta T = 0.$$ (38)

Here all $\Gamma$'s and Q's are preconditioned and all temperature derivatives are calculated using the old populations. Appendix A summarizes the details of the derivatives. This final system is linear on all the variables and allows iterative solution.

3.5 Initial conditions and convergence

To facilitate the solution, different starting models can be used, and different techniques can be used to assist in the convergence. For the starting solution we can use a previously converged model, or use the results from a 1D CMFGEN (Hillier & Miller 1998) model. Since we use preconditioning, and not linearization, it not necessary for the starting solution to be "close'' to the final solution. Thus we can even start the model using a LTE initial guess which adds just a few iterations to the model.

To facilitate convergence we initially hold the temperature constant, until the changes of level populations are reduced below some given value - usually 100%. The temperature corrections are continued until they are reduced below some specified level, usually 1%. The code then continues to iterate on the populations until the changes are reduced below some specified limit.

Ng acceleration on the level populations is started when the changes drop below some given level, usually 20-30%. Starting the acceleration earlier leads to poorer convergence, a feature noted previously by Hillier (1990) and Rybicki & Hummer (1991). Ng acceleration of the radiation field, with strong electron scattering, is applied all the time if more than 25 iterations are needed.

In some cases the convergence rate is reduced due to stabilization or an oscillation in the corrections of a single level at one or two depth points. Figure 2 shows the convergence of a 1D model (Model 1, see Sect. 4.3) compared to the same model calculated with CMFGEN. No acceleration was applied. After the first 25 iterations, the only level whose population is changing is He  I 1s 2p ¹P$^\circ$ at a single depth. During the remaining iterations changes in this level governs the rate of convergence. This behavior is often encountered in other models, and is also seen in some CMFGEN models (particularly those using a diagonal operator and/or models in which not all levels are treated in the linearization). We are working on an algorithm which detects and avoids these kinds of problems. Our implementation of the preconditioning converges somewhat slower than the diagonal linearization scheme, but the difference is not significant and cannot be considered to be a major disadvantage.

	Figure 2: Convergence rate of Model 1 compared to the same model calculated with CMFGEN. The two codes were run from the same initial guess. In CMFGEN a diagonal linearization scheme was used. After iteration 25 of Model 1, the maximal relative correction to the level populations is determined by the He  I 1s 2p 1P$^\circ$ level.
Open with DEXTER

4 Test of the entire model

In this section we present four types of tests. First we tested the performance of the short characteristics method and the ALO in one and two dimensions using the tough case of gray opacity and pure scattering. Then a complete 1D non-LTE model was run and the results were compared with the model, using the same parameters, computed by CMFGEN (Hillier & Miller 1998). Finally, we ran a full 2D model to test the performance of the code in 2D. We also use this model to discuss the effects of the deviation from the spherical symmetry. The parameters of the non-LTE models presented below are specified as described in Table 1. The models with gray opacity are described in the text.

Table 1: Parameters of the test models.

4.1 Test of the solution of ERT. 1D Gray atmospheres

Before we could advance to the solution of the SEE we first performed tests to check the formal solver and the construction and performance of the ALO. We constructed a spherically symmetric envelope with pure gray scattering opacity $\chi_{\rm tot}(r) =\sigma= C r^{-2}$ and no thermal emission ( $\eta_{\rm therm} = 0.0$, and thus J=S, which is equivalent to solving the radiative transfer problem for a gray opacity source under the assumption of radiative equilibrium). The constant C is chosen so the total radial optical depth $\tau_{\rm tot} = \int_{r_{\min}}^{r_{\max}} \chi(r) {\rm d}r$is equal to a given value $T_{\rm tot}$ and we set the total flux H to a fixed value. The solution is started with $S = \rm const.$ and the code runs until the fractional corrections to the source function become less than $10^{-5}\%$. We ran the same test using the long characteristic code (Hillier 1996) and also compared the result with the semi analytical formula derived by Hummer & Rybicki (1971). Figure 3 shows this comparison for $T_{\rm tot}=15.0$and H=1.0 for models with a different number of radial grid points. The units of the flux H are arbitrary because the source function and the radiation field are simply proportional to this value. In all cases, the J as a function of r is close to the other solutions. The higher number of grid points leads to shorter characteristics. Consequently, the optical depth along them is calculated with better precision, the angle integration is more accurate, and the parabolic approximation of the source function S is better. As a result we obtain a more accurate solution as shown in Fig. 3.

Figure 3: Comparison between the average intensity in a spherically symmetric wind with gray pure scattering opacity $\chi \sim r^{-2}$ calculated with our short characteristics code and the long characteristics code (Hillier 1996). Solid line: 60 depth points; dashed line: 120 points; dash-dotted line: 200 points. The difference between the two solutions decreases with increasing number of points as expected if the main source of errors is the calculation of the optical depth along the short characteristics and the angle integration.

Open with DEXTER

Another indicator of the performance of the code is the conservation of the total flux. In a spherically symmetric atmosphere in radiative equilibrium, the quantity r² H(r), where H(r) is the total flux at radius r, is constant. In Fig. 4, the quantity r² H(r) is compared to its exact value H₀   R_*² for models with different numbers of points in the grid. As for the average intensity, the error in the flux is reduced when the number of points in the grid is increased, which again reflects the improvement of the approximations used to calculate the radiation field.

We conclude that our implementation of the short characteristic method performs well in spherical symmetry and we move to 2D test cases.

	Figure 4: Conservation of the flux for spherical models with a different number of radial points. Solid line - 100 depth points, dashed line - 200 points, dash-dotted line - 300 points. The ordinate is the difference between the calculated and the exact flux multiplied by r2. As expected, the difference between the calculated and exact fluxes decreases with increasing numbers of points.
Open with DEXTER

   
4.2 2D Gray atmosphere

The density structure of a spherically symmetric wind is defined by its mass-loss rate and the velocity field by means of the equation of continuity. But if the wind deviates from the symmetry, its density structure is not so simply obtained. In this paper we concentrate on a method to model axisymmetric winds, but do not discuss the origin of the asymmetry. We assume the wind to have a simple sine variation of the mass-loss rate $\dot M(\beta)$with the polar distance $\beta$.

$$\dot M(\beta) = {\dot M \over M_{\rm contr}} \left( 1 + \left[M_{\rm contr} -1 \right]\sin^{\alpha_1}\beta\right)$$ (39)

where $M_{\rm contr} = \dot M_{\rm equator}/\dot M_{\rm polar}$ is the contrast between the equatorial and polar mass-loss rate and $\dot M$ is the mass-loss rate which correspond to a spherically symmetric wind. The coefficient $\alpha_1$ governs the concentration of the gas towards the equator. Higher values of this parameter correspond to a more disk-like structure, while $\alpha_1=0$ corresponds to a spherically symmetric wind. Figure 6 shows the iso surfaces for $\tau=1$, which reflects the density structure, for two values of this parameter. The radial dependence of the density for each $\beta$ is maintained as in a spherically symmetric wind with the corresponding $\dot M(\beta)$.

As in the 1D case, we first compare the short characteristics solver to another short characteristics code, which was extensively tested (Busche & Hillier 2000). The same simple pure scattering opacity with r^-3 is used but with angular dependence of the mass-loss rate from Eq. (39) with $\dot M_{\rm contr}=2.0$ and concentration parameter $\alpha _1 = 2$. The left panel of Fig. 5 shows the average intensity calculated with the present short characteristics code divided by the J of a spherically symmetric wind with the same opacity as on the equator. The right panel shows the difference between the average intensity calculated with the two short characteristics codes. As it is seen the two codes give the same result within a few percent which, as above, we attribute mainly to the calculation of the optical depth and the interpolation techniques which are different in the two codes.

Figure 5: Left: the ratio of the average intensity in a wind with angular dependence of the opacity with contrast 2 and the spherically symmetric wind with the radial distribution of the opacity as on the equator. The pole is at $\beta =0$ (lower left angle). Right: comparison of the average intensity of the left panel with the same model calculated with the code described by Busche & Hillier (2000).

Open with DEXTER

The models discussed here have the polar wind optically thinner than the equatorial wind. As a result, the radial component of the flux at the pole is increased with respect to the spherical wind and the equatorial flux is reduced (Fig. 6, left). As the spherical symmetry is broken, the $\beta$ component of the flux (tangential flux) is not zero. The right panels of Fig. 6 show the variation of this component in the wind. The behavior of the tangential flux depends not only on the contrast between the polar and equatorial densities, but also on the actual density distribution. Right panels of Fig. 6 show this component of the flux for two models. They have the same mass-loss rate contrast $\dot M_{\rm contr}=100$ but a different concentration of the material towards the equator (see Eq. (39) below). The iso surfaces (for $\tau = 1.0$) of the two models are shown at the upper right corner of the right panels of Fig. 6. As it is seen, the more toroidal star (Fig. 6, upper panels) has tangential flux from the pole toward the equator (positive flux) at all radial distances and polar angles. The more disk-like star (Fig. 6, lower panels) has flux which changes its direction at different points in the wind. In general, one can attribute a force due to the radiation pressure for each component of the flux. The presence of $H_\beta$, which changes for different distributions of the material, might affect the hydrodynamics of the wind.

	Figure 6: Left: radial flux normalized to the same value for spherical wind. Right: tangential flux. The two models have the same contrast $\dot M_{\rm contr}=100.0$ but $\alpha _1 = 2$ for the upper panels and $\alpha _1 = 8$ for the lower. The figures at the upper right corner of the right panels show the iso surface of $\tau = 1.0$ for the two models. See (39) for definition of the parameters.
Open with DEXTER

   
4.3 1D non-LTE model test

Once we tested the solver of the radiation transport equation and the construction of the $\Lambda^*$ we were able to test the entire code. As above we first calculated 1D models and compared the results with a model with the same parameters but calculated with CMFGEN. The test was performed in two steps. First we checked the solution of the statistical equilibrium equation with a fixed temperature and then we checked the temperature correction procedure.

For the first test we chose the velocity structure as described in Hillier & Miller (1998) with parameters described in Table 1, Model 0 and the number of radial points N_r=60 . The solution was started with the temperature distribution of the converged CMFGEN model with the same parameters, but the starting level population and electron density was set to LTE. We calculated the same model with CMFGEN setting the code to use the Sobolev approximation and atoms with the same number of levels. The departure coefficients of the first five H I levels are shown in Fig. 7 together with the results obtained with CMFGEN. The two codes give results which differ by less than 3%. Again, we checked the behavior of the difference by increasing the number of grid points. The departure coefficient for the first level of H I for the model with N_r=120 radial points is shown in Fig. 7 with diamonds. As one might expect, the difference between the two models decreases with increasing N_r.

	Figure 7: Comparison between departure coefficients of the first five hydrogen levels for Model 0 (Table 1) calculated with CMFGEN (dashed lines) and the new code (solid lines). The curves for different levels are marked. The diamonds show the departure coefficient for the first H I level for a model with the same parameters but with twice as many radial points.
Open with DEXTER

The same model, but with added carbon levels (Model 1 in Table 1), was used to test the temperature correction procedure (Fig. 8). As in the case of the departure coefficients, the temperature distribution obtained with the new code was in good agreement with the temperature obtained with CMFGEN. As expected the emitted spectrum calculated with the new code's level populations also agrees well with the one calculated with CMFGEN. Figure 9 shows some of the lines of HI, HeI, HeII and CIV usually used for wind diagnostics. The spectra calculated by the two codes differ by less than 1%. For the test we switched off the dielectronic recombination and level dissolution which are treated in CMFGEN but are not yet implemented into the new code.

	Figure 8: Comparison between temperature distribution of Model 1 (Table 1) calculated with CMFGEN (dashed line) and the new code (solid line). The maximal difference is less than 3%.
Open with DEXTER

	Figure 9: Comparison between the emitted spectrum of Model 1 (Table 1) calculated with CMFGEN (thick gray line) and the new code (thin black line). The maximal difference is less than 1%.
Open with DEXTER

4.4 2D non-LTE tests

We introduced axial symmetry into the models by two means. First we modified the mass-loss rate at each polar distance as described in 4.2, Eq. (39). In the same way we assume an axisymmetric terminal velocity $V_\infty$ but we keep the flow radial, with the radial velocity obeying the standard velocity law:

 

$$V_r(r,\beta) V_\infty(\beta) (1-R_0(\beta)/r)^\omega$$ = (40)

$$V_\beta(r,\beta)$$ = 0 (41)

$$V_\phi(r,\beta)$$ = 0 (42)

$$V_\infty(\beta) V_\infty(\max) \left(V_{\rm contr}+(1-V_{\rm contr})\cos^{\alpha_2}\beta\right).$$ = (43)

As in Eq. (39), $V_{\rm contr} = V_{\rm equator}/V_{\rm pole}$ and the coefficient $\alpha_2$ together with $\alpha_1$ defines the thickness of the equatorial disk. Higher values correspond to a more defined disk. We use $\omega$ for the power index instead of the standard notation $\beta$ so there is no confusion with the polar distance coordinate. The code can handle deviations from spherical symmetry due to both mass-loss rate and velocity field, but in this first test we kept the velocity law symmetric ( $V_{\rm contr}=1$) and let only the mass-loss rate deviate from symmetry (Table 1). The continuity law, in that case, still has the usual solution

$$\rho(r,\beta) = \frac{\dot M(\beta)}{r^2 V_r(r,\beta)}\cdot$$ (44)

Figure 10 shows the temperature structure of Model2. The temperature close to the stellar core depends mainly on the optical depth and therefore the mass-loss rate. The higher the mass-loss rate, the higher the central temperature. This effect is seen in Fig. 10. The equatorial regions ( $\beta=90^\circ$) are hotter than the pole. The temperature of the outer part of the wind does not depend so strongly on the density but does depend on the number of coolers. In our simple case, the main coolant in the wind is C IV. In the polar regions carbon is ionized to C V, while in the equatorial regions it is mainly in C IV (Fig. 12). As a result, the wind temperature rises toward the pole.

	Figure 10: Temperature structure of 2D wind. The rapid rise of the temperature toward the pole is caused by the ionization of C IV which reduces the cooling.
Open with DEXTER

	Figure 11: Helium ionization. Left: He II; Right He III.
Open with DEXTER

The helium ionization balance shows a similar pattern. As often occurs in an O star's wind, He III recombines to He II in the region where the velocity increases to supersonic values. Subsequently, the ionization increases and the wind is maintained completely ionized. The more pronounced dip in the helium ionization creates a more pronounced rise in the temperature. Again, a strong ionization front in the tangential direction is formed in the outer part of the wind (Fig. 11). As we said above, we do not attempt to model a real wind in this paper. However even these simple calculations show that behavior of the ionization points to a possible correlation between the line strength, the line profile and the ionization stage, which could be used as a diagnostic tool.

	Figure 12: Carbon ionization. Left: C IV, Right: C V. The increase in the ionization of C IV near the pole ($\beta =0$) causes the rise in the temperature (as seen in Fig. 10).
Open with DEXTER

5 Conclusions

In this paper we have described a code we have developed to model non-spherically symmetric stellar winds. The code is based on a short characteristics solver for the radiation transport, preconditioning of the equations of statistical equilibrium combined with linearization of the electron thermal balance equation. Extensive testing of the code has been undertaken in order to verify its correctness. Importantly, the tests show that the short characteristics method, combined with ALO, preconditioning of the statistical equilibrium equations and the temperature correction based on the electron thermal balance, can model 2D winds.

The main limitation of the current version of the code is the Sobolev approximation which we will improve with the implementation of the new solver described in Zsargó et al. (2006). However we have designed the code with flexibility in mind, and it is relatively straightforward to insert alternative radiation field solvers. As expected, the increase in the number of dimensions increases the time per iteration. Fortunately, the short characteristics method allows for relatively easy parallelization which will also be implemented in the next version.

Even a simple model of an axisymmetric wind, as presented here, shows important differences in its temperature and ionization structure compared with a spherically symmetric wind. Such differences are expected, and are seen in other calculations (e.g., Kraus & Lamers 2003). The differences will cause changes in the emitted spectrum which in turn will provide observable diagnostics for the deviation from symmetry. Computations with the newly developed code will also be important for investigating the importance of the influence of the tangential flux on the hydrodynamics of stellar winds. Calculations show, for example, that the tangential flux can suppress the formation of wind compressed disks (Owocki et al. 1996).

The code, and implemented methods, are very flexible and allow the application of the code to many different objects with complex velocity fields such as accretion disks, rotating winds, interacting winds in binary systems and asymmetric outbursts in novae or supernovae.

Acknowledgements

This work was partially supported by Mexico CONACyT grants 40864-F and 42809 (L.G.). Support from NSF grants AST-9987390 and AST-0507328 is gratefully acknowledged.

   
Appendix A: Definition of the thermal processes and their thermal derivatives

We have previously shown that the net rate in a bound-free transition can be written as

 

$$\Gamma_{il} n_i 4\pi \int_{\nu_{0_{il}}}^\infty \frac{\alpha_{il}}{h\nu} ~ \Delta \tilde{J}^\dagger_\nu ~ {\rm d}\nu$$ =

$$- n_l \left(\frac{n_i^*}{n_l^*}\right) ~ 4\pi ~ \int_{\nu_{0_il}}... ...2} + \Delta \tilde{J}^\dagger_\nu\right)~ {\rm e}^{-\frac{h\nu}{kT}}~{\rm d}\nu$$

$$+ \sum_k n_k ~ 4\pi ~ \int_{\nu_{0_{il}}}^\infty \frac{\alpha_{il... ...left(\frac{n_i^*}{n_l^*}\right)~ {\rm e}^{-\frac{h\nu}{kT}}\right] ~ {\rm d}\nu$$

$$+ \sum_m n_m ~ 4\pi ~ \int_{\nu_{0_{il}}}^\infty \frac{\alpha_{il... ...left(\frac{n_i^*}{n_l^*}\right)~ {\rm e}^{-\frac{h\nu}{kT}}\right] ~ {\rm d}\nu$$ (A.1)

where

   

$$\frac{2 h \nu^3}{c^2} ~ C_I ~ N_{\rm e} ~ T^{-\frac{3}{2}} \frac{1}{g_k} \sum_{i_k} \alpha_{ik} g_i \exp\left(\frac{E_{ik}-h\nu}{kT}\right)$$                            A_k = (A.2)

$$\frac{2 h \nu^3}{c^2} N_{\rm e} \alpha_m \exp\left(-\frac{h \nu}{kT}\right)$$ B_m = (A.3)

$$\left[n_i-n_l\left(\frac{n_i^*}{n_l^*}\right){\rm e}^{-\frac{h\nu}{kT}}\right].$$ C_il = (A.4)

The definitions for A_k, B_m and C_il follow directly from the expressions given for the bound-free rates in (Mihalas 1978). The summation in A_k arises because a single ion (e.g. H⁺) can recombine to many levels. For simpler notation we rewrite these equations so that the terms $\frac{2 h \nu^3}{c^2}$ and $\exp(-\frac{h\nu}{kT})$ are separate from the sums:

$$\frac{2 h \nu^3}{c^2} a_k \exp\left(-\frac{h \nu}{kT}\right)$$ (A.5)

$$\frac{2 h \nu^3}{c^2} b_m \exp\left(-\frac{h \nu}{kT}\right)$$ (A.6)

and

$$\sum_i \left(\frac{n_i^*}{n_k^*}\right) \alpha_{ik}$$ (A.7)

$$N_{\rm e} \alpha_{m} .$$ (A.8)

Now we rewrite (A.1) in the form

$$\Gamma_{il} = n_i \Gamma_{il}^{1} + n_l \Gamma_{il}^{2} + \sum_k n_k\Gamma_{il}^{3k} + \sum_m n_m\Gamma_{il}^{4m} .$$ (A.9)

With these abbreviations we can rewrite the four terms of the net rate from (A.9) as

    

$$\Gamma_{il}^{1} 4\pi \int_{\nu_{0il}}^\infty \frac{\alpha_{il}}{h\nu} \Delta \tilde{J}^\dagger_\nu {\rm d}\nu$$ = (A.10)

$$\Gamma_{il}^{2} -\left(\frac{n_i^*}{n_l^*}\right) 4\pi \int_{\nu_{0_il}}^\infty \... ...c^2} + \Delta \tilde{J}^\dagger_\nu\right) {\rm e}^{-\frac{h\nu}{kT}}{\rm d}\nu$$ = (A.11)

$$\Gamma_{il}^{3k} 4\pi \int_{\nu_{0il}}^\infty \frac{\alpha_{il}}{h\nu} \frac{2 h \... ...Lambda}^* ~ a_k\left(\nu\right) ~ C_{il} {\rm e}^{-\frac{h \nu}{kT}} {\rm d}\nu$$ = (A.12)

$$\Gamma_{il}^{4m} 4\pi \int_{\nu_{0il}}^\infty \frac{\alpha_{il}}{h\nu} \frac{2 h \... ...mbda}^* ~ b_m\left(\nu\right) ~ C_{il} {\rm e}^{-\frac{h \nu}{kT}} {\rm d}\nu .$$ = (A.13)

The two terms $\Gamma_{il}^{3k}$ and $\Gamma_{il}^{4m}$ contains the upper level l as well. We can combine these two terms with $\Gamma_{il}^{2}$ and write the recombination rate as

 

$$\Gamma_{il}^{2}+\Gamma_{il}^{3l}+\Gamma_{il}^{4l} = -4\pi \int_{\... ...t] + \Delta \tilde{J}^\dagger_\nu\right) {\rm e}^{-\frac{h\nu}{kT}} {\rm d}\nu.$$ (A.14)

Now we have all the necessary evaluations for the net transition rates. For the temperature correction we also need their derivatives. First we calculate the derivative of A_k and B_m.

Having in mind that

$$\frac{\partial}{\partial T} {\rm e}^{-\frac{h\nu}{kT}} = \frac{h\nu}{kT^2}{\rm e}^{-\frac{h\nu}{kT}}$$ (A.15)

and

$$\frac{\partial}{\partial T} \left(\frac{n_i^*}{n_l^*}\right) ... ...\right) \left(\frac{3}{2}+\frac{\varepsilon_{il}}{kT}\right)/T$$ (A.16)

we can write the derivative of the $\Gamma$'s as follow:

$$\frac{\partial}{\partial T}\Gamma_{il}^{1} = 0$$ (A.17)

 

$$\frac{\partial}{\partial T}\Gamma_{il}^{2} 4\pi \int_{\nu_{0il}}^\infty \frac{\alpha_{il}}{h\nu} \left(\frac{2 h \nu^3}{c^2} + \Delta \tilde{J}^\dagger_\nu \right)$$ =

$$\times \left[-\left(\frac{n_i^*}{n_l^*}\right) \left(\frac{3}{2}+... ...{il}}{kT}\right)/T+\frac{h\nu}{kT^2}\right]{\rm e}^{-\frac{h\nu}{kT}}{\rm d}\nu$$ (A.18)

$$\frac{\partial}{\partial T}\Gamma_{il}^{3k} 4\pi \int_{\nu_{0il}}^\infty \frac{\alpha_{il}}{h\nu} \frac{2 h \nu^3}{c^2} \tilde{J}^\dagger_\nu$$ =

$$\times \left[\dot a_k C_{il}+a_k \dot C_{il}+a_k C_{il} \frac{h\nu}{kT^2}\right]{\rm e}^{-\frac{h\nu}{kT}}{\rm d}\nu$$ (A.19)

$$\frac{\partial}{\partial T}\Gamma_{il}^{4m} 4\pi \int_{\nu_{0il}}^\infty \frac{\alpha_{il}}{h\nu} \frac{2 h \nu^3}{c^2} \tilde{J}^\dagger_\nu$$ =

$$\times \left[\dot b_m C_{il}+b_m \dot C_{il}+b_m C_{il} \frac{h\nu}{kT^2}\right]{\rm e}^{-\frac{h\nu}{kT}}{\rm d}\nu$$ (A.20)

where the $\dot a$, $\dot b$ and $\dot C$ are the temperature derivatives

$$\begin{eqnarray*}\dot a_k = \frac{\partial a_k}{\partial T} &=& \sum_i \alpha_{i... ... \left(-\frac{3}{2}+\frac{h\nu -\varepsilon_{il}}{kT}\right)/T . \end{eqnarray*}$$

The collisional transition rates $\Xi$ were defined as

$$\Xi_{ij} = n_i C_{ij} - n_j C_{ji} = n_i N_{\rm e} q_{ij}\left(T\right) \left(1-\frac{n_j}{n_i} \frac{n_i^*}{n_l^*}\right)$$ (A.21)

so its temperature derivative is straight forward

 

$$\frac{\partial}{\partial T} \Xi_{ij} = n_i N_{\rm e} q_{ij}\left[... ... \ln{q_{ij}}}{\partial T}\right)\frac{n_j}{n_i} \frac{n_i^*}{n_l^*}\right]\cdot$$ (A.22)

The derivatives are calculated for the level population at the previous iteration so we assume n_i, n_j and $N_{\rm e}$ do not depend on the temperature.

Finally we deduce the heating-cooling balances. The balance between heating and cooling by bound-free transitions is more laborious. But, if we define $\Gamma_{il}\left(\nu\right)$ as in (35), the balance $\Delta Q_{\rm bf}$ becomes

$$\Delta Q_{\rm bf} = \sum_{i,l} \int_{\nu_{il}}^\infty \Gamma_{il}(\nu)~\left(h\nu-h\nu_{il}\right) {\rm d}\nu$$ (A.23)

and its temperature derivative is

$$\frac{\partial}{\partial T} \Delta Q_{\rm bf} = \sum_{i,l} \i... ...{il}(\nu)}{\partial T} ~\left(h\nu-h\nu_{il}\right) {\rm d}\nu$$ (A.24)

with all the quantities already defined in (A.1) and (A.20).

The same idea gives us the collisional thermal balance $\Delta Q_{\rm coll}$ and its derivative.

$$\Delta Q_{\rm coll} \sum_{i,j} \Xi_{ij} \varepsilon_{ij}$$ =

$$\frac{\partial}{\partial T} \Delta Q_{\rm coll} \sum_{i,j} \frac{\partial}{\partial T}\Xi_{ij} \varepsilon_{ij}$$ = (A.25)

where the quantities are defined in (A.21) and (A.22) and $\varepsilon_{ij}$ is the transition energy. This equation if valid for both bound-bound and bound-free transitions where the only difference is in the transition probability $q_{ij}\left(T\right)$ and its derivative.

The last term in (18) is the free-free transitions temperature balance, given by

$$\Delta Q_{\rm ff} = N_{\rm e} \sum_m n_m n_m 4\pi\int_0^\infty \chi_{\rm ff}^m (J^\dagger_\nu - B_\nu(T)) {\rm d}\nu$$ (A.26)

where $\chi_{\rm ff}^m = \alpha_{\rm ff}^m(1-{\rm e}^{-\frac{h\nu}{kT}})$ is the free-free transition cross section corrected for the stimulated emission. The temperature derivative is:

$$\frac{\partial }{\partial T}\Delta Q_{\rm ff}^m = N_{\rm e} n_m 4... ...nu(T)) -\chi_{\rm ff}^m \frac{\partial }{\partial T}B_\nu(T) \right] {\rm d}\nu$$ (A.27)

where and $\frac{\partial }{\partial T}\chi_{\rm ff}^m$ is

$$\frac{\partial }{\partial T}\chi_{\rm ff} \left(1-{\rm e}^{-\frac{h\nu}{kT}}\right) \alpha_{\rm ff} \left(\... ...{\alpha_{\rm ff}} + \frac{h\nu}{kT^2}\right) - \alpha_{\rm ff}\frac{h\nu}{kT^2}$$ = (A.28)

$$\chi_{\rm ff} \omega - \alpha_{\rm ff}\frac{h\nu}{kT^2} \cdot$$ = (A.29)

We treat, the free-free opacity and emissivity as well as electron scattering opacity as background processes in the solution of the radiative transfer Eq. (22) (Rybicki & Hummer 1991). Then we use $J^\dagger_\nu$ from the previous iteration to evaluate $\Delta Q_{\rm ff}$ and its derivative.

These equations define all necessary members of the system (38), and its solution give the unknown level populations n_i and the correction to the temperature $\Delta T$.

References

• Auer, L. H., & Koenigsberger, G. 1994, ApJ, 436, 859 [NASA ADS] [CrossRef]
• Auer, L. H., & Paletou, F. 1994, A&A, 285, 675 [NASA ADS]
• Bjorkman, J. E., & Cassinelli, J. P. 1993, ApJ, 409, 429 [NASA ADS] [CrossRef]
• Bjorkman, J. E., & Wood, K. 2001, ApJ, 554, 615 [NASA ADS] [CrossRef]
• Busche, J. R., & Hillier, D. J. 2000, ApJ, 531, 1071 [NASA ADS] [CrossRef]
• de Koter, A., Schmutz, W., & Lamers, H. J. G. L. M. 1993, A&A, 277, 561 [NASA ADS]
• Dullemond, C. P., & Turolla, R. 2000, A&A, 360, 1187 [NASA ADS]
• Hatchett, S., & McCray, R. 1977, ApJ, 211, 552 [NASA ADS] [CrossRef]
• Hill, G. M., Moffat, A. F. J., & St-Louis, N. 2002, MNRAS, 335, 1069 [NASA ADS] [CrossRef]
• Hillier, D. J. 1990, A&A, 231, 116 [NASA ADS]
• Hillier, D. J. 1996, A&A, 308, 521 [NASA ADS]
• Hillier, D. J. 2003, in ASP Conf. Ser., 199
• Hillier, D. J., & Miller, D. L. 1998, ApJ, 496, 407 [NASA ADS] [CrossRef]
• Hubeny, I., & Lanz, T. 1995, ApJ, 439, 875 [NASA ADS] [CrossRef]
• Hummer, D. G., & Rybicki, G. B. 1971, MNRAS, 152, 1 [NASA ADS]
• Juvela, M. 1997, A&A, 322, 943 [NASA ADS]
• Kraus, M., & Lamers, H. J. G. L. M. 2003, A&A, 405, 165 [EDP Sciences] [NASA ADS] [CrossRef] (In the text)
• Kubát, J., Puls, J., & Pauldrach, A. W. A. 1999, A&A, 341, 587 [NASA ADS]
• Kunasz, P., & Auer, L. H. 1988, J. Quant. Spectrosc. Radiat. Transf., 39, 67 [NASA ADS] [CrossRef]
• Lucy, L. B. 1999a, A&A, 344, 282 [NASA ADS]
• Lucy, L. B. 1999b, A&A, 345, 211 [NASA ADS]
• Lucy, L. B. 2003, A&A, 403, 261 [EDP Sciences] [NASA ADS] [CrossRef]
• Maeder, A. 1999, A&A, 347, 185 [NASA ADS]
• Maeder, A., & Meynet, G. 2003, in IAU Symp., 267
• Mihalas, D. 1978, Stellar atmospheres, 2nd edition (San Francisco: W. H. Freeman and Co.), 650
• Mihalas, D., Auer, L. H., & Mihalas, B. R. 1978, ApJ, 220, 1001 [NASA ADS] [CrossRef]
• Olson, G. L., Auer, L. H., & Buchler, J. L. 1986, J. Quant. Spectrosc. Radiat. Transf., 35, 431 [NASA ADS] [CrossRef]
• Owocki, S. P., Cranmer, S. R., & Blondin, J. M. 1994, ApJ, 424, 887 [NASA ADS] [CrossRef]
• Owocki, S. P., Cranmer, S. R., & Gayley, K. G. 1996, ApJ, 472, L115 [NASA ADS] [CrossRef] (In the text)
• Park, Y.-S., & Hong, S. S. 1998, ApJ, 494, 605 [NASA ADS] [CrossRef]
• Rybicki, G. B., & Hummer, D. G. 1978, ApJ, 219, 654 [NASA ADS] [CrossRef]
• Rybicki, G. B., & Hummer, D. G. 1991, A&A, 245, 171 [NASA ADS]
• van Noort, M., Hubeny, I., & Lanz, T. 2002, ApJ, 568, 1066 [NASA ADS] [CrossRef]
• Whitney, B. A., Wood, K., Bjorkman, J. E., & Cohen, M. 2003a, ApJ, 598, 1079 [NASA ADS] [CrossRef]
• Whitney, B. A., Wood, K., Bjorkman, J. E., & Wolff, M. J. 2003b, ApJ, 591, 1049 [NASA ADS] [CrossRef]
• Zsargó, J., Hillier, D. J., & Georgiev, L. N. 2006, A&A, 447, 1093 [EDP Sciences] [NASA ADS] [CrossRef]