2 Model atmospheres

The non-LTE spectrum formation is calculated for an expanding atmosphere under the standard assumptions of spherical symmetry, stationarity and homogeneity. The model calculations are in line with our previous work (Koesterke et al. 1992; Hamann et al. 1992; Koesterke & Hamann 1995; Leuenhagen & Hamann 1994; Leuenhagen et al. 1996; Hamann & Koesterke 1998). However, for the inclusion of iron group line-blanketing we had to modify our code extensively. In this section, we give a summarizing description of the method with special emphasize to those new features which concern the iron line-blanketing.

   
2.1 Model parameters and atomic data

A model atmosphere is specified by the luminosity and radius of the stellar core at the inner boundary, and by the chemical composition and the density- and velocity structure of the envelope. For the stellar core the radius $R_\star$ at Rosseland optical depth $\tau_{\rm R}=20$ and the stellar temperature $T_\star$ are prescribed. $T_\star$ is related to the stellar luminosity $L_\star$ by Stefan-Boltzmann's law

$$% T_\star = \left( \frac{L_\star}{4\pi~\sigma_{\rm SB}~R_\star^2} \right)^\frac{1}{4}\cdot$$ (1)

 

 

Table 1: Summary of the model atom. The iron group ions (Fe) are described by a relatively small number of super-levels, each representing a large number of true energy levels (called "sub-levels'', see Sect. 3.1).

Ion	Levels		Ion	Super-levels	Sub-levels
He I	17		Fe III	1
He II	16		Fe IV	18	30 122
He III	1		Fe V	19	19 804
C I	2		Fe VI	18	15 155
C II	32		Fe VII	16	11 867
C III	40		Fe VIII	10	8669
C IV	54		Fe IX	11	12 366
C V	1		Fe X	1
O II	3
O III	33
O IV	25
O V	36
O VI	15
O VII	1
Si III	10
Si IV	7
Si V	1

For the expansion of the envelope a velocity law of the form

$$% v(r)=v_\infty\;\left(1-\frac{r_0}{r}\right)$$ (2)

with the terminal velocity $v_\infty$ is adopted. This outer velocity law is augmented by an inner part, which corresponds to an exponentially decreasing density distribution in the nearly hydrostatic domain. r₀ is suitably determined to connect both domains smoothly. The density in the atmosphere is related to the velocity by the equation of continuity $\dot{M} = 4\pi~r^2~\rho(r)~v(r)$ and therefore requires the specification of the mass loss rate $\dot{M}$.

Density-inhomogeneities (clumping) are accounted for in the limit of small-scale clumps with a density enhanced by a factor D = 1/f_V over the mean density $\rho$(see Schmutz 1995; Hillier 1996; Hamann & Koesterke 1998). The inter-clump medium is supposed to be void. The radiation transport is calculated for the spatially averaged opacity, whereas the statistical equations are solved for the enhanced density in the clumps.

For models with the same stellar temperature $T_\star$, the strength of emission lines depends mainly on the so-called transformed radius $R_{\rm t}$, which is defined as

$$% R_{\rm t} = R_\star \left[\frac{v_\infty}{2500 ~ {\rm km}~{... ...10^{-4} ~ {M_\odot}~{\rm yr^{-1}}}\right]^{2/3} \right. \cdot$$ (3)

Models with the same $R_{\rm t}$ are found to show similar line equivalent widths (Schmutz et al. 1989). This invariance holds due to the fact that the line emission in WR stars is dominated by recombination processes (see Hamann & Koesterke 1998).

The chemical composition is given by mass fractions $X_{\rm He}$, $X_{\rm C}$, $X_{\rm O}$, $X_{\rm Si}$ and $X_{\rm Fe}$ of helium, carbon, oxygen, silicon and iron group elements. The model atoms contain the ionization stages He I-He III, C I-C V, O II-O VII, Si III-Si V and Fe III-Fe X. Except for He I and Si III, the lowest and highest ionization stages are restricted to a few auxiliary levels. A summarizing description of the model atoms is given in Table 1.

   
2.2 Radiation transfer

Due to the iron group elements, almost the whole spectral range is crowded by lines. Therefore, we abandon the distinction between continuum- and spectral line transfer used in our previous code. Analogous to the work of Hillier & Miller (1998), the radiation transfer is now calculated on one comprehensive frequency grid, which covers the whole relevant frequency range with typically three points per Doppler width wherever line opacities are present, and a wider spacing in pure-continuum regions.

The equation of radiative transfer in a spherically-expanding atmosphere is formulated in the co-moving frame of reference, neglecting aberration and advection terms (Mihalas et al. 1976a). The angle dependent transfer equation then becomes a partial differential equation for the intensity $I_\nu$,

 

$$% { \mu\frac{\partial I_\nu}{\partial r} + \frac{1-\mu^2}{r} ~ \f... ...} + \mu^2 \frac{{\rm d}v}{{\rm d}r} \right) \frac{\partial I_\nu}{\partial x} }$$

$$\eta_\nu - \kappa_\nu I_\nu = \eta_\nu^{\rm true} + \kappa^{\rm Th} J_\nu - (\kappa_\nu^{\rm true} + \kappa^{\rm Th}) ~ I_\nu,$$ = (4)

where $\mu$ is the angle cosine. The dimensionless frequency x is in Doppler units, while the velocity v and the velocity gradient $\frac{{\rm d}v}{{\rm d}r}$ must be expressed in units of the corresponding reference velocity. Thomson scattering is accounted for under the assumption of coherence and isotropy, i.e. the total opacity splits into a "true'' part and a "Thomson'' part $\kappa_\nu = \kappa_\nu^{\rm true} + \kappa^{\rm Th}$, and the total emissivity into a true part and the corresponding scattering part $\eta_\nu = \eta_\nu^{\rm true} + \kappa^{\rm Th} J_\nu$. Whereas the "true'' terms on the right-hand side of Eq. (4) can be calculated directly from the non-LTE population numbers, the angle-averaged intensity $J_\nu$ in the scattering term includes the intensities of different angles. Therefore, $J_\nu$ must either be taken from the preceding iteration, or the complete set of equations must be solved for all angles at once.

In the present work, the numerical solution of Eq. (4) is achieved by a short-characteristic method, described in detail in Koesterke et al. (2002). However, as the angle-dependent equation is computationally expensive, and because of the drawback with the Thomson-scattering term, we employ the "method of variable Eddington factors'' (Auer & Mihalas 1970). This means, the angle-dependent transfer equation is only solved from time to time (typically every six ALI iteration cycles). The resulting intensities $I_\nu$ at each radius and frequency grid point are integrated up over angles with weight factors $\mu^n$, yielding the nth moments of the radiation field J, H, K and N for n = 0, 1, 2 and 3, respectively. Only the Eddington factors f and g are finally exploited from that step, which are defined as

$$% f = \frac{K_\nu}{J_\nu} ~~~{\rm and}~~~ g = \frac {N_\nu} {H_\nu + \epsilon J_\nu}\cdot$$ (5)

The short-characteristic solution formalism guarantees always positive $I_\nu$ and thus positive $J_\nu$. However, for the flux-like moments $H_\nu$ and $N_\nu$ negative or zero values may occur. This is the reason to introduce the above definition of g as a "mixed'' Eddington factor (as originally proposed by Hillier & Miller 1998). The $\epsilon$ (which may be set to different values for different radius points) is adjusted to guarantee a positive denominator at all frequencies. Moreover, it is chosen large enough to ensure that $g~(\frac{v}{r} - \frac{{\rm d}v}{{\rm d}r}) < \frac{v}{r}$, which is mandatory for the system of moment equations (Eqs. (6) and (7)) being of hyperbolic type.

Based on these Eddington factors, the radiative transfer is solved in each ALI iteration cycle by means of the moment equations. By integration of Eq. (4) over $\mu$(weighted with $\mu^0$ and $\mu^1$ respectively), one obtains the zeroth and first moment equation as

$$% 0.\!:~\frac{\partial \tilde{H}_\nu}{\partial r} + \left( \f... ... = \kappa_\nu^{\rm true} ( \tilde{S}_\nu - \tilde{J}_\nu )\\$$ (6)

$$% 1.\!:~\frac{\partial (q\tilde{K}_\nu)}{-q\partial r} - \le... ...tial \tilde{H}_\nu}{\partial x} = \kappa_\nu ~\tilde{H}_\nu ,$$ (7)

where q is the "sphericality factor'' (Auer 1971, but in the form defined by Mihalas & Hummer 1974). The tilde indicates quantities multiplied by r². On the right hand side of Eq. (6) the scattering terms cancel out, leaving only the true source function $S_\nu = \eta_\nu^{\rm true} / \kappa_\nu^{\rm true}$.

After the moments $K_\nu$ and $N_\nu$ are substituted in Eqs. (6) and (7) with the help of the Eddington factors (Eq. (5)), these are solved by a differencing scheme as proposed by Mihalas et al. (1976b).

   
2.3 Statistical equations and accelerated lambda iteration

In the ALI formalism, the consistent solution of the equations of statistical equilibrium and the radiative transfer equation is achieved iteratively by solving both sets of equations in turn. However, in order to obtain convergence it is necessary to "accelerate'' the iteration by incorporating some "approximate'' radiative transfer into the statistical equation, which at least accounts for the locally trapped radiation in optically thick situations.

In the present work we use the concept of Hamann (1985, 1986). However, our definition of the diagonal "approximate lambda operators'' must be modified and extended with respect to the iron-line opacities, as we will describe in the following.

As in Hamann (1986), the atomic population numbers $\vec{n}$ are calculated at each depth point from the equations of statistical equilibrium, which are of the form

$$% \vec{n}\cdot {\cal P}(n_{\rm e}, \vec{{J}}(\vec{n}))=\vec{b}.$$ (8)

The rate matrix ${\cal P}$ depends on the electron density $n_{\rm e}$ and on the radiation field $\vec{{J}}(\vec{n})$. In the formal solution of the transfer equation (as described in Sect. 2.2), the radiation field $\vec{{J}}^{\rm fs}$ is obtained from the old populations $\vec{n}_{\rm old}$. However, as the essential point of the ALI technique, the $\vec{{J}}^{\rm fs}$ are modified by a correction term which pre-estimates the local response on the new population numbers $\vec{n}$ just being calculated. The resulting $\vec{{J}}(\vec{n})$ enters the coefficient matrix ${\cal P}$, which then becomes dependent on the total set of population numbers $\vec{n}$ itself, i.e. the system of rate equations becomes non-linear in $\vec{n}$.

This set of nonlinear equations is solved numerically by the application of a hybrid technique, a combination of the Broyden and the Newton algorithm (Hamann 1987; Koesterke et al. 1992). For this purpose the Jacobian matrix $\partial \cal{P}/\partial \vec{n}$ must be calculated occasionally, and the derivatives $\partial \vec{{J}}/\partial \vec{n}$ must be provided for all relevant transitions.

The calculation of the $\vec{{J}}(\vec{n})$ and their derivatives is modified in comparison to Hamann (1985). For line transitions we keep the core saturation approach but utilize the "diagonal operator'' of Rybicki & Hummer (1991) for the core integration, additionally accounting for the interaction with blending iron opacities. The continuum and iron transition rates are also calculated using the diagonal operator but in the "standard'' approach, i.e. the full frequency integral is performed (see Sect. 3.3 for iron transitions).

The diagonal operator provides the local response of the radiation field on the population numbers as a by-product of the radiative transfer. By applying it to the discretized moment Eqs. (6) and (7), we obtain the response on the true source function, fully accounting for coherent scattering.

In discretized form $J_\nu(r)$ is represented by a vector $\vec{J}_k$ on the radial depth grid at the corresponding frequency index k, and the first moment vector $\vec{H}_k$ is defined on the radius interstices. After elimination of $\vec{H}_k$ by substitution of Eq. (7) into Eq. (6) one obtains a second order equation of the form

$$% \vec{T}_k \vec{J}_k = \vec{V}_k \vec{H}_{k-1 } + \vec{U}_k \vec{J}_{k-1} + \vec{S}_k,$$ (9)

with the tridiagonal coefficient matrix $\vec{T}_k$, the diagonal matrix $\vec{U}_k$, the lower bidiagonal matrix $\vec{V}_k$, the vector of source functions $\vec{S}_k$, and the zero and first moment vectors  $\vec{J}_{k-1}$ and  $\vec{H}_{k-1}$ at the blueward frequency index k-1. This linear system of equations can be solved for $\vec{J}_k$ by multiplication with the inverse of the tridiagonal matrix  $\vec{T}_k^{-1}$.

At this point it is possible to extract the diagonal $\vec{D}_k \equiv {\rm diag}(\vec{T}_k^{-1})$, which describes the local response of the radiation field on the source function, i.e.

$$% \frac{\partial J_k}{\partial S_k} = D_k,$$ (10)

with the diagonal element D_k at the considered depth point. For the optically thick case, where the radiation transfer is in fact local, D_k approaches unity.

For a line transition between the energy levels l and u, the radiation field enters the radiative rates via the scattering integral

$$% {J}_{lu} = \int J_\nu~\phi_{lu}(\nu)~{\rm d}\nu,$$ (11)

with the normalized profile function $\phi_{lu}(\nu)$. For the calculation of ${J}_{lu}(\vec{n})$ we apply Eq. (10) only in the optically thick line core, and obtain

$$% {J}_{lu}(\vec{n}) {=} {J}_{lu}^{\rm fs} {+} \int_{\nu_{\rm... ...-} S_\nu(\vec{n}^{\rm old}) \right) \phi_{lu}(\nu)~{\rm d}\nu$$ (12)

with the core-confining frequencies $\nu_{\rm red}$ and $\nu_{\rm blue}$ as defined in Hamann (1985). Due to the complex interactions with blending lines, it is essential to calculate $S_\nu(\vec{n})$ in a consistent way to the radiation transport. Therefore, $S_\nu(\vec{n})$ must be calculated on the same frequency grid and on basis of the same opacities (including the iron group) as in the formal solution.

The derivatives $\partial {J}_{lu}/\partial \vec{n}$ are calculated directly from Eq. (12)

$$% \frac{\partial {J}_{lu}}{\partial \vec{n}} = \int_{\nu_{\r... ...{\partial S_\nu}{\partial \vec{n}}~ \phi_{lu}(\nu)~{\rm d}\nu.$$ (13)

In order to save numerical effort, only the derivatives with respect to the lower and upper levels of the considered line transition are calculated, while the dependence on background opacities and emissivities is neglected here.

The temperature stratification has to be derived from the equation of radiative equilibrium in the co-moving frame

$$% { \int_{0}^{\infty} \eta_\nu - \kappa_\nu J_\nu~{\rm d}\nu = 0. }$$ (14)

In previous work this equation was incorporated into the system of rate equations as an additional constraint for the electron temperature $T_{\rm e}$. This approach has the disadvantage that the whole set of rate equations is strongly coupled by the additional constraint. As the speed and accuracy of the numerical solution decrease for an increasing number of equations, this limits the total number of energy levels taken into account. In our new approach Eq. (14) is solved separately from the rate equations. The different chemical elements are then only weakly coupled via the radiation field, and the method is mainly limited by the number of energy levels per element. Equation (14) is solved by a temperature correction procedure, which is similar to the method of Unsöld (1955) and Lucy (1964) (see Mihalas 1978, p. 174). For the adaption of this method to the non-LTE case, mean values weighted by the Planck function are substituted by mean values weighted by the non-LTE source function. Additional means of the Eddington factors are introduced. After convergence, the resultant model flux (including the mechanical work which is performed by the radiation field, cf. Hillier & Miller 1998) is constant within $\pm 0.5\%$. The two terms in Eq. (14) then cancel to a relative accuracy of $\approx$10^-11 in the outer part of the wind and $\approx$10^-7 in the photosphere.

A sufficient model accuracy is reached when the correction between consecutive ALI iterations drops below 1%. The convergence properties depend strongly on the start model. In practice it is often possible to choose start models which are already close to the solution. In such a case about 40 iterations are needed to obtain convergence. For a model which is started from LTE, this number can increase by a factor of 10. In the case of the line blanketed WC model in Sect. 4 (70 depth points, 396 levels, 36 210 frequencies) 24 iterations are needed from a "good'' start model, whereas 300 iterations are necessary for an LTE start. The computing time per ALI iteration cycle ranges from 150 to 350 s on a Compaq Alpha XP1000/667 workstation.