A&A 435, 669-698 (2005)
DOI: 10.1051/0004-6361:20042365

Atmospheric NLTE-models for the spectroscopic analysis of blue stars with winds

II. Line-blanketed models

J. Puls¹ - M. A. Urbaneja² - R. Venero³ - T. Repolust¹ - U. Springmann⁴ - A. Jokuthy¹ - M. R. Mokiem⁵

1 - Universitäts-Sternwarte München, Scheinerstrasse 1, 81679 München, Germany
2 - Institute for Astronomy, University of Hawaii at Manoa, 2680 Woodlawn Drive, Honolulu, Hawaii 96822, USA
3 - Facultad de Ciencias Astronómicas y Geofísicas, Universidad Nacional de La Plata, Paseo del Bosque s/n, B1900FWA La Plata, Argentina
4 - BT (Germany) GmbH & Co. oHG, Barthstr. 22, 80339 München, Germany
5 - Astronomical Institute "Anton Pannekoek'', Kruislaan 403, 1098 SJ Amsterdam, The Netherlands

Received 15 November 2004 / Accepted 6 February 2005

Abstract
We present new or improved methods for calculating NLTE, line-blanketed model atmospheres for hot stars with winds (spectral types A to O), with particular emphasis on fast performance. These methods have been implemented into a previous, more simple version of the model atmosphere code F ASTWIND (Santolaya-Rey et al. 1997) and allow us to spectroscopically analyze large samples of massive stars in a reasonable time-scale, using state-of-the-art physics. Although this updated version of the code has already been used in a number of recent investigations, the corresponding methods have not been explained in detail so far, and no rigorous comparison with results from alternative codes has been performed. This paper intends to address both topics.

In particular, we describe our (partly approximate) approach to solve the equations of statistical equilibrium for those elements that are primarily responsible for line-blocking and blanketing, as well as an approximate treatment of the line-blocking itself, which is based on a simple statistical approach using suitable means of line opacities and emissivities. Both methods are validated by specific tests. Furthermore, we comment on our implementation of a consistent temperature structure.

In the second part, we concentrate on a detailed comparison with results from two codes used in alternative spectroscopical investigations, namely CMFGEN (Hillier & Miller 1998) and WM-Basic (Pauldrach et al. 2001). All three codes predict almost identical temperature structures and fluxes for $\lambda >$ 400 Å, whereas at lower wavelengths a number of discrepancies are found. Particularly in the He II continua, where fluxes and corresponding numbers of ionizing photons react extremely sensitively to subtle differences in the models, we consider any uncritical use of these quantities (e.g., in the context of nebula diagnostics) as unreliable. Optical H/He lines as synthesized by FASTWIND are compared with results from CMFGEN, obtaining a remarkable coincidence, except for the He I singlets in the temperature range between 36 000 to 41 000 K for dwarfs and between 31 000 to 35 000 K for supergiants, where CMFGEN predicts much weaker lines. Consequences of these discrepancies are discussed.

Finally, suggestions are presented as to adequately parameterize model-grids for hot stars with winds, with only one additional parameter compared to standard grids from plane-parallel, hydrostatic models.

Key words: methods: numerical - line: formation - stars: atmospheres - stars: early-type - stars: mass-loss

   
1 Introduction

The quantitative spectroscopy of massive stars with winds has made enormous progress due to the development of NLTE (non-local thermodynamic equilibrium) atmosphere codes that allow for the treatment of metal-line blocking and blanketing. With respect to both spectral range (from the extreme ultraviolet, EUV, to the infrared, IR) and metallicity of the analyzed objects (from SMC-abundances to Galactic center stars), a wide range in parameters can now be covered. Presently, five different codes are in use which have been developed for specific objectives, but due to constant improvements they can be applied in other contexts as well. In particular, these codes are CMFGEN (Hillier & Miller 1998), the "Potsdam-group'' code developed by W.R. Hamann and collaborators (for a status report, see Gräfener et al. 2002), the "multi-purpose model atmosphere code'' PHOENIX (Hauschildt & Baron 1999), WM-Basic (Pauldrach et al. 2001) and FASTWIND, which will be described here (see also Santolaya-Rey et al. 1997; Herrero et al. 2002, for previous versions).

The first three of these codes are the most "exact'' ones, since all lines (including those from iron-group elements) are treated in the comoving frame (CMF), which of course is a very time-consuming task. Moreover, since the first two of these codes have originally been designed for the analysis of the very dense winds from Wolf-Rayet stars, the treatment of the photospheric density stratification is approximative (constant photospheric scale-height). For several analyses this problem has been resolved by "coupling'' CMFGEN with the plane-parallel, hydrostatic code TLUSTY developed by Hubeny & Lanz (1995) (e.g., Bouret et al. 2003).

The multi-purpose code PHOENIX is mainly used for the analysis of supernovae and (very) cool dwarfs, but also a small number of hotter objects have been considered, e.g., the A-type supergiant Deneb (Aufdenberg et al. 2002). Due to this small number a detailed comparison with corresponding results is presently not possible, and, therefore, we will defer this important task until more material becomes available.

In contrast to all other codes that use a pre-described mass-loss rate and velocity field for the wind structure, the model atmospheres from WM-Basic are calculated by actually solving the hydrodynamical equations (with the radiative line-pressure being approximated within the force-multiplier concept, cf. Pauldrach et al. 1986; Castor et al. 1975) deep into the photosphere. Thus, this code provides a more realistic stratification of density and velocity, particularly in the transonic region (with the disadvantage that the slope of the velocity field cannot be manipulated if the wind does not behave as theoretically predicted). Since WM-Basic aims mainly at the prediction of EUV/UV fluxes and profiles, the bound-bound radiative rates are calculated using the Sobolev approximation (including continuum interactions), which yields "almost'' exact results except for those lines which are formed in the transonic region (e.g., Santolaya-Rey et al. 1997). Moreover, line-blocking is treated in an effective way (by means of opacity sampling throughout a first iteration cycle, and "exactly'' in the final iterations), so that the computational time is significantly reduced compared to the former three codes.

FASTWIND, finally, has been designed to cope with optical and IR spectroscopy of "normal'' stars with $T_{\rm eff}\ga 8500$ K, i.e., OBA-stars of all luminosity classes and wind strengths.

Since the parameter space investigated for the analysis of one object alone is large, comprising the simultaneous derivation of effective temperature  $T_{\rm eff}$, gravity $\log g$, wind-strength parameter $Q = \dot M/(R_\starv_\infty)^{1.5}$ (cf. Sect. 9), velocity field parameter $\beta$, individual abundances (most important: helium-abundance $Y_{\rm He}$) and also global background metallicity z, much computational effort is needed to calculate the large number of necessary models. This is one of the reasons why the samples which have been analyzed so far by both CMFGEN and WM-Basic are not particularly large, comprising typically five to seven objects per analysis (e.g., Martins et al. 2004; Hillier et al. 2003; Bouret et al. 2003, for recent CMFGEN-analyses; and Fullerton et al. 2000; Garcia & Bianchi 2004; Bianchi & Garcia 2002, for recent WM-Basic analyses).

Although the number of fit-parameters gets smaller when the wind-strength becomes negligible, a difference between the results from "wind-codes'' and plane-parallel, hydrostatic model atmospheres still remains: independent of the actual mass-loss rate, there will always be an enhanced probability of photon escape from lines in regions close to the sonic point and above, if a super-sonic velocity field is present. An example for the consequences of this enhanced escape is the He II ground-state depopulation in O-stars (Gabler et al. 1989), even though it is diminished by line-blocking effects compared to the original case studied with pure H/He atmospheres (see also Sect. 4.7).

With the advent of new telescopes and multi-object spectrographs, the number of objects that can be observed during one run has significantly increased (e.g., FLAMES attached to the VLT allows for observation of roughly 120 objects in parallel). An analysis of those samples will result in more reliable parameters due to more extensive statistics but remains prohibitive unless the available codes are considerably fast.

This is the motivation that has driven the development of FASTWIND. We have always considered speed to be of highest priority. The required computational efficiency is obtained by applying appropriate physical approximations to processes where high accuracy is not needed (regarding the objective of the analysis - optical/IR lines), in particular concerning the treatment of the metal-line background opacities.

Meanwhile, a number of analyses have been performed with our present version of FASTWIND, with significant sample sizes, of the order of 10 to 40 stars per sample (e.g., Massey et al. 2004; Repolust et al. 2004; Urbaneja et al. 2003; Massey et al. 2005; Trundle et al. 2004; Urbaneja 2004). Although the code has been carefully tested and first comparisons with results from CMFGEN and TLUSTY have been published (Herrero et al. 2002), a detailed description of the code and an extensive comparison have not been presented so far. Particularly the latter task is extremely important, because otherwise it is almost impossible to compare the results from analyses performed using different codes and to draw appropriate conclusions. An example of this difficulty is the discrepancy in stellar parameters if results from optical and UV analyses are compared. Typically, UV-spectroscopy seems to result in lower values for  $T_{\rm eff}$ than a corresponding optical analysis, e.g., Massey et al. (2005). Unless the different codes have been carefully compared, no one can be sure whether this is a problem related to either inadequate physics or certain inconsistencies within the codes.

This paper intends to answer part of these questions and is organized as follows: in Sect. 2 we give a brief overview of the basic philosophy of the code, and in Sect. 3 we describe the atomic data used as well as our treatment of metallicity regarding the flux-blocking background elements. Sections 4 and 5 give a detailed description of our approach to obtain the fast performance desired: Sect. 4 details the approximate NLTE solution for the background elements (which is applied if no consistent temperature structure is aimed at), and Sect. 5 describes our present method to tackle the problem of line-blocking. Both sections include important tests supporting the validity of our approach, particularly after a comparison with results from WM-Basic. Section 6 covers the problem of level inversions and how to deal with them, and Sect. 7 comprises the calculation of a consistent temperature structure. In Sect. 8, a detailed comparison with results from a grid of CMFGEN models is performed, and Sect. 9 suggests how to parameterize model-grids adequately and reports on first progress. In Sect. 10 we present our summary and an outlook regarding future work.

   
2 The code

The first version of the code (unblocked atmosphere/line formation) was introduced by Santolaya-Rey et al. (1997, hereafter Paper I), and has been significantly improved since. We distinguish between two groups of elements, the so-called explicit ones and the background elements.

The explicit elements (mainly H, He, but also C, N, O, Si, Mg in the B-star range, see below) are those used as diagnostic tools and are treated with high precision, i.e., by detailed atomic models and by means of CMF transport for the bound-bound transitions. In order to allow for a high degree of flexibility and to make use of any improvements in atomic physics calculations, the code is atomic data-driven with respect to these ions, as explained in Paper I: the atomic models, all necessary data and the information on how to use these data are contained in a user-supplied file (in the so-called DETAIL input form, cf. Butler & Giddings 1985) whereas the code itself is independent of any specific data.

The background ions, on the other hand, are those allowing for the effects of line-blocking/blanketing. The corresponding data originate from Pauldrach et al. (2001,1998) and are used as provided, i.e., in a certain, fixed form.

FASTWIND follows the concept of "unified model atmospheres'' (i.e., a smooth transition from a pseudo-hydrostatic photosphere to the wind) along with an appropriate treatment of line-broadening (Stark, pressure-) which is a prerequisite for the analysis of O-stars of different luminosity classes covering a variety of wind densities. Particularly and as already described in Paper I, the photospheric density consistently accounts for the temperature stratification and the actual radiation pressure, now by including both the explicit and the background elements.

The corresponding occupation numbers and opacities (of the background-elements) can be derived in two alternative ways:

a)

when the temperature stratification is calculated by means of NLTE Hopf parameters (see below), we apply an approximate NLTE solution for all background elements following the principal philosophy developed by Abbott & Lucy (1985), Schmutz (1991), Schaerer & Schmutz (1994) and Puls et al. (2000), where important features have now been improved (cf. Sect. 4). Particularly, the equations of approximate ionization equilibrium have been re-formulated to account for the actual radiation field as a function of depth and frequency, and a consistent iteration scheme regarding the coupling of the rate equations and the radiation field has been established to avoid the well-known convergence problems of a pure Lambda Iteration (Sect. 4.6).

b)

when the T-stratification is calculated from first principles, the complete set of rate equations is solved almost "exactly'' for the most abundant background elements (C, N, O, Ne, Mg, Si, S, Ar, Fe, Ni, if not included as explicit ions), employing the Sobolev approximation for the net radiative rates (with actual illuminating radiation field). The remaining background elements, on the other hand, remain to be treated by the approximation as outlined in a).

In order to account for the effects of line-blocking, we use suitable means for the line opacities, averaged over a frequency interval of the order of 1000... 1500 km s^-1, and appropriate emissivities (Sect. 5).

Finally, the temperature stratification can be calculated in two different ways. If one is exclusively interested in an optical analysis, the concept of NLTE-Hopf parameters (cf. Paper I) is still sufficient, if the background elements are accounted for in a consistent way, i.e., have been included in the particular models from which these parameters are derived. Since this method is flux-conservative, the correct amount of line-blanketing is "automatically'' obtained. Note that for optical depths $\tau_{\rm R}$ $\la$ 0.01 a lower cut-off temperature is defined, typically at $T_{\rm min} = 0.6 T_{\rm eff}$.

Alternatively, the new version of FASTWIND allows for the calculation of a consistent temperature, utilizing a flux-correction method in the lower atmosphere and the thermal balance of electrons in the outer one (Sect. 7). As has been discussed, e.g., by Kubát et al. (1999), the latter method is advantageous compared to exploiting the condition of radiative equilibrium in those regions where the radiation field becomes almost independent of  $${T_{\rm e}}$$. Particularly for the IR-spectroscopy, such a consistent T-stratification is important, since the IR is formed above the stellar photosphere in most cases and depends critically on the run of  $${T_{\rm e}}$$ in those regions, where our first method is no longer applicable.

   
3 Atomic data and metallicity

Explicit elements.

In order to obtain reliable results also in the IR, we have significantly updated our H- and He-models compared to those described in Paper I. Our present H and He II models consist of 20 levels each (vs. 10 and 14 in the previous version, respectively), and He I includes levels until n=10, where levels with $n=8{\ldots} 10$ have been packed (previous version: 8 levels, packed from 5...8). Further information concerning cross-sections etc. can be found in Jokuthy (2002). Present atomic models for metals have been accumulated from different sources, mainly with respect to an analysis of B-stars, i.e., for ionization stages  II and  III, except for Mg ( I, II) and Si ( II, III, IV). Information on our Si atomic model can be found in Trundle et al. (2004), and on the other metals incorporated so far (C, N, O, Mg) in Urbaneja (2004).

Background elements.

The atomic data for background elements come from Pauldrach et al. (2001,1998), who have given a detailed description of the various approaches and sources. These data comprise the elements from hydrogen to zinc (except Li, Be, B and Sc which are too rare to affect the background opacity) with ionization stages up to VIII. The number of connecting lines (lower and upper level present in the rate equations) is of the order of 30 000, and the number of lines where only the lower level is present is $4.2 \times 10^6$. The former group of lines is used to solve the rate equations, whereas the latter is used to derive the metal-line background opacities (cf.  Sect. 5). In addition to bound-free cross-sections and gf-values, there is also detailed information about the collision-strengths for the most important collisional bound-bound transitions in each ion.

Metallicity.

The abundances of the background elements are taken from the solar values provided by Grevesse & Sauval (1998, and references therein). For different "global'' metallicities, $z=Z/Z_\odot$, these abundances are scaled proportionally with respect to mass ratios, e.g., by 0.2 for the SMC and by 0.5 for the LMC (although these values are certainly disputable, e.g. Massey et al. 2004, and references therein).

A particular problem (independent of the actual value of z) appears in those cases when the He/H ratio becomes non-solar. In this case, we retain the specific relative mass fractions of the other elements, which of course has a significant effect on the number ratios. Although this procedure is not quite right, it preserves at least the overall mass fraction of the metals, particularly the unprocessed iron group elements, which are most important for the line-blocking. Further comments on the validity of this procedure have been given by Massey et al. (2004). We briefly mention a comparison to evolutionary calculations from Schaerer et al. (1993) performed by P. Massey (priv. comm.):

For the 120 $M_\odot$ track at Z=0.008 (roughly the LMC metallicity), Zstays essentially unchanged in the core until the end of core H burning, even though the mass fractions of C and N increase while O decreases: at a number ratio $Y_{\rm He}$ = 2 (i.e., the mass ratio Y has changed from 0.265 to 0.892), the value for Z has changed insignificantly from 0.0080 to 0.0077, and even more interestingly, the mass fraction of the sum of C, N, O, and Ne has essentially changed in the same way (0.0075 to 0.0070), even though the actual mass fraction of N has more than doubled.

   
4 Background elements: Approximate NLTE occupation numbers

To save significant computational effort, the occupation numbers of the background elements are calculated by means of an approximate solution of the NLTE rate equations. Such an approach has been successfully applied in a variety of stellar atmosphere calculations, e.g., to derive the radiative acceleration of hot star winds (Lucy & Abbott 1993; Abbott & Lucy 1985) and for the spectroscopy of hot stars (Schaerer & Schmutz 1994; Schmutz 1991) and Supernova remnants (Mazzali 2000; Lucy 1999; Mazzali & Lucy 1993). Puls et al. (2000) have used this method for an examination of the line-statistics in hot star winds, by closely following a procedure discussed by Springmann (1997) which in turn goes back to unpublished notes by L. Lucy.

One might argue that such an approximate approach can poorly handle all the complications arising from sophisticated NLTE effects. However, in the following we will show that the approximate treatment is able to match "exact'' NLTE calculations to an astonishingly high degree, at least if some modifications are applied to the original approach. Moreover, the calculated occupation numbers will not be used to synthesize line-spectra, but serve "only'' as lower levels for the line-opacities involved in the blocking calculations.

Actually, the major weakness of the original approach is the assumption of a radiation field with frequency independent radiation temperatures  $T_{\rm rad}$. Since particularly the difference in radiation temperatures at strong ionization edges is responsible for a number of important effects, we have improved upon this simplifiction by using consistent radiation temperatures (taken from the solution of the equations of radiative transfer). As we will see in the following, this principally minor modification requires a number of additional considerations.

4.1 Selection of levels

One of the major ingredients entering the approximate solution of the rate equations is a careful selection of participating atomic levels. In agreement with the argumentation by Abbott & Lucy (1985) only the following levels are used:

• the ground-state level;
• all meta-stable levels (from equal and different spin systems), denoted by "M'';
• all excited levels which are coupled to the ground-state via one single permitted transition where this transition is the strongest among all possible downward transitions; in the following denoted as subordinate levels "N'';
• all excited levels coupled to one of the meta-stable levels $m \in M$ in a similar way (subordinate levels "S'').

In the above definition, the term "strongest'' refers to the Einstein coefficients A_ji. All other levels are neglected, since their population is usually too low to be of importance and cannot be approximated by simple methods.

4.2 Ionization equilibrium

In order to allow for a fast and clearly structured algorithm, we allow only for ionizations to and recombinations from the ground-state of the next higher ion, even if this is not the case in reality. Due to this restriction and by summing over all line-processes an "exact'' rate equation connecting two neighboring ions is derived which exclusively consists of ionization/recombination processes. In the following, we will further neglect any collisional ionization/recombination processes, which is legitimate in the context considered here, namely in the NLTE-controlled atmospheric regime of hot stars. (In the lowermost, LTE dominated part of the atmosphere, $\tau_{\rm R} > 2/3$, we approximate the occupation numbers a priori by LTE conditions).

	Figure 1: Ratio of ionization to recombination rate coefficients: relative error between "exact'' ratios (Eq. (12)) and approximate ones (Eq. (15), with $\beta = 1$ and s = 2) as a function of  $T_{\rm rad}/T_{\rm e}$, for different combinations of ($\beta , s$ ). The error decreases for even higher ionization energies.
Open with DEXTER

At first, let us consider an ion with only one spin system, e.g., a hydrogenic one. In this case, the ionization equilibrium becomes

$$% \sum_i n_{i} R_{i \kappa} = n_{\kappa} \sum_i \left(\frac... ...n_\kappa} \right)^* R_{\kappa i} = \sum_i n_i^* R_{\kappa i},$$ (1)

with n_i the occupation numbers of the lower ionization stage, $n_{\kappa}$ the (ground-state) occupation number of the higher ion, the asterisks denoting LTE-conditions (at the actual electron density, $(n_i/n_{\kappa})^* = n_e \Phi(T)$, cf. Mihalas 1975, Sect. 5) and ionization/recombination rate coefficients

$$% R_{i \kappa} = \int_{\nu_i}^\infty \frac{4 \pi J_\nu}{h \nu} a(\nu)~ {\rm d}\nu$$ (2)

$$% R_{\kappa i} = \int_{\nu_i}^\infty \frac{4 \pi}{h \nu} \l... ... J_\nu \right) {\rm e}^{-h \nu/kT_{\rm e}} a(\nu)~{\rm d}\nu.$$ (3)

$J_\nu$ is the mean intensity, $a(\nu)$ the ionization cross-section and all other symbols have their usual meaning. Once more, within our above approximation (ionization to ground state only), Eq. (1) is "exact'' and does not depend on any assumption concerning the bound-bound processes (radiative or collisional; optically thick or thin) since the corresponding rates drop out after summation.

By introducing the recombination coefficient $\alpha_i$ defined in the conventional way,

$$% n_\kappa n_e \alpha_i = n_i^* R_{\kappa i},$$ (4)

the ionization equilibrium can be reformulated

$$% \sum_i n_{i} R_{i \kappa} = n_{\kappa} n_{\rm e} \sum_i \alpha_i,$$ (5)

and we extract all quantities referring to the ground-state of the lower ion,

$$% n_{\kappa} n_{\rm e} = \frac{1}{\sum_i \alpha_i} n_1 R_{1 \... ...rac{\sum_{i>1} n_i R_{i \kappa}}{n_1 R_{1 \kappa}}\right)\cdot$$ (6)

Finally, inserting the ground-state recombination coefficient $\alpha_1$ (cf. Eq. (4)) on the rhs, we obtain the ionization equilibrium expressed as the ratio of two neighboring ground-states,

$$% \frac{n_\kappa}{n_1} = \left(\frac{n_\kappa}{n_1} \right)^*... ...n_{i}}{n_{1}} \frac{R_{i \kappa}}{R_{1 \kappa}} \right\}\cdot$$ (7)

Note that this ratio depends on the radiation field, the actual electron density and temperature and on the excitation within the lower ion, which will be discussed in the next subsection. Note also that all that follows is "only'' a simplification of this equation.

So far, our derivation and the above result are identical to previous versions of the approximate approach. From now on, however, we will include the frequency dependence of the radiation field. To this end, we describe the ionization cross-sections by the Seaton approximation (Seaton 1958), which is not too bad for most ions,

$$% a(\nu) = a_i \left( \beta \left( \frac{\nu_i}{\nu}\right)^s + (1 - \beta)\left( \frac{\nu_i}{\nu}\right)^{s+1}\right)\cdot$$ (8)

Writing the mean intensity as $J_{\nu}(r) = W(r) B_{\nu}(T_{\rm rad}(\nu, r))$ with dilution factor W and neglecting the stimulated emission in the recombination integral (valid for all important ionization edges), we obtain (radial dependence of all quantities suppressed in the following)

 

$$R_{i \kappa} = \frac{8 \pi W}{c^2} \left(\frac{k T_{{\rm r},i}}{h}\right)^3 a_i~ \mathcal{F}(x_{-{\rm r},i};\beta,s)$$ (9)

$$R_{\kappa i} = \frac{8 \pi}{c^2} \left(\frac{k T_{\rm e}}{h}\right)^3 a_i~ \mathcal{F}(x_{{\rm e},i};\beta,s)$$ (10)

with $x_{\rm r,i}=h\nu_i/kT_{{\rm r},i}$, $x_{{\rm e},i}=h\nu_i/kT_{\rm e}$ and

$$% \mathcal{F}(x;\beta,s) = \beta~ x^s~ \Gamma(3-s,x) + (1-\beta)~x^{1+s}~\Gamma(2-s,x).$$ (11)

We have assumed $T_{\rm rad}(\nu) =: T_{{\rm r},i}$ to be constant over the decisive range of the ionizing continuum $\nu \ga \nu_i$, since only those frequencies close to the edge are relevant. In other words, each transition is described by a unique radiation temperature. In the above equation, the incomplete Gamma-function  $\Gamma(a,x)$ has been generalized to also include negative parameters, $a \le 0$. The ratio of ground-state ionization/recombination rate coefficients is thus given by

$$% \frac{R_{1 \kappa}}{R_{\kappa 1}} = W~\left(\frac{\displa... ...{F}(x_{\rm r,1};\beta,s)} {\mathcal{F}(x_{\rm e,1};\beta,s)},$$ (12)

i.e., is independent of the actual value of the cross-section at the threshold, a_i. Although this expression is rather simple, it requires the somewhat time-consuming evaluation of the incomplete Gamma-functions. To keep things as fast as possible, we generally use the parameter set ( $\beta=1, s=2$) instead of the actual parameters which results in a particularly simple function  $\mathcal{F}$,

$$% \mathcal{F}(x;1,2)=x^2 \exp(-x).$$ (13)

Note that these parameters do not correspond to the hydrogenic cross-section, which would be described by s=3. Using this parameter set, the ionization/recombination rates simplify to

 

$$R_{i \kappa} = \frac{8 \pi W}{c^2} \left(\frac{k T_{{\rm r},i}}{h}\right) a_i \nu_i^2 {\rm e}^{-h \nu_i/k T_{\rm r,i}}$$

$$R_{\kappa i} = \frac{8 \pi}{c^2}\left(\frac{kT_{\rm e}}{h}\right) a_i \nu_i^2 {\rm e}^{-h \nu_i/kT_{\rm e}},$$ (14)

and the ratio $R_{1 \kappa}/R_{\kappa 1}$ becomes

$$% \frac{R_{1 \kappa}}{R_{\kappa 1}} = W ~ \frac{\displaystyl... ...style {T_{\rm r,1}}} - \frac{1}{T_{\rm e}} \right)\right]\cdot$$ (15)

We have convinced ourselves that this approximation leads to acceptable errors of the order of 10%, cf. Fig. 1. Furthermore, we define the following quantities, where $\zeta$ is the ratio of ground-state to total recombination coefficient,

$$% \frac{\alpha_1}{\sum_i \alpha_i} = \zeta, \sum_{i \in {\rm ... ...1}} \frac{R_{i \kappa}}{R_{1 \kappa}} \equiv C_{{\rm M(N,S)}}.$$ (16)

Any ratio $\alpha_i/\alpha_j$ (particularly, the case j=1 and thus $\zeta$) is independent of the temperature and depends exclusively on atomic quantities, namely cross-section, transition frequency and statistical weight, a fact which follows from Eqs. (4) and (14):

 

$$% \frac{\alpha_i}{\alpha_j} \left(\frac{n_i}{n_j} \right)^* \frac{R_{\kappa i}}{R_{\kappa j}}... ...nu_i}{\nu_j} \right)^2 \exp \left[ -\frac{h(\nu_i - \nu_j)}{kT_{\rm e}} \right]$$ =

$$\frac{a_i}{a_j} \left(\frac{\nu_i}{\nu_j} \right)^2 \frac{g_i}{g_j}\cdot$$ = (17)

Collecting terms, our approximate ionization equilibrium finally reads

 

$$% \frac{n_\kappa}{n_1} \left(\frac{n_\kappa}{n_1} \right)^* W~\frac{\displaystyle {T_{\r... ...left(\frac{1}{\displaystyle {T_{\rm r,1}}} - \frac{1}{T_{\rm e}} \right)\right]$$ =

$$\times\zeta~(1+C_{\rm N}+C_{\rm M}+C_{\rm S})$$

$$\left(\frac{n_\kappa}{n_1} \right)^*_{\displaystyle {T_{\rm r,1}}... ...rm e}}{\displaystyle {T_{\rm r,1}}}}~ \zeta ~(1+C_{\rm M}+C_{\rm N}+C_{\rm S}),$$ = (18)

where the second variant uses the LTE ratio evaluated at the actual electron-density and radiation temperature of the ionizing continuum.

   
4.3 Excitation

The remaining step concerns the term in the bracket above, i.e., the approximate calculation of the excitation inside the lower ion (which, of course, is also required to calculate the partition functions). For consistency, frequencies (energies) are still defined with respect to the ionization threshold, i.e., line frequencies have to be calculated from $\nu_{ij} = \nu_i - \nu_j >0$ instead of the usual definition (upper - lower) which would refer to excitation energies.

   
4.3.1 Meta-stable levels

We begin with the occupation numbers of meta-stable levels which can be populated via excited levels or via the continuum (see also Abbott & Lucy 1985).

Population via excited levels.

Denoting the excited level by j, considering the fact that this excited level is fed by the ground state (otherwise it would not exist in our level hierarchy) and neglecting collisional processes, the population can be approximated by

$$% \frac{n_m}{n_1} = \frac{W \left(\frac{n_j}{n_1} \right)^*_{... ...\rm r,}mj}}}}, \qquad m \in {\rm M} \quad (j \in {\rm N} > m)$$ (19)

(see also Eq. (26) with $\delta \approx 0$), where the dilution factor cancels out. In the following, we have to distinguish between two cases: the meta-stable level lies either close to the ground-state or close to the excited level, cf. Fig. 2.

	Figure 2: Population of meta-stable levels via excited ones (see text).
Open with DEXTER

Case a: low lying meta-stable level. The transition frequencies of both transitions are fairly equal, $\nu_{1j} \approx \nu_{mj}$, i.e., $${T_{{\rm r,1}j}}\approx \displaystyle {T_{{\rm r,}mj}}$$, and we find

$$% \frac{n_m}{n_1} \approx \frac{g_m}{g_1} \exp\left(-\frac{... ...frac{n_m}{n_1} \right)^*_{\displaystyle {T_{{\rm r,1}j}}}\cdot$$ (20)

Note that the population is controlled by the radiation field  $${T_{{\rm r,1}j}}$$, i.e., from frequencies much larger than the "excitation energy'', $\nu_{1m}$.

Case b: high lying meta-stable level. Now we have $\nu_{mj} < \nu_{1j} \approx \nu_{1m}$,

$$% \frac{n_m}{n_1} \!\approx\! \frac{g_m}{g_1} \exp \left( -\... ...ft(\frac{n_m}{n_1} \right)^*_{\displaystyle {T_{{\rm r,1}m}}},$$ (21)

and the population depends on the radiation field at (or close to) the excitation energy.

Population via continuum.

The third case comprises a population via the continuum which will only be treated in a crude approximation, where a correct evaluation will be given later (Sect. 4.5). If we neglect for the moment the influence of any meta-stable and excited levels, we find from Eq. (18) with $\zeta \rightarrow1, C_{{\rm N(M,S)}} \rightarrow0$

 

$$% \frac{n_m}{n_1} = \frac{\left(\frac{n_{\kappa}}{n_1} \right)} {\left(\frac{n_{\kappa}}{n_m} \right)} \approx \frac{\left(\frac{n_{\kappa}}{n_1} \right)^*_{\displaystyle {T_{\... ...,m}}}} \sqrt{\frac{\displaystyle {T_{{\rm r},m}}}{\displaystyle {T_{\rm r,1}}}}$$

$$\frac{g_m}{g_1}~ \frac{\displaystyle {T_{\rm r,1}}}{\displaystyle... ...T_{\rm r,1}}} - \frac{\nu_m}{\displaystyle {T_{{\rm r},m}}} \right)\right]\cdot$$ = (22)

Note that all three cases converge to the identical result

$$% \frac{n_m}{n_1} \rightarrow\left(\frac{n_m}{n_1}\right)^*_{T_{\rm rad}} \quad {\rm for} \quad T_{\rm rad}= {\rm const}.$$ (23)

which is quoted by Abbott & Lucy (1985).

In order to continue our calculation of $C_{\rm M}$, we find from Eqs. (14) and (17)

$$% \frac{R_{m \kappa}}{R_{1 \kappa}} = \frac{\displaystyle {T... ... \frac{\nu_1}{\displaystyle {T_{\rm r,1}}} \right)\right]\cdot$$ (24)

Multiplying by n_m/n₁ we find that for the three cases Eqs. (20)-(22)

$$% \left(\frac{n_m R_{m \kappa}}{n_1 R_{1 \kappa}} \right)_{m... ... \displaystyle\frac{\alpha_m}{\alpha_1}, \end{array} \right.$$ (25)

respectively. As mentioned before, the result for the third case (population over continuum) is only a crude approximation, which is also evident from the fact that it depends only on atomic quantities but not on any radiation temperature.

4.3.2 Subordinate levels

Due to our definition of subordinate levels their population can be approximated by a two-level-atom Ansatz (between ground-state j=1 and subordinate level  $i \in {\rm N}$ or between meta-stable level  $j \in {\rm M}$ and subordinate level  $i \in {\rm S}$), such that the population can be expressed by

$$% \frac{n_i}{n_j} \!=\! W(1-\delta)\left( \frac{n_i}{n_j} \r... ...e {T_{\rm e}}}, \quad i \in {\rm N (S)},\; j \in {\rm 1 (M)}$$ (26)

where $\delta$ is the parameter expressing the competition between thermalization ( $\delta \rightarrow1$) and local escape (in the Sobolev approximation),

$$% \delta = \frac{\epsilon}{\epsilon(1-\beta) + \beta}\cdot$$ (27)

$\epsilon$ is the usual LTE parameter in a two-level atom,

$$% \epsilon = \frac{C_{ji}}{A_{ji} + C_{ji}},$$ (28)

with collisional de-excitation rate C_ji and Einstein-coefficient A_ji. $\beta$ is the local escape probability in the Sobolev approximation,

$$% \beta = \frac{1}{2} \int_{-1}^{1} \frac{1 - {\rm e}^{-\tau_S(\mu)}}{\tau_S(\mu)}~ {\rm d}\mu,$$ (29)

and the illuminating radiation field is approximated by

$$% \beta_{\rm c} I_{\rm c} = \frac{1}{2} \int_{\mu_{\ast}}^{1}... ...}\mu \approx W B_{\nu}(\displaystyle {T_{{\rm r},ij}})~\beta.$$ (30)

Note that our approximation (26) neglects any coupling to the continuum inside the resonance zone. By means of Eq. (17), the individual terms comprising $C_{\rm N}$ can be calculated from

$$% \left(\frac{n_{i}}{n_{1}} \frac{R_{i \kappa}}{R_{1 \kappa}} \right)_{i ~\in~ {\rm N}} \frac{\alpha_i}{\alpha_1} \frac{T_{{\rm r},i}}{T_{\rm r,1}}~ \exp... ...k} \left(\frac{\nu_i}{T_{{\rm r},i}} - \frac{\nu_1}{T_{\rm r,1}} \right)\right]$$ =

$$\times \left[W (1\!-\!\delta_{1i}) \exp \left(-\frac{h \nu_{1i}}{... ...ght) \!+ \!\delta_{1i} \exp \left(-\frac{h \nu_{1i}}{kT_{\rm e}}\right) \right]$$ (31)

whereas the components of $C_{\rm S}$ are described by

$$% \left(\frac{n_{i}}{n_{1}} \frac{R_{i \kappa}}{R_{1 \kappa}} \right)_{i ~\in ~{\rm S}} \left(\frac{n_{m}}{n_{1}}\right)_{m \in {\rm M}} \cdot \left(\frac{n_{i}}{n_{m}} \frac{R_{i \kappa}}{R_{1 \kappa}} \right)_{i \in {\rm S}}$$ =

$$\left(\frac{n_m}{n_1}/\frac{g_m}{g_1}\right) \times \frac{\alpha_... ...k} \left(\frac{\nu_i}{T_{{\rm r},i}} - \frac{\nu_1}{T_{\rm r,1}} \right)\right]$$ =

$$\times \left[W (1\!-\!\delta_{mi}) \exp \left(-\frac{h \nu_{mi}}{... ...ght) \!+\! \delta_{mi} \exp \left(-\frac{h \nu_{mi}}{kT_{\rm e}}\right) \right]$$ (32)

(with (n_m/n₁) taken from Eqs. (20)-(22), respectively). Obviously, the population of subordinate levels is controlled by at least three different radiation temperatures (ionization from the considered level, ionization from the connected lower level and excitation due to line processes).

4.4 Limiting cases

In the following, we will consider some limiting cases which have to be reproduced by our approach.

Constant radiation temperature, no collisional excitation

are the assumptions underlying the description by Springmann (1997) and Puls et al. (2000) on the basis of Lucy's unpublished notes. With $T_{\rm rad}= {\rm const}$, meta-stable levels are populated via

$$% \left(\frac{n_{m}}{n_{1}} \frac{R_{m \kappa}}{R_{1 \kappa}}... ...{\alpha_m}{\alpha_1} \qquad (\nu_m - \nu_1 + \nu_{1m} = 0!),$$ (33)

independent of the actual feeding mechanism. With $\delta = 0$ (only radiative line processes), we thus obtain for the population of subordinate levels (both $i \in {\rm N}$ and $i \in {\rm S}$!)

$$% \left (\frac{n_{i}}{n_{1}} \frac{R_{i \kappa}}{R_{1 \kappa}} \right)_{i ~\in~ {\rm N, S}} = W \frac{\alpha_i}{\alpha_1}\cdot$$ (34)

Thus, for constant radiation temperatures, it does not play a role in how the meta-stable levels are populated, and whether subordinate levels are connected to the ground-state or to a meta-stable level. Only the corresponding recombination coefficient is of importance and the fact that subordinate levels suffer from dilution (since they are fed by a diluted radiation field), whereas for meta-stable levels this quantity cancels out (cf. Abbott & Lucy 1985). In total, our simplified ionization equilibrium then becomes

$$% \frac{n_\kappa}{n_1} = \left(\frac{n_\kappa}{n_1} \right)... ...sum_{i ~\in~ {\rm N,S}} \frac{\alpha_i}{\alpha_1} \right)\cdot$$ (35)

If we define $\eta$ as the fraction of recombination coefficients for all meta-stable levels,

$$% \eta=\frac{\sum_{i \in {\rm M}} \alpha_i}{\sum_{i} \alpha_i}$$ (36)

we find

$$% \sum_{i ~\in ~{\rm M}} \frac{\alpha_i}{\alpha_1} = \frac{1... ...ha_1}}{\sum_i \frac{\alpha_i}{\alpha_1}} = \frac{\eta}{\zeta}$$ (37)

and

$$% \sum_{i~ \in~ {\rm N,S}} \frac{\alpha_i}{\alpha_1} = \frac{1-\eta-\zeta}{\zeta},$$ (38)

and the ionization equilibrium can be described in a very compact way,

$$% \frac{n_\kappa}{n_1} = \left(\frac{n_\kappa}{n_1}\right)^*_... ...rm rad}}}\; \bigl(~\zeta + \eta + W (1-\eta - \zeta) ~\bigr),$$ (39)

which indeed is the result of the previous investigations mentioned above. If we further prohibit all ionizations from meta-stable and subordinate levels, i.e. allow for

Ionization/recombination only from and to the ground-state,

we find with $\zeta =1$ and $\eta =0$

$$% \frac{n_\kappa}{n_1} \left(\frac{n_\kappa}{n_1}\right)^*_{T_{\rm rad}} W \sqrt{\frac{T_{\rm e}}{T_{\rm rad}}}$$ =

$$\frac{2 g_\kappa}{g_1} \frac{1}{n_{\rm e}} \left(\frac{2 \pi m_{\... ...ft(-\frac{h \nu_1}{kT_{\rm rad}}\right) W \sqrt{\frac{T_{\rm e}}{T_{\rm rad}}},$$ = (40)

which is a well-known result and also valid for the case where all lines are optically thick and in detailed balance, e.g., Abbott (1982). The

LTE-case is recovered independently from the specific values of $\zeta$ and $\eta$ in the lowermost atmosphere, when the dilution factor approaches unity, W=1, and the radiation field becomes Planck, $T_{\rm rad}\rightarrow\displaystyle {T_{\rm e}}$. In this case, the ionization balance becomes

$$% \frac{n_\kappa}{n_1} = \left(\frac{n_\kappa}{n_1}\right)^*_... ...ta - \zeta)) = \left(\frac{n_\kappa}{n_1}\right)^*_{T_{\rm e}}$$ (41)

and for the excitation we have

$$% \frac{n_i}{n_1} = \left( \frac{n_i}{n_1} \right)^*_{\displaystyle {T_{\rm e}}}, \quad i \in {\rm M, N, S}.$$ (42)

   
4.5 Different spin systems

The last problem to overcome is the presence of different spin systems, a problem already encountered for He I. Our approximation is to consider the different systems as completely decoupled (except if strong inter-combination lines are present, see below), since a coupling via collisional inter-combination is effective only at high densities (i.e., in or close to LTE, which is treated explicitly in our procedure anyway).

Then for each of the separate multiplets, the ionization equation can be calculated independently. The different subsystems are defined in the following way:

$\bullet$

the first subsystem includes all levels coupled to the ground-state plus those meta-stable levels fed from higher-lying (subordinate) levels (case a/b in Sect. 4.3.1). In this way, we include also systems of different spin which are connected to the ground-state system via strong inter-combination lines, a condition which is rarely met. In total, the ground-state subsystem includes the levels $i \in$ 1, N, M', S', where M' comprises all case a/b meta-stable levels and S' those excited levels which are coupled to M'. For reasons of convenience, we will denote this set of levels by (1, N');

$\bullet$

a second group of j subsystems comprises

• systems of different spin decoupled from the ground-state;
• "normal'' meta-stable levels populated via continuum processes (poorly approximated so far) and excited levels coupled to those.

Both groups can be treated in a similar way and are also identified in a similar manner, namely from the condition that the lowest state of these systems is meta-stable and not fed from higher-lying levels. Each subsystem comprises the "effective'' ground state $m_j \in$ M'' (either different spin or fed by continuum) and coupled levels, $i \in$ S_j''.

Once more, j is the number of meta-stable levels per ion that are not fed by higher-lying levels. Note that for a single spin-system with meta-stable levels, there are now 1 +j different subsystems if jcontinuum fed meta-stable levels were present. Note also that by using this approach we neglect a possible coupling of two or more non-ground-state multiplets via strong inter-combination lines, if there were any.

Because of the assumed decoupling, for each subsystem we can write down the appropriate ionization equation. For the ground-state system, we have

  

$$\frac{n_\kappa}{n_1} \left(\frac{n_\kappa}{n_1} \right)^*_{\displaystyle {T_{\rm r,1}}... ...\sqrt{\frac{T_{\rm e}}{\displaystyle {T_{\rm r,1}}}}~ \zeta_1 ~(1+C_{{\rm N}'})$$ = (43)

$$\zeta_1 \frac{\alpha_1}{\sum_{(1,{\rm N}')} \alpha_i},\qquadC_{N'} = \sum_{i \in {\rm N'}} \frac{n_{i}}{n_{1}} \frac{R_{i \kappa}}{R_{1 \kappa}}$$ = (44)

where, again, N' comprises the "old'' levels $\in$ N, M' and S'. Note the difference between $\zeta_1$ and $\zeta$ from Eq. (16).

For each of the j additional subsystems, we obtain by analogy

 

$$\frac{n_\kappa}{n_{m_j}} = \left(\frac{n_\kappa}{n_{m_j}} \right)... ...j}}}}~\frac{\alpha_{m_j}}{\alpha_1}~ \zeta_{m_j}~\left(1+C_{{\rm S}''_j}\right)$$ (45)

$$\zeta_{m_j} = \frac{\alpha_1}{\sum_{({m_j},{\rm S}''_j)} \alpha_i... ...{i \in {\rm S''_j}} \frac{n_{i}}{n_{m_j}} \frac{R_{i \kappa}}{R_{{m_j} \kappa}}$$ (46)

and S''_j comprises all levels coupled to m_j. The individual components of  $C_{{\rm N}'}$ and  $C_{{\rm S}''_j}$ are calculated as described in Sect. 4.3. Dividing Eq. (43) by Eq. (45), we find for the ratios  (n_mj/n₁) (required, e.g., for calculating the partition functions),

$$% \frac{n_{m_j}}{n_1} = \frac{(n_\kappa/n_1)^*_{\displaystyl... ...m_j}} ~ \left(\frac{1+C_{{\rm N}'}}{1+C_{{\rm S}''_j}}\right),$$ (47)

or, explicitly written,

$$% \frac{n_{m_j}}{n_1} = \frac{g_{m_j} \displaystyle {T_{\rm r... ..._{m_j}} \left(\frac{1+C_{{\rm N}'}}{1+C_{{\rm S}''_j}}\right).$$ (48)

The last equation is the "correct approximation'' for continuum-fed meta-stable levels. On the one hand, if the ion consists of the ground-state plus a number of meta-stable levels alone, we would have $C_{{\rm N}'} = C_{{\rm S}''_j}= 0$, $\zeta_1 = 1$ and $\zeta_{m_j} = \alpha_1 / \alpha_{m_j}$. In this case, Eqs. (48) and (22) would give identical results, which shows that both approaches are consistent under the discussed conditions. But as already pointed out, Eq. (22) is highly approximative if a variety of levels are involved, and the occupation numbers always should be calculated according to Eq. (48).

The major difference to our former approach (one spin system only) is the following. In approach "one'', the ground-state population, $n_{\kappa}/n_1$, is affected by all meta-stable levels, whereas in approach "two'' only those meta-stable levels coupled to the ground-state system via higher levels have an influence.

Constant radiation temperature, no collisional excitation.

Concerning the limiting case where $T_{\rm rad}= {\rm const}$ and $\delta = 0$, Eq. (39) remains valid if we account for the different "normalization'', i.e., if we replace $\zeta$ by $\zeta_1$ and include into $\eta$ only those meta-stable levels that are populated via excited levels:

$$% \frac{n_\kappa}{n_1} = \left(\frac{n_\kappa}{n_1}\right)^*_... ...}}\; \bigl(~\zeta_1 + \eta_1 + W (1-\eta_1 - \zeta_1) ~\bigr)$$ (49)

with

$$% \zeta_1= \frac{\alpha_1}{\sum_{(1,{\rm N}')} \alpha_i}, \qq... ..._{i \in {\rm M'}} \alpha_i}{\sum_{(1,{\rm N}')} \alpha_i}\cdot$$ (50)

Inside the individual sub-systems we then obtain

$$% \frac{n_\kappa}{n_{m_j}} = \left(\frac{n_\kappa}{n_{m_j}}\r... ...a') ~\bigr), \; \zeta'=\frac{\alpha_{m_j}}{\alpha_1}\zeta_{mj}$$ (51)

which immediately indicates the correct thermalization for W=1 and  $T_{\rm rad}\rightarrow\displaystyle {T_{\rm e}}$. After dividing Eqs. (49) by (51), we find for the population of  (n_mj/n₁) in the same limit

$$% \frac{n_{m_j}}{n_1} = \frac{g_{m_j}}{g_1} \exp \left( -\f... ... + W (1-\eta_1 - \zeta_1)} {\zeta' + W (1-\zeta')}\right)\cdot$$ (52)

First, we obtain the correct population in LTE when  $W \rightarrow1$. Second, the difference to our crude approximation in Sect. 4.3.1 becomes obvious: the quasi-LTE ratio (23) has to be multiplied by the last factor in the above equation to obtain consistent populations. This factor (which can be lower or higher than unity) becomes unity only when $W \rightarrow1$ (i.e., in the lower atmosphere) or for $\zeta_1 = \zeta' = 1$, i.e., in those cases where only ground-state and meta-stable levels are present, as already discussed above.

   
4.6 Accelerated Lambda iteration

To overcome the well-known problems of the Lambda-iteration when coupling the rate-equations with the equation of radiative transfer, we apply the concept of the Accelerated Lambda Iteration (ALI, for a review see Hubeny 1992) to obtain a fast and reliable convergence of the solution. Since our rate-equations have been formulated in a non-conventional way and since the radiation field is expressed in terms of local, frequency-dependent radiation temperatures, the procedure has to be modified somewhat, and we will describe the required re-formulations as follows (for a comparable implementation see also de Koter et al. 1993).

At first, assume that only one bound-free opacity is present, i.e., the radiation field is controlled by the opacity of the considered transition i (no overlapping continua present). In this case, the usual ALI formulation for the mean intensity $J_{\nu}^n$ at iteration cycle n is given by

$$% J_{\nu}^n \rightarrow J_{\nu}^{n-1} + \Lambda_{\nu}^{\ast}\left(S_i^n - S_i^{n-1}\right)$$

$$\Delta J_{\nu}+ \Lambda_{\nu}^{\ast}S_i^n \quad {\rm with} \quad \Delta J_{\nu}= J_{\nu}^{n-1} - \Lambda_{\nu}^{\ast}S_i^{n-1},$$ = (53)

where S_i is the continuum source-function for transition i and  $\Lambda_{\nu}^{\ast}$the corresponding Approximate Lambda Operator (ALO), calculated in parallel with the solution of the continuum transfer following the method suggested by Rybicki & Hummer (1991, Appendix A).

Substituting this expression into the rate equations, we find for the corresponding effective ionization/recombination rate coefficients

$$R_{i \kappa}\rightarrow\int_{\nu_i}^\infty \frac{4 \pi a_\nu}{h \nu} \Delta J_{\nu}~ {\rm d}\nu$$ (54)

$$R_{\kappa i}\rightarrow\int_{\nu_i}^\infty \frac{4 \pi a_\nu}{h \... ..._{\nu}^{\ast}) + \Delta J_{\nu}\right) {\rm e}^{-h \nu/kT_{\rm e}} ~{\rm d}\nu,$$ (55)

i.e., the problematic, optically thick part of the radiation field has been canceled analytically. Again neglecting stimulated emission (the $\Delta J_{\nu}$-term in the recombination rate coefficient above), approximating $S_i^{n-1} = B_{\nu}/b_i^{n-1}$ with the Planck-function $B_{\nu}$ and NLTE-departure coefficient b_i^n-1, and using the radiation temperature at the threshold, $T_{\rm rad,i}$ along with Seaton parameters $\beta=1, s=2$, we have in analogy to Eq. (14)

$$R_{i \kappa}\rightarrow\frac{8 \pi}{c^2} a_i \nu_i^2 \left\{ W \f... ...\ast}}{b_i^{n-1}} \frac{k T_{\rm e}}{h} {\rm e}^{-h \nu_i/k T_{\rm e}} \right\}$$

$$\hspace*{5mm} := R_{i \kappa}^{n-1} - R_{i \kappa}'$$ (56)

$$R_{\kappa i}\rightarrow\frac{8 \pi}{c^2} \frac{kT_{\rm e}}{h} a_i... ...\rm e}} \left(1 - \Lambda_{\nu}^{\ast}\right) ~:=~ R_{\kappa i}- R_{\kappa i}'.$$ (57)

In those cases where an overlapping continuum is present, i.e., if different transitions contribute to the opacity, the ALO has to be modified according to

$$% \Lambda_{\nu}^{\ast}\rightarrow\beta_i(\nu) \Lambda_{\nu}^{... ...} \quad \beta_i = \frac{\chi_i(\nu)}{\chi_{\rm tot}(\nu)}\cdot$$ (58)

$\chi_i$ is the opacity of the considered transition, $\chi_{\rm tot}$ the total opacity and $\beta_i$ is assumed to be constant between two subsequent iterations (cf. Paper I). The opacities used for the radiative transfer are calculated from their actual Seaton parameters  $(\beta,~ s)$, whereas the uniform values $(\beta=1,~ s=2)$ are applied "only'' to evaluate the approximate ionization/recombination rates.

	Figure 3: Approximate NLTE vs. the exact case: He ionization fractions ( from top to bottom: He III, He II, He I) for pure H/He atmospheric models at $T_{\rm eff}$ = 40 000 K ( left panel: dwarf with $\log g$ = 4.5 and thin wind; right panel: supergiant with $\log g$ = 3.7 and thick wind). Bold: exact solution for helium; dotted: He in approximate NLTE (see text).
Open with DEXTER

Since the Lambda Iteration fails only in the optically thick case, we apply the ALI-scheme exclusively for ground state transitions. Thus, by substituting the effective rate coefficients  $R_{1 \kappa}$ and  $R_{\kappa 1}$into Eqs. (4), (5), we have

$$% \frac{n_\kappa}{n_1} = \left(\frac{n_\kappa}{n_1} \right)^*... ...c{n_i}{n_1} \right)^* \frac{R_{\kappa i}}{R_{\kappa_1}}}\cdot$$ (59)

Using again Eq. (17) and the definitions given in Eq. (44), we finally obtain

$$% \frac{n_\kappa}{n_1} = \left(\frac{n_\kappa}{n_1} \right)^*... ...pa}}} {1 - \zeta_1 ~\frac{R_{\kappa 1}'}{R_{\kappa 1}}}\cdot$$ (60)

In the case of $\Lambda_{\nu}^{\ast}\equiv 0$ (implying $R_{1 \kappa}' = R_{\kappa 1}' = 0$), we immediately recover the original result, Eq. (43), since

$$% \left(\frac{n_\kappa}{n_1} \right)^* \frac{R_{1 \kappa}}{R_... ...r,1}}} W \sqrt{\frac{T_{\rm e}}{\displaystyle {T_{\rm r,1}}}}$$ (61)

by means of Eq. (15). If, on the other hand, the ALO is significant (i.e., close to unity), we find

$$\frac{R_{1 \kappa}'}{R_{1 \kappa}} = \frac{\Lambda_{\nu}^{\ast}}{... ...ystyle {T_{\rm e}}} - \frac{1}{\displaystyle {T_{\rm r,1}}^{n-1}}\right)\right]$$ (62)

$$\frac{R_{\kappa 1}'}{R_{\kappa 1}} = \Lambda_{\nu}^{\ast}.$$ (63)

Thus, the reformulated ALI-scheme collapses to a simple correction of the original Eq. (43) for the ground-state population,

$$\frac{n_\kappa}{n_1} = \frac{n_\kappa}{n_1}(\Lambda_{\nu}^{\ast}\... ...t(\displaystyle {T_{\rm r,1}}^{n-1},~b_1^{n-1}\right), \quad \mbox{with factor}$$ (64)

$$C_{\rm A} = \frac{1 - \displaystyle{\frac{\Lambda_{\nu}^{\ast}}{(... ...le {T_{\rm r,1}}^{n-1}} \right)\right]} {1 - \zeta_1 \Lambda_{\nu}^{\ast}}\cdot$$

The consistency of this scheme is easily proven, because after convergence we would get (cf. Eq. (43))

$$% \frac{1}{b_1} \left(\frac{n_1^*}{n_1}\right) = \frac{n_\kappa}{n_1} \left(\frac{n_1}{n_\kappa} \right)^*$$ =

$$W ~ \frac{\displaystyle {T_{\rm r,1}}}{T_{\rm e}}~ \exp \left[ -\... ...T_{\rm r,1}}} - \frac{1}{T_{\rm e}} \right)\right] \zeta_1 ~(1 + C_{{\rm N}'}),$$ = (65)

so that the "ALO-correction factor'' $C_{\rm A}$ becomes unity. Throughout the iteration the correction factor can take values smaller or larger than unity, leading to a fast and reliable convergence.

   
4.7 Test calculations

In order to check the accuracy of our approximate approach, we will present two different test calculations. The first test aims at an investigation of the methods outlined above, unaffected by additional complications such as line-blocking/blanketing. To this end, we have computed a pure H/He atmosphere at $T_{\rm eff}$ = 40 000 K, for two different sets of parameters: the first model (A4045 with $\log g$ = 4.5) corresponds to a dwarf with thin wind, the second (F4037 with $\log g$ = 3.7) to a supergiant with thick wind.

For both models we have calculated an "exact'' solution as described in Paper I, namely by solving for the H/He occupation numbers from the complete rate equations, with all lines in the CMF and a temperature stratification calculated from NLTE Hopf-parameters. In order to test our approach, we calculated two additional models, with an exact solution for hydrogen only, whereas helium has been treated by means of our approximate approach. (In the standard version of our code, helium is always treated exactly.)

	Figure 4: Approximate NLTE (dotted) vs. the exact case (bold): He departure coefficients for model F4037. Upper panel: He II ground-state departure coefficient. Lower panel: He I triplet and singlet "ground''-states (upper and lower curves, respectively).
Open with DEXTER

Figure 3 shows the very good agreement of the resulting ionization fraction for helium in both cases. The small differences at large optical depths (i.e., for LTE conditions) are due to the different atomic models for helium used in both the exact and the approximate solution. (The data-base applied to the approximate solution comprises a lower number of levels for both He I and He II, so that the partition functions are somewhat smaller than in the exact case, and consequently also the ionization fractions. The occupation numbers of the levels in common are identical though).

The excellent agreement of the He II ground state departure coefficient as a function of depth (Fig. 4, upper panel) is most intriguing. The crucial feature is the depopulation of the He II ground-state close to the sonic point, which is a sophisticated NLTE-effect arising in unified model atmospheres and depends on a delicate balance between the conditions in the He II ground-state, the n=2 ionization edge and the He II Ly$_\alpha$ line (which itself depends on the radiation field at 303 Å and the escape probabilities), cf. Gabler et al. (1989). The comparison between exact and approximate solution shows clearly that our approach, accounting for frequency-dependent radiation temperatures and important line transitions, is actually able to cope with such complicated problems.

In the lower panel of the figure, we have displayed the "ground''-state departure coefficients of He I, for the triplet and singlet system (upper and lower curves, respectively). Although the precision is not as good as for the He II ground-state, He I at 40 000 K is an extremely rare ion, and the major features (depopulation of the singlet ground-state, no depopulation for the triplet ground-state) are reproduced fairly well.

	Figure 5: Approximate NLTE (grey) vs. the results of a solution of the complete rate equations, using the Sobolev line transfer (black): ionization fractions of important metals for model F4037. The ionization stages  III, IV, V (dotted, dashed and dashed-dotted, respectively) are displayed.
Open with DEXTER

The second test investigates the behaviour of the metals. We compare the results from the approximate method with results from an "almost'' exact solution, for model F4037. As we will see in Sect. 7, the introduction of a consistent temperature structure calculated in parallel with the solution of the rate equations forced us to consider the most important elements (in terms of their abundance) in a more precise way than described so far, at least if we calculate the temperature from the electron thermal balance. In this case it is extremely important that the occupation numbers from all excited levels are known to a high precision in order to account for the cooling/heating by bound-bound collisions in a concise way. Unfortunately, this latter constraint cannot be fulfilled by our approximate method, because not all excited levels are considered, and small deviations from the exact solution (which are negligible for the effects of line-blocking, see below) can have disastrous effects on the total cooling/heating rates.

	Figure 6: As Fig. 5, but for model A4045.
Open with DEXTER

Thus, for the most abundant elements the complete set of rate-equations has to be solved for in any case, and this solution (which uses a Sobolev line transfer, cf. Sect. 7) is compared to our approximate one in Fig. 5, for the ionization stages III to  V of some important metals, namely C, O, Si, Ar, Fe and Ni. Note that the comparison includes the effects of line-blocking on the radiation field, where this radiation field has been calculated either from the exact occupation numbers or from the corresponding approximate values. Our comparison demonstrates

• The transition between LTE and NLTE (taking place at $\tau_{\rm R}$ > 2/3 in our approximate approach) is described correctly.
• The approximate treatment works particularly well for elements with complex electronic structure (Ar, Fe, Ni), i.e., our treatment of meta-stable levels is reasonable.
• If there are differences, they occur predominantly in the outer wind.

In almost all considered cases, the principal run of the approximate ionization fractions agrees reasonably or even perfectly with the exact result. The only exception is oxygen where the major/minor stages ( IV/ III) appear reversed in the outer wind (no problems have been found for nitrogen and neon which are not displayed here). These differences in the outer wind (see also C and Si) are partly due to two effects. On the one hand, our approach becomes questionable in those cases when all line transitions are optically thin so that the two-level-atom approach fails to describe the excitation-balance of subordinate levels. If only this effect were responsible this would imply (as suggested by our referee) that the discrepancy should become worse for thinner winds. Thus, we performed a similar comparison for model A4045, which has a considerably lower wind density than model F4037, by a factor of almost 100. The corresponding ionization fractions are shown in Fig. 6. Note that the transition point between photosphere and wind is located at lower values of  $\tau_{\rm R}$, compared to model F4037, due to the weaker wind. Interestingly, the discrepancies between approximated and "exact'' ionization fractions in the outer wind have remained at the same level as for model F4037, and in the case of oxygen the situation is almost perfect now. Consequently, the effect discussed above cannot alone be responsible for the observed discrepancy, and we attribute it to a combination of various factors inherent to our approximative approach.

For our models, however, this is of minor importance, since we are not aiming at a perfect description of the occupation numbers in the outer wind unless we actually need it, i.e., when a consistent temperature structure is derived. In this case, the occupation numbers are calculated exactly anyway.

Different occupation numbers influence the radiation field, which in turn influences the occupation numbers, and so on. This is the second important process which might affect our final approximate solution. Figure 7 compares the emergent fluxes (expressed as radiation temperatures) for the converged models of F4037, calculated by both alternative approaches.

Due to the excellent agreement between the ionization fractions in the line/continuum-forming part of the atmosphere, the fluxes also agree very well. The maximum differences, located between 200 to 400 Å, are of the order of $\pm$1000 K, which translates to a typical difference in population of $\pm$0.15 dex in the outer wind.

Globally, however, the differences in flux are so small that we can consider the two results as equivalent. Thus, the radiation field calculated in parallel with the line-blocking background elements is insensitive to the chosen approach (exact vs. approximate occupation numbers) which primarily differs in the precision (and presence) of subordinate levels.

   
5 Approximate line-blocking

The most time-consuming part of the computation of realistic stellar atmospheres is the calculation of the radiation field, realizing the multitude of overlapping lines with considerable opacity (see also the discussion by Puls & Pauldrach 1990; Pauldrach et al. 2001).

For CMFGEN as well as for the wind-code developed by the Potsdam group (for a recent status report, see Gräfener et al. 2002), this problem has been tackled by performing a comoving-frame solution for the complete EUV/UV range. Obviously, this approach is very time-consuming. A quick calculation shows that the number of frequency points which must be treated is of the order of 900 000, if a range between 200 and 2000 Å and a typical resolution of 0.8 km s^-1 is considered (i.e., ten points covering a thermal width of 8 km s^-1).

In the approach followed by WM-Basic, on the other hand, an observer's frame solution is performed which requires "only'' a few thousand frequency points to be considered. The conservation of work, however, immediately implies that in this case a lot of time has to be spent on the resolution of the resonance zones of the overlapping lines, a problem which is avoided a priori in a CMF calculation.

In order to solve the problem on a minimum time-scale, both a Monte-Carlo solution (e.g., Schaerer & de Koter 1997; Schaerer & Schmutz 1994), and a statistical approach are feasible.

Since the number of metal lines to be treated is very large, the information about the exact position of individual lines inside a (continuum transfer) frequency grid interval becomes less important for obtaining a representative mean background. As shown by Wehrse et al. (1998), the Poisson Point Process is well suited to describe such a line ensemble, particularly because it is very flexible and can be described by relatively few parameters.

The additional introduction of a Generalized Opacity Distribution Function by Baschek et al. (2001) serves two purposes. First, additional analytical insight is given into the effects of the vast amount of blocking lines on the mean opacity in differentially moving media with line overlap. Second, it is a fast tool to derive such mean backgrounds numerically. In particular, it is able to "solve problems that have been inaccessible up to now as e.g. the influence of very many, very weak lines'' (Baschek et al. 2001), and to describe the transition from a static to a moving configuration, since it is equally efficient in both cases.

In our opinion, this approach is very promising, and work adapting and applying the corresponding method is presently under way in our group. Since it will take some time to finalize this approach (the most cumbersome problem is the formulation of consistent emissivities), we have followed a somewhat simplified approach, which relies on similar arguments and has been developed by carefully comparing with results from "exact'' methods, mostly with the model grid calculated with WM-Basic as described by Pauldrach et al. (2001).

	Figure 7: As Fig. 5. Comparison of radiation temperatures of converged models.
Open with DEXTER

Again, the principal idea is to define suitably averaged quantities that represent a mean background and that can be calculated easily and fast. The multitude of lines will be approximated in terms of a pseudo continuum (split into a "true'' absorption and a scattering component), so that the radiative transfer can be performed by means of a standard continuum solution, for relatively few frequency points (see below). Strongest emphasis has been given to the requirement that any integral quantity calculated from the radiation field (such as the photo-integrals) has to give good approximations compared to the exact case, because these quantities (and not the frequential ones) are most decisive for a correct description of the level populations and, in turn, for the blocked radiation field.

5.1 Mean opacities

To this end, we define a "coarse grid'' with spacing $2 N v_{\rm th}^{\rm m}$, where  $v_{\rm th}^{\rm m}$ is a typical thermal velocity (say, of oxygen) including micro-turbulence, and 2 N is an integer of the order of 100. (The reason to define 2N here instead of N will soon become clear.) Under typical conditions, this grid has a resolution of 1000-1500 km s^-1 and is used to calculate appropriate averaged opacities. With respect to a simplified approach, a mean constructed in analogy to the Rosseland mean is perfectly suited, i.e., an average of the inverse of the opacity,

$$% \frac{1}{\langle\chi^{\rm tot}\rangle} =: \frac {\displayst... ...t}}}} {\displaystyle{\int_{2N v_{\rm th}^{\rm m}}{\rm d}\nu}},$$ (66)

since it has the following advantages:

1.

if no lines are present, the pure continuum opacity is recovered;

2.

if one frequency interval is completely filled with non-overlapping, strong lines of equal strength, the average opacity also approaches this value; whereas

3.

in those cases when the interval has "gaps'' in the opacity, these gaps lead to a significant reduction of the mean, i.e., allow for an appropriate escape of photons. Note that any linear average has the effect that one strong line alone (of typical width  $2 v_{\rm th}^{\rm m}$) would give rise to a rather large mean opacity (just a factor of N weaker than in case b) and, thus, would forbid the actual escape;

4.

the average according to Eq. (66) is consistent with the standard Rosseland mean in the lowermost atmosphere (as long as $\partial B_{\nu}/\partial T$ is roughly constant over one interval), i.e., it is consistent with the diffusion approximation applied as a lower boundary condition in the equation of (continuum) radiative transfer.

Because of the large number of contributing lines (typically $5 \times 10^5$(O-type) to 10⁶ (A-type) lines if only ions of significant population are considered), the calculation of this mean has to be fast.

First, assume that any velocity field effects (leading to Doppler-shift induced line overlaps) are insignificant, i.e., assume a thin wind, so that line blocking is essential only in the subsonic regions of the wind. The generalization to approximate line-overlap in the wind will be described later on.

Instead of evaluating the "exact'' profile function, for each line we use a box car profile of width  $2 v_{\rm th}^{\rm m}$. The frequential line opacity is, thus, given by

 

$$\chi_{\rm L}(\nu) = \left\{ \begin{array}{l} \chi_{\rm L}\quad {\... ... \nu ~\le~ \nu_0 - \Delta \nu_{\rm D}\\ 0 \qquad {\rm else} \end{array}\right.$$ (67)

$$\chi_{\rm L}= \frac{1}{2 \Delta \nu_{\rm D}} \frac{\pi e^2}{m_{\r... ...\frac{n_l}{g_l}, \quad \Delta \nu_{\rm D} = \frac{\nu_0 v_{\rm th}^{\rm m}}{c},$$ (68)

where stimulated emission has been neglected again. Due to this definition, at least the frequency-integrated line opacity is correctly recovered. The coarse frequency grid is now divided into N sub-intervals of width $\Delta \nu = 2 \Delta \nu_{\rm D}$. Inside each of these sub-intervals ("channels'') we sum up any line opacity which has an appropriate rest-wavelength. Thus, we account (approximatively) for any intrinsic (i.e., not wind-induced) line overlap. Inside each channel i, we thus have a (total) frequential opacity

$$% \chi_{\nu,i}^{\rm tot} = \sum_{j} \chi_{{\rm L}j}+ \chi_\nu... ...uad \chi_\nu^{\rm cont}= \chi_\nu^{\rm c,true}+ \sigma_{\rm e}$$ (69)

if lines j are located inside channel i and the continuum opacity is assumed to be constant inside each coarse grid interval. $\chi_\nu^{\rm c,true}$ is the contribution by true absorption processes, and  $\sigma_{\rm e}$ the contribution by electron scattering. After replacing the integrals by appropriate sums and since all channels have the same width, the mean opacity (on the coarse grid) is simply given by

$$% \frac{1}{\langle\chi^{\rm tot}\rangle} \approx \frac{\displ... ...rac{1}{N} \sum_{i=1}^N \frac{1} {\chi_{\nu,i}^{\rm tot}}\cdot$$ (70)

For later purposes we split this mean opacity into the contribution from lines and continuum, respectively, where the line-contribution is given by

$$% \langle\chi_{\rm L}\rangle= \frac{\displaystyle{N}}{\displa... ..._{i} \frac{1}{\chi_{\nu,i}^{\rm tot}}}} - \chi_\nu^{\rm cont}$$ (71)

and we have

$$% \langle\chi^{\rm tot}\rangle= \langle\chi_{\rm L}\rangle+ \chi_\nu^{\rm cont}.$$ (72)

Note that both mean opacities, $\langle\chi_{\rm L}\rangle$ and $\langle\chi^{\rm tot}\rangle$, are frequency dependent as a function of coarse grid index. In accordance with our reasoning from above, Eq. (71) implies that

1.

if $\chi_{{\rm L}j}= 0$ for all lines inside one interval, the correct result $\langle\chi_{\rm L}\rangle= 0$ is obtained;

2.

if the same total line opacity $\chi_{\rm L}(\nu)$ is present inside each channel, this value will also be obtained for the mean, $\langle\chi_{\rm L}\rangle= \chi_{\rm L}(\nu)$;

3.

if only one (strong) line is present, the mean line opacity is given by $\langle\chi_{\rm L}\rangle\approx \chi^{\rm cont}/(N-1)$, i.e., it will be much smaller than the continuum opacity, since most of the flux can escape via the (N-1) unblocked channels (according to our present assumption that Doppler-induced line overlap is negligible).

The opacities constructed in this way are used also to calculate the photospheric line pressure, in analogy to the description given in Paper I (Eq. (3)), however including the line contribution (cf. Fig. 11).

5.2 Emissivities

In order to calculate the corresponding emissivities, we assume that each transition can be described by a two-level atom, where the lower occupation number is known from the solution ("exact'' or approximate) of the rate equations.

Although this assumption is hardly justified for (weak) recombination lines, it is a fair representation for most of the stronger transitions arising from either the ground-state or a meta-stable level, particularly if the level population itself is calculated from a multi-level atom.

It might be argued that the two-level atom approach is superfluous for those connecting transitions which are calculated from an exact NLTE solution, since the occupation numbers for both levels and, thus, the source-functions are already known. The maximum number of these lines is of the order of 30 000, and therefore much lower than the total number of lines we are using for our line-blocking calculations (cf. Sect. 3). For the latter transitions, however, only the lower level is present in the atomic models, so that the corresponding source-functions have to be approximated in any case.

Moreover, treating all lines (including the connecting transitions) in a two-level way has the additional advantage that the contribution of scattering and thermal processes can be easily split, which allows us to simulate their impact by means of a pseudo-continuum, so that the standard continuum transfer can be applied without any modification.

To keep things simple and as fast as possible and to be in accordance with our assumption of box car profiles, we replace the scattering integral inside the two-level source-function by mean intensities, i.e., we write

$$% S_{\rm L}= \rho J_\nu + \delta B_\nu, \quad \rho = 1- \delta,$$ (73)

where $\delta$ has been defined in Eq. (27) and is evaluated for the line-specific thermalization parameter and escape probability. The total source-function (in channel i, before averaging) is then given by

$$% S_{\nu,i} = \frac{\eta_\nu^{\rm c,true}+ \sigma_{\rm e}J_\n... ...\delta_j B_\nu)}{\chi_\nu^{\rm cont}+ \sum_j \chi_{{\rm L}j}},$$ (74)

with $\eta_\nu^{\rm c,true}$ being the thermal component of the continuum emissivity. Note that the frequential line-opacity $\chi_{{\rm L}j}$ includes the "profile function'' $(2 \Delta\nu_{\rm D})^{-1}$, cf. Eq. (67).

In the following, we will investigate how to average the above quantities in order to be consistent with our definition of  $\langle\chi^{\rm tot}\rangle$ and  $\langle\chi_{\rm L}\rangle$. With respect to the equation of transfer, which will be finally solved on the coarse grid, we find that after integration over the subgrid-channels

$$% \frac{1}{\langle\chi^{\rm tot}\rangle} \frac{{\rm d}}{{\rm ... ...e I_\nu\rangle = \langle S_\nu\rangle - \langle I_\nu\rangle ,$$ (75)

with z being the depth variable along the impact parameter p in the usual (p,z)-geometry. Strictly speaking, the first term in the above equation (i.e., the mean inverse opacity) is given by

$$% \frac{1}{\langle\chi^{\rm tot}\rangle} = \frac{\int \frac{1... ...nu {\rm d}\nu}{\int \frac{{\rm d}}{{\rm d}z} I_\nu {\rm d}\nu}$$ (76)

(where the denominator is equivalent to ${\rm d}/{\rm d}z \langle I_\nu\rangle$, and all integrals extend over the range  $2 N v_{\rm th}^{\rm m}$), i.e., a different definition applies when compared to the corresponding quantity in Eq. (66). Our crucial approximation is to equate both definitions, i.e., inside each coarse grid cell (of width $\approx$1000...1500 km s^-1) we assume that

$$\begin{eqnarray*}\int \frac{1}{\chi_\nu} \frac{{\rm d}}{{\rm d}z} I_\nu {\rm d}\... ...ox \quad \int \frac{{\rm d}\nu}{\chi_\nu} ~ / ~ \int {\rm d}\nu. \end{eqnarray*}$$

We frankly admit that this approximation can be justified only if a) the spatial gradient of the specific intensity is a slowly varying function of frequency (i.e., deep in the atmosphere); or b) the opacities are similar for most of the sub-channels, i.e., either no lines are present at all or the (summed) line-opacities do not vary much. Most important, this approximation still works in those cases when only few channels are populated by large opacities and the rest is filled by a weak background due to the inverse relation between opacity and intensity: on the lhs, the high opacity channels do not contribute to the fraction because of the correspondingly low intensities in both the nominator and the denominator, whereas on the rhs these channels drop out at least in the nominator because of the low value of $1/\chi$.

There are, of course, a number of cases where the above approximation is poor. With respect to the results presented below and since we are not aiming at a perfect, highly resolved description of the radiation field in the line-blocking EUV/UV regime, the errors introduced by the above approximation (and the following one, which is of similar quality) are acceptable though.

In order to proceed with appropriate expressions for the emissivity, the mean source-function, $\langle S_\nu\rangle$, is given by

$$% \langle S_\nu\rangle \!=\! \frac{\eta_\nu^{\rm c,true}}{\la... ...i_{{\rm L}j}\delta_j}{\chi_\nu^{\rm tot}}\right\rangle ~B_\nu,$$ (77)

where the Planck-function $B_{\nu}$ is assumed to be constant within one macro-grid interval. For those averages which are multiplied by $\langle J_\nu\rangle$, we have employed an approximation similar to the one discussed above. If we denote the opacity-dependent means of the third and fourth term by  $f_{\rm nth}$ and  $f_{\rm th}$, respectively (i.e., non-thermal/thermal), the equation of radiative transfer for the averaged quantities becomes

$$\frac{1}{\langle\chi^{\rm tot}\rangle} \frac{{\rm d}}{{\rm d}z} \... ...\ranglef_{\rm th}~B_\nu}{\langle\chi^{\rm tot}\rangle} - \langle I_\nu\rangle ,$$ (78)

and can be solved in the conventional way (pure continuum transport). The resulting quantities for the radiation field are to be understood as average quantities, in the sense that integral quantities such as $\int J_\nu {\rm d}\nu$ or $\int H_\nu {\rm d}\nu$ are described correctly, at least in most cases. The coefficients  $f_{\rm th}$ and  $f_{\rm nth}$ can be calculated by summing over the sub-channels,

$$f_{\rm th}= \frac{1}{N} \sum_i \frac{(\sum_j \chi_{{\rm L}j}\delta_j)_i}{(\chi_\nu^{\rm cont}+ \sum_j \chi_{{\rm L}j})_i}$$ (79)

$$f = f_{\rm th}+ f_{\rm nth}= \frac{1}{N} \sum_i \frac{(\sum_j \chi_{{\rm L}j})_i}{(\chi_\nu^{\rm cont}+ \sum_j \chi_{{\rm L}j})_i} \quad < 1$$ (80)

$$f_{\rm nth}= f - f_{\rm th},$$ (81)

and after some simple algebraic manipulations

$$% \langle\chi_{\rm L}\rangle= \langle\chi^{\rm tot}\rangle~(f_{\rm th}+ f_{\rm nth}).$$ (82)

With this equation it is easy to show that the mean source-function (77) allows for a correct thermalization, if $\eta_\nu^{\rm c,true}\rightarrow\chi_\nu^{\rm c,true}B_\nu$ and $\langle J_\nu\rangle \rightarrow B_\nu$. In this case, the mean source-function becomes $\langle S_\nu\rangle = B_\nu$, q.e.d.

	Figure 8: FASTWIND vs. WM-Basic (grey): comparison of emergent fluxes for two dwarf and two supergiant models at 35 and 45 kK (for parameters, cf. Pauldrach et al. 2001). In order to allow for a meaningful comparison, the high resolution frequency grid provided by WM-Basic has been re-mapped while keeping the corresponding flux-integrals conserved.
Open with DEXTER

We now need to incorporate the effects of the velocity field into our approach. Due to the method to average the opacity, we cannot simply shift the lines with respect to the stellar frame. Consider, e.g., one strong line to be present without any other interfering lines. In "reality'' and in the observer's frame, the absorption part of this line becomes broader as a function of velocity, i.e., the larger the velocity the more flux is blocked (of course, a significant part is reemitted due to scattering). If we simply shift our line(s) as a function of velocity, almost nothing would happen, since, as shown above, the mean opacity/radiation field remains almost unaffected by one strong line, due to the possible escape via the (N-1) unblocked sub-channels. Thus, in order to simulate the physical process, we proceed in a different way. When the velocity shift becomes larger than twice the average "thermal'' width (including micro-turbulence), we combine (in proportion to the local velocity) more and more subchannels to increase the relative weight of the line in the mean opacity. In particular, the line width (more precisely, the width of the sub-channels) is set to the value

$$% \Delta \nu = {\rm max}\left(2 \Delta \nu_{\rm D}, \frac{{v}(r)}{\lambda_0}\right)\cdot$$ (83)

Although this procedure is highly approximative, it allows us to deal with the effects of "line-shadowing'' and prevents any premature escape of photons when the lines begin to overlap.

Figure 9: Comparison of ionizing photon number for the model grid provided by Pauldrach et al. (2001). Left panel: logarithm of Zanstra-integrals, $\log Q_x$, (Eq. (84)), for H I, He I, O II and He II. Bold/dashed: dwarfs/supergiants as calculated by WM-Basic; dotted/dashed-dotted: results from FASTWIND. Right panel: ratio of corresponding Zanstra integrals, $\Delta \log Q_x = \log Q_x^{\rm WMB} - \log Q_x^{\rm FW}$ (WMB: WM-Basic, FW: FASTWIND), for dwarfs (triangles) and supergiants. For the supergiant model at 30 000 K ("S30''), we have used the FASTWIND model without photospheric line-pressure in both figures. The corresponding results for the "correct'' model, i.e., including photospheric line-pressure, are indicated by the "+'' sign (see text).

Open with DEXTER

   
5.3 Tests/comparison with wm-Basic

Before we test our approximate approach by comparing with alternative calculations, we mention two important consistency checks we have performed.

a)

The calculated models (and spectral energy distributions/line profiles) are almost independent of the actual value of coarse grid cells, N, at least if varied within a reasonable range. (We checked for values between 0.5 Nto 2 N, for N=60.)

b)

As long as the IR/radio-range is not considered, our simpler models with a temperature-structure calculated from Hopf-parameters and all background elements in approximate NLTE agree very well with complex models including a consistent T-structure. This check verifies analyses performed with previous versions of FASTWIND, e.g., Herrero et al. (2002) and Repolust et al. (2004).

In the following, we will compare the fluxes from our models with those calculated by WM-Basic by means of the O-star grid presented by Pauldrach et al. (2001). These tests should give a reasonable agreement, since both codes use the same atomic data base for the background-elements. A comparison with results from CMFGEN will be discussed later on.

The parameters of the corresponding models (calculated without X-rays) can be found in Pauldrach et al. (2001, Table 5). Our models have been constructed as closely as possible to the approach inherent to WM-Basic, i.e., including a consistent temperature stratification (which will be described in Sect. 7) and Sobolev line-transfer. For the velocity-field, we have used $\beta = 0.9$, which results in a stratifiction very close to the one predicted by WM-Basic (see below). The computation time on a 2 GHz processor machine is of the order of 15 to 20 min per model (typically 40 to 50 iteration cycles for a final convergence below 0.003 in all quantities, if the temperature is updated each 2nd cycle).

The grid comprises 6 "dwarfs'' and 5 "supergiants'' in the range between 30...55 kK ("D30''..."D55'' and "S30''..."S50'', respectively), and we have concentrated on the grid with solar abundances, in order to deal with more prominent effects related to line-blocking/blanketing. Figure 8 compares the emergent fluxes for some typical cases, two dwarf and two supergiant models at 35 and 45 kK. In order to allow for a meaningful comparison, we have re-mapped the high resolution frequency grid provided by WM-Basic while keeping the corresponding flux-integrals conserved.

	Figure 10: Comparison of velocity/density structure for model S30. Grey: WM-Basic; bold/dashed: FASTWIND with/without photospheric line-pressure, respectively.
Open with DEXTER

Overall, the agreement is rather good; in particular, the range above 400 Å is reproduced very well, except for some strong absorption/re-emission features which are missing in our mean-opacity approach. We have convinced ourselves that in all cases the IR-fluxes (not displayed here) also agree perfectly, i.e., the IR flux-excess induced by the wind is reproduced equally well in both codes. Major differences are "only'' present in two regions: most models differ in fluxes below 200 Å, although the strength of He II Lyman-jump itself is very similar. Mostly, this problem is related to the enormous bound-free opacity provided by O IV (and Fe V or C IV for the hotter or cooler objects, respectively) leading to an optically thick wind from the outermost radius point on (in our case, $R_{\rm max}= 120~R_\star$), so that the flux is rather badly defined in this frequency region. As we will see from a comparison with models calculated with CMFGEN (Fig. 15), these models predict a third alternative for $\lambda < 200$ Å, and even the Lyman-jump is different. As a result, we consider the ionizing fluxes in this wavelength range as not particularly reliable. Moreover, the influence of X-rays becomes decisive, implying that any tool for nebula diagnostics should use these numbers with care.

The second inconsistency is found in the region between 300 to 400 Å. Although this range poses no problem for supergiants, the flux-blocking predicted by FASTWIND for dwarfs between 35 to 45 kK is larger than calculated by WM-Basic, with a maximum discrepancy around 35 kK. CMFGEN again produces somewhat different results in this range: agreement with WM-Basic is found for dwarfs, whereas the fluxes emitted from supergiants are larger compared to both FASTWIND and WM-Basic.

This dilemma becomes particularly obvious if we consider the corresponding Zanstra-integrals,

$$% Q_x = \int_{\nu_x}^\infty \frac{H_\nu}{h\nu} {\rm d}\nu,$$ (84)

which are proportional to the emitted number of ionizing photons. In the left panel of Fig. 9, we compare the logarithm of Q_x, evaluated for H I, He I, O II and He II, whereas the corresponding ratios, $\Delta \log Q_x = \log Q_x^{\rm WMB} - \log Q_x^{\rm FW}$ (WMB: WM-Basic, FW: FASTWIND), are displayed in the right panel. Obviously both codes predict the same numbers in the hydrogen Lyman and in the He I continuum. As already discussed, the situation is much less satisfactory for the He II continuum, where the differences are particularly significant for supergiants. Note, however, that the principal dependence of  $Q_{\rm HeII}$ on spectral type and luminosity class, which shows the largest variation throughout the spectrum (lower left panel), is much more consistent than one might expect on basis of the right panel alone. In the O II continuum ($\lambda <$ 352 Å), finally, the differences for the dwarfs at intermediate spectral type are evident.

Note that in this wavelength range the line-density is very large, and differences in the treatment of the weakest background opacities might explain the established disagreement. An argument in support of this hypothesis is given by the fact that FASTWIND recovers the results by WM-Basic perfectly if a line-list is used which has significantly fewer (overlapping) weaker lines in the considered interval. For a final statement, however, more tests are certainly required. Note that a comparison with CMFGEN addressing this point will not solve the problem, since the number of lines included in this code is mostly lower than described here, because CMFGEN uses only those lines where the occupation numbers of both levels are known, in contrast to our approach which uses also lines where the upper level is lying too high to be included into the rate equations.

One last point we would like to mention concerns model S30. In a first comparison, we immediately encountered the problem that this model provided fluxes which showed significantly less agreement at all frequencies than the other models (indicated by the plus-signs in Fig. 9). Comparing the models themselves, it turned out that temperature, density and velocity structure showed a severe mismatch in photospheric regions (cf. Fig. 10, grey vs. black curves). After some tests, we found that both models agree well if the photospheric line pressure is neglected in FASTWIND (grey vs. dashed curves in Fig. 10). Most likely, this problem is related to the treatment of the line pressure in WM-Basic. Whereas the continuum forces are calculated from correctly evaluated opacities, the line pressure, independent of location, is calculated in terms of the force-multiplier concept, utilizing the Sobolev approximation. Particularly, $g_{\rm rad}^{\rm line} \propto t^{-\alpha}$, with "depth parameter'' $t\propto \rho/({\rm d}v/{\rm d}r)$. Thus, $g_{\rm rad}$ decreases rapidly in photospheric regions when the density is large and the velocity gradient small.

In those cases where the (static) line pressure is non-negligible in photospheric regions, the chances are high that the above approximation leads to a too large effective gravity, i.e., too high densities. Actually, this problem has been known for a long time and has been discussed in detail in Pauldrach et al. (1986, particularly Fig. 6c#. The reason that this problem occurs only in S30 is that the Eddington factor is considerably higher than for almost all other models ( $\Gamma = 0.52$). Thus, the photospheric line pressure has much more impact than for models with either high gravity or low $\Gamma$. Moreover, at an effective temperature of 30 kK, Fe IV with its enormous number of lines spread throughout the spectrum is the dominant (or almost dominant) ionization stage in the "middle'' photosphere, thus contributing a much larger amount of static line pressure than for hotter temperatures, where Fe V or even Fe VI are contributing.

We have also compared our (cooler) models (from our grids as described in Sect. 9 and from additional A-star models) with corresponding Kurucz models, where in most cases a very good agreement regarding the photospheric radiative acceleration has been found, e.g. Fig. 11. Only for models cooler than 9000 K does a mismatch become obvious, where "our'' radiation pressure is too low, due to a number of missing Fe II lines in the optical (improvements are under way).

Figure 11: Comparison of total photospheric radiative acceleration for model S30 (bold) vs. results from an analogous hydrostatic Kurucz-model (dotted). Note that the gravitational acceleration for this model is 1000 cm s-2, i.e., the radiative acceleration is very close to this value and, thus, of extreme importance (cf. Fig. 10). The deviations at largest depths are due to the fact that this model becomes (spherically) extended in the lowermost photosphere, an effect which cannot be treated in a plane-parallel approach (cf. Paper I).

Open with DEXTER

For a meaningful comparison concerning our approximate line-blocking, in Fig. 9 we have used the results from our S30 model without photospheric line-force, whereas the results from the "actual'' model (including $g_{\rm rad}^{\rm line}$) are indicated by "+''. Independently, however, Fig. 10 (left panel) also shows the validity of our treatment of the transition zone from photosphere to wind (cf. Paper I), since in this region both velocity fields agree perfectly. (Remember that WM-Basic solves the hydrodynamical equations in a consistent way.)

   
6 Treatment of inverted levels

One of the more complex problems when solving the coupled equations of statistical equilibrium and radiative transfer is the presence of population inversions, which often occur in the outermost layers of hot expanding stellar atmospheres. The amount of the overpopulation (i.e., n_u/g_u > n_l/g_l) is usually small, but even in this case it invokes a number of problems concerning the solution of the radiative transfer equation. Particularly with respect to the usual concept of using source functions, a problem occurs in the transition zone between "normal'' population and overpopulation, where the source function formally diverges. In addition, factors like $\exp(-\tau)$ may produce numerical problems for $\tau<0$. In a number of codes, this problem is "solved'' by setting the upper level into LTE with respect to the lower one or by other approximations. Since level inversions are particularly present between levels responsible for IR-lines and since FASTWIND aims at a reliable solution also in these cases, we cannot afford such approximations and have to solve the "exact'' case which in turn has an influence on the degree of overpopulation itself. In this section, we briefly describe how we have solved the problem in FASTWIND both with respect to the Sobolev approach and within the CMF-transport.

6.1 Treatment of inversions in the Sobolev approximation

Since the Sobolev approach uses only local quantities, a divergence of the source function is not possible, except for the extremely unlikely case that upper and lower populations, normalized to the appropriate statistical weight, are numerically identical. Thus, we can retain the standard concept (optical depth and source function) and follow the approach described in Taresch et al. (1997): in the case of a level inversion, the interaction function $U(\tau_{\rm S},\beta_{\rm P})$ is split into two parts in order to avoid numerical problems,

$$% \bar{U}=U_{1}+U_{2}, \qquad {\rm with} \qquad U_{1}=1-\beta.$$ (85)

$\beta$ is the usual escape probability in Sobolev approximation (Eq. (29)), which for the case of inversion is given by

$$% \beta=\frac{\exp{\vert\tau_{\rm S}\vert}-1}{\vert\tau_{\rm S}\vert}$$ (86)

and U₂ has been described in Taresch et al. (1997, Eq. (A13)).

For $\vert\beta_{\rm P}\vert \rightarrow \infty$, i.e., dominating continuum, U₂ approaches zero. In the case of dominant line processes, on the other hand, and $\tau_{\rm S} < 0$, U₂ approaches ($\beta-1)$ and $\bar{U}$ goes to zero. Thus, we recover the "classical'' result by Sobolev, where the influence of continua has been neglected.

In our approach, we have significantly extended the grid used by Taresch et al. (1997) from which U₂ is calculated by means of interpolation. Due to the different behaviour of this function in different regions of the $(\tau_{\rm S}, \beta_{\rm P})$ plane, (four) different tables with different degrees of resolution have been calculated. The boundaries of the complete grid comprise the area between $-6 \le \log\vert\tau_{\rm S}\vert < 2.8$ and $-6 \le \log\vert\beta_{\rm P}\vert < 6$. Beyond the boundaries, U₂ is calculated analytically (by either considering the appropriate limits or using a first order expansion). In particular,

$$U_2 = 0.5 \vert\tau_{\rm S}\vert, \qquad \log \vert\tau_{\rm S}\vert < -6 \quad \mbox {independent of } \beta_{\rm P}$$ (87)

$$U_2 = \beta-1, \qquad \log \vert\beta_{\rm P}\vert < -6 \quad \mbox {independent of } \tau_{\rm S}$$ (88)

$$U_2 = \frac{1}{\sqrt{2\pi\vert\beta_{\rm P}\vert}}, \quad \log \vert\beta_{\rm P}\vert > 6, \quad \log \vert\tau_{\rm S}\vert > -6$$ (89)

and the limits for large $\log \vert\tau_{\rm S}\vert >2.8$ and $-6 \le \log\vert\beta_{\rm P}\vert \le 6$ result from appropriate extrapolations from the pre-calculated tables.

6.2 Treatment of inversions in the comoving frame

In the CMF solution, the problem of source-function divergence is inevitable when a population inversion occurs and the standard formalism is used. Even if the local quantities are not diverging, there will be an implicit divergence between the two depth-points before and at the beginning of overpopulation, which, due to the applied discretization, will not be handled consistently. To avoid this problem, it is more suitable to work directly with emissivities and opacities rather than with optical depths and source functions. Thus, in the case of inversion, we solve the two coupled equations of radiative transfer in the comoving frame according to

$$\frac{\partial u}{\partial z}-\frac{\partial v}{\partial x} = -\chi v$$ (90)

$$\frac{\partial v}{\partial z}-\frac{\partial u}{\partial x} = \eta-\chi u,$$ (91)

where u and v are the usual Feautrier variables, x is the frequency measured from the center of the line in Doppler width units, and z is the depth variable along the impact parameter. The opacity is $\chi = \chi_{c}(z) + \chi_{L}(z,x)$ and the emissivity is $\eta=\eta(z,x)$.

In order to discretezise the equations with respect to z and x, a fully implicit scheme is used. As was shown by Mihalas et al. (1975, Appendix B) this method is unconditionally stable.

6.3 Tests

A number of tests have been performed concerning both the Sobolev and the CMF implementations. Most importantly, we have also tested models where the above discretization of the CMF equations with respect to z has been used for all transitions, not only for the "inverted'' ones, and found satisfactory agreement with our standard implementation using a discretization with respect to $\tau$.

After convincing ourselves that the algorithms are working in principal, we have tested our improved methods by comparing them with older results (where in case of inversion the upper level and the line source function were set to zero). This comparison has been performed for the O-star grid described in the previous section. The results were very satisfying, and a number of convergence problems originating from the older treatment of inverted populations are no longer present.

The differences in the resulting H/He line profiles (both in the optical and in the IR) turned out to be rather small, since for our grid parameters these lines are formed below those regions where the inversion sets in. However, a consistent treatment might be important for winds with more extreme mass-loss rates and for a number of metallic IR transitions with an inversion already occurring in photospheric regions.

   
7 Temperature stratification

As has been previously mentioned, the present version of FASTWIND allows for the calculation of a consistent temperature stratification, utilizing a flux-correction method in the lower wind and the thermal balance of electrons (cf. Kubát et al. 1999) in the outer part. The region where both methods are connected is somewhat dependent on mass-loss, but typically lies at $\tau_{\rm R}$ = 0.5. Although the implementation of this method is straightforward, and the contribution of individual processes have been discussed in detail by Drew (1985,1989), three points are worth mentioning.

In order to calculate the appropriate heating/cooling rates resulting from collisional bound-bound transitions, the population of excited levels is as important as the population of ground and meta-stable ones. This can readily be seen from the fact that the net heating rate from a collisional transition between lower level l and upper level u can be expressed as

$$% Q_{ul}-Q_{lu} = \left(n_u C_{ul} - n_l C_{lu}\right) h\nu_{lu} = n_l C_{lu} h \nu_{lu} \left(\frac{b_u}{b_l} -1\right),$$ (92)

with collisional rates C_ul and NLTE departure coefficients b_l, b_u. Thus, the ratio of departure coefficients controls whether a certain transition heats or cools the plasma and its deviation from unity controls the degree of energy transfer. Heating results from transitions with an upper level being overpopulated with respect to the lower one, and cooling vice versa. Thus, the occupation numbers of all ionic levels have to be known with some precision, and we have to modify our approach when the electron thermal balance is used to calculate the temperature profile. The approximate NLTE solution as described in Sect. 4 does not yield the required occupation numbers of excited levels (except those which are directly connected to the ground or meta-stable level), and any brute force approximation would give incorrect heating/cooling rates.

To overcome this dilemma we incorporated a detailed solution of the statistical equilibrium at least for those elements with large contributions to the net heating rates (positive or negative). After some experiments it turned out that the inclusion of the most abundant background elements C, N, O, Ne, Mg, Si, S, Ar, Fe, Ni (plus the explicit elements, of course) is sufficient to stabilize the results. For these elements then, the complete rate-equations are solved with line transitions treated in Sobolev approximation, whereas for the remaining ones the approximate NLTE solution is employed.

The second point to be mentioned regards the flux-conservation of the final models. The conventional approach to calculate the energy balance, formulated in terms of radiative equilibrium, satisfies this constraint by construction, at least in principle. (Most numerical codes, including CMFGEN and FASTWIND, calculate mean intensity and flux on different grids, which somewhat destroys the coupling between radiative equilibrium and flux conservation). On the other hand, our formulation in terms of the electron thermal balance is decoupled from the latter requirement, at least regarding any explicit dependence. Note, however, that there is an implicit coupling via the rate equations, assuring that the constraints of electron thermal balance and radiative equilibrium are physically equivalent (cf. Hillier & Miller 1998; Hillier 2003, where further discussion concerning both methods and their correspondance is given). Thus, we can use the achieved flux-conservation as an almost independent tool to check whether our models have been constructed in a consistent way. In most of the cases considered so far we have found a perfect conservation, but in the worst cases (below 5% of all models) a violation up to 1.5% is possible.

The third point to be discussed is mainly relevant for our specific approach of modeling stellar atmospheres. Presently, and in accordance with the majority of similar codes, we do not update the photospheric density stratification once it has been calculated. Since the photospheric structure equations are solved for the gas-pressure P and the density is calculated from the ratio P/T, the density is only as good as the initial "guess'' for the temperature stratification. Moreover, an implicit dependence of the final temperature distribution on this initial guess is created.

Thus, it is still important to obtain a fair approximation for the latter quantity, which in our models is accomplished via the corresponding NLTE Hopf-parameters (see Paper I) which have to account for line-blanketing effects. Meanwhile, we have accumulated a large set of these parameters from our model-grid calculations (and, for cooler temperatures, from corresponding Kurucz-models). If, on the other hand, the initial (photospheric) temperature stratification were not appropriate, both occupation numbers and line profiles would be affected by the erroneous density (although the flux would be conserved, see above).

In Fig. 12 we show some of our results in comparison with results calculated by means of WM-Basic, a code that also uses the electron thermal balance. Obviously, the differences are tiny and visible only for the temperature bumps of supergiants, which are predicted by WM-Basic to be more prominent. Note, however, that our solution is more consistent with the results from CMFGEN (see Fig. 14), which will be presented in the next section.

Comparing the computation time of models with and without consistent temperature structure, we find a typical difference of a factor of two. Interestingly, the number of iterations becomes only moderately larger (because of the fast convergence of the temperature when using the electron thermal balance, see Kubát et al. 1999), and most of the additional time is spent solving the NLTE equations for the important back-ground elements.

	Figure 12: FASTWIND (bold) vs. WM-Basic (grey): comparison of temperature stratification for some of the models described in Sect. 5.3.
Open with DEXTER

	Figure 13: Extreme temperature-bump around 22 000 K: FASTWIND (bold) vs. WM-Basic (grey, dashed) and CMFGEN (crosses). See text.
Open with DEXTER

We finish this section with an interesting finding and warning. After having calculated a large number of models with our code, in certain domains of $T_{\rm eff}$ we have found temperature bumps of extreme extent. In contrast to "normal'' bumps (arising from line-heating in the outer photosphere) which are of the order of 2000 K or less for O-stars (Fig. 12), corresponding values at lower effective temperatures might reach 5000 K, as shown for an exemplary dwarf-model at $T_{\rm eff}$ = 22 400 K in Fig. 13.

This behaviour has been confirmed by calculations performed by WM-Basic and CMFGEN, kindly provided by T. Hoffmann and F. Najarro. This finding allows for two conclusions. First, the effect is "real'', at least in terms of the applied physics (see below), and second, the results from different codes using different techniques are strongly converging, which increases our confidence in the results.

After some investigation, it turned out that the feature under discussion originates from bound-bound heating by C III(which is a major ion at these temperatures), contributed by few transitions connected to the ground-state (singlet), to the meta-stable level (lowermost triplet state) and the transition between ground and meta-stable level at roughly 1909 Å. Note that the latter transition has been identified to be of significant importance for the energy-balance in the wind of P Cyg, in that case as a cooling agent (cf. Drew 1985, Fig. 3). In our case, however, the C III ground-state is strongly underpopulated in the transonic region (because of the same effect under-populating the He II ground-state in hot stars, cf. Gabler et al. 1989), so that the bracket in Eq. (92) becomes very large and the heating-rate enormous, also because of the large collisional strengths of these transitions. If, on the other hand, the contributions by C III are neglected, a temperature bump of only moderate size is created.

The lesson we learn from this exercise is two-fold. First, only a few lines (from one ion) can lead to a considerable heating in stellar atmospheres, at least theoretically. Since this heating takes place in the outer photosphere it will have a significant effect on the spectra, and we can check this prediction observationally. However, we have also to consider that the degree of heating (i.e., the extent of the temperature bump) depends strongly on the corresponding collision strengths of the responsible transitions (as a function of temperature), and before relying on our results we have to carefully check for possible uncertainties.

   
8 Comparison with CMFGEN

In this section, we will compare the results from our models with corresponding results from CMFGEN, with particular emphasis on the optical H/He profiles which cannot be compared to results from WM-Basic, due to lack of comoving frame transport and adequate line-broadening. For this purpose, we have used the CMFGEN-simulations by Lenorzer et al. (2004), who have provided a grid of dwarf, giant and supergiant models (no clumping) in the O-/early B-star range. The corresponding FASTWIND models have been calculated with identical parameters, and the explicit elements (H/He) have been treated with comoving frame transport. Thus, the only "physical'' difference in both calculations concerns the photospheric density stratification, which is approximated by a constant scale-height in CMFGEN, but described consistently by FASTWIND (cf. Sect. 5).

	Figure 14: As Fig. 12, but for FASTWIND (bold) vs. CMFGEN (grey, dashed). The stellar parameters are similar to the models displayed in Fig. 12, with $T_{\rm eff}$(6V) = 35 861 K, $T_{\rm eff}$(3V) = 43 511 K, $T_{\rm eff}$(5Ia) = 35 673 K and $T_{\rm eff}$(2Ia) = 44 642 K. CMFGEN results from the model grid as calculated by Lenorzer et al. (2004).
Open with DEXTER

	Figure 15: As Fig. 8, but for FASTWIND vs. CMFGEN (grey). Effective temperatures as in Fig. 14. Only the EUV part is plotted, at larger wavelength the results are extremely similar.
Open with DEXTER

Figure 16: As in Fig. 9, right panel, but for CMFGEN vs. FASTWIND (positive values result from Zanstra-integrals being larger in CMFGEN). Triangles: dwarfs; asterisks: supergiants. Note that the x-axis extends only until $T_{\rm eff}$ = 50 000 K. The three objects denoted with arrows in the lowermost panel (He II) correspond to the dwarf models "3V'' and "4V'' and the supergiant model "1Ia'', respectively. The resulting differences for $\Delta \log Q_{\rm HeII}$ are given beneath the arrows.

Open with DEXTER

The corresponding temperature profiles are displayed in Fig. 14, for two dwarf and two supergiant models with parameters similar to our comparison with WM-Basic. Remember that the temperature structure is derived from radiative equilibrium in CMFGEN, whereas FASTWIND uses the thermal balance of electrons in the outer atmosphere. Overall, the differences are small, and the extent of the temperature bumps are comparable. The only disagreement is found in the outer wind, where FASTWIND uses an artificial cut-off ( $T_{\rm min}=$ 0.4  $T_{\rm eff}$) in order to prevent numerical problems at lower effective temperatures. We have convinced ourselves that this cut-off has no further consequences for the models as described here, which neglect adiabatic cooling in the outer wind.

Figure 15 compares the corresponding EUV-fluxes, in analogy to Fig. 8. As already discussed in Sect. 5.3, the largest differences occur in the He II-continua. This effect can be seen even more clearly in Fig. 16, lowest panel. Regarding the supergiants, the deviation is contrary to our comparison with WM-Basic. The WM-Basic He II-fluxes were mostly lower than those from FASTWIND, whereas the CMFGEN-fluxes are larger, particularly at the edges, so that the corresponding Zanstra integrals become larger as well. Thus, the FASTWIND results for  $Q_{\rm HeII}$ lie roughly in the middle of the results from CMFGEN and WM-Basic, at least for the supergiants. Again, note the extreme sensitivity of the model predictions in this frequency range. The reader is warned about any uncritical use of corresponding results, e.g., with respect to nebula modeling.

Regarding the dwarf models, both codes give more or less identical results for the He II-continua for $T_{\rm eff}$ < 36 000 K, whereas at hotter temperatures extreme differences are found for the two models at $T_{\rm eff}$ = 41 000 K and 43 500 K, respectively. In contrast both to our predictions and those from WM-Basic, the CMFGEN-models do not show any He II-edge at all, cf. Fig. 15, model "3V''.

Concerning the O II-continua (actually, for the complete range within 300 Å  $< \lambda <$ 400 Å), the hotter models ( $T_{\rm eff}$ > 35 000 K) show a higher flux-level in CMFGEN, for both the supergiants and the dwarfs. We have already commented on this problem in Sect. 5.3 and speculated that this behaviour is related to missing line-opacity. (O II itself plays no role at these temperatures.) Of course, we cannot exclude a problem in our approximate treatment of line-blocking. In accordance with the comparison with WM-Basic, the agreement of the H I- and He I-continua is almost perfect.

Figures 17 and 19 display the strategic H/He lines in the optical ( CMFGEN-profiles in magenta). Regarding the dwarfs, the agreement of almost all lines is excellent. The only differences are found for the line cores of He II 4686, which are shallower in CMFGEN at almost all temperatures, and for the He I singlets for models "4V'' to "6V'' with $T_{\rm eff}$ lying in the range between 41 000 K and 36 000 K, respectively. (Note that for model "4V'' He I 4387 agrees well whereas He I 4922 and He I 6678 differ.) Most prominent are the differences for models "5V'' and "6V'' (the same is true for the giant models not displayed here), where all singlet lines predicted by CMFGEN are almost a factor of two smaller in equivalent width than those predicted by FASTWIND. Most interestingly, however, the triplet lines agree perfectly throughout the grid.

So far, the origin of this discrepancy could not be identified; particularly, the atomic data used (incl. broadening functions) are very similar, and also the ionizing continua (important for the singlet-formation) agree very well, as shown above. One might speculate that there is a connection to the flux differences around the He II resonance line at 304 Å or to possible discrepancies at the He I resonance line(s), but this has to be checked carefully (investigations under way). Further comments on this discrepancy will be given after we have discussed the results for the supergiants.

The corresponding profiles are displayed in Fig. 19, upper panel. There, the situation is somewhat different to the dwarf case. The deviations of the He I singlets are not as extreme as before. Significant disagreement is found only for He I 4922 and 6678 (no problem for He I 4387) in model "5Ia'' (36 000 K), where these singlets are weak anyway. For model "6Ia'' the differences are moderate, much less than the factor of two in equivalent width encountered above. Noticeable differences are found for other lines though. At first, the hydrogen Balmer line wings predicted by CMFGEN are much stronger, which would lead to lower gravities if an analysis of observed spectra were performed. Second, both H$\alpha$and He II 4686 show stronger wind emission which would lead to lower mass-loss rates compared to FASTWIND. Note however that the wind emission in both lines is a strongly increasing function of mass-loss (e.g., Puls et al. 1996), and an analysis of observed spectra would result in $\dot M$-differences not exceeding the 20 to 30% level.

The difference in the Balmer line wings points to a problem mentioned above, namely the assumption of a constant photospheric scale height in CMFGEN. In order to obtain an impression of how far this approximation (as well as the somewhat artificial transition from photosphere to wind) has an influence on the resulting models and profiles, Lenorzer et al. (2004) have calculated an additional set of "low-gravity'' supergiants, where the gravity has been lowered by typically 0.1 to 0.2 dex (model series "_lg'') with respect to their "standard'' grid of supergiants. Due to this manipulation, at least part of the effect of photospheric radiation pressure  $g_{\rm rad}$ is accounted for (although this quantity is not constant throughout the photosphere), since the profiles provide a measure of the effective gravity (i.e., $g_{\rm grav}-g_{\rm rad}$) alone.

Figure 17: FASTWIND (black) vs. CMFGEN (magenta): comparison of strategic H/He lines in the optical for the dwarf-models from the grid by Lenorzer et al. (2004). For both models, the lines have been degraded to a resolution of 10 000 and rotationally broadened with $v\sin i=80$ km s-1. He I  $\lambda \lambda$ 4387, 4922 and 6678 are singlet lines, and He I  $\lambda \lambda$ 4471 and 4713 are triplets. The horizontal and vertical lines in the bottom right corner indicate the scale used and correspond to 20 Å in wavelength and 0.5 in units of the continuum, respectively (extending from 0.65 to 1.15).

Open with DEXTER

Figure 18: Wind-strength parameter Q as an optical depth invariant: H/He profiles for the model of $\alpha$ Cam as determined by Repolust et al. (2004), with $\dot M$ = 6.04 $\times$ $10^{-6}~M_{\odot}/{\rm yr}$ and $R_\star$ = 32.5 $M_\odot$. Overplotted in magenta are the corresponding profiles for a model with identical Q-parameter (Eq. (93)) but different mass-loss rate and radius ($\dot M$ = 3.3 $\times$ $10^{-6}~M_{\odot}/{\rm yr}$ and $R_\star$ = 21.7 $M_\odot$).

Open with DEXTER

In Fig. 19, lower panel, we compare the FASTWIND profiles (identical to those from the upper panel, since our "high gravity'' models do include the photospheric  $g_{\rm rad}$) with these low-gravity models calculated by CMFGEN. Consequently, the photospheric densities should be much more similar than in the previous case, at least in those regions where the Balmer line wings are formed. Indeed, the differences in H$_\gamma$ and H$_\beta$ have now vanished, and the H$\alpha$ emission is also very similar, except for the hottest models on the blue side of the profile. In some cases, the discrepancy for He II 4686 has become weaker as well. The He I triplets have not changed (they seem to be almost independent of the photospheric density in CMFGEN), whereas a strong influence on the He I singlets is found. In the "critical'' temperature region, they have become significantly weaker, and a strong discrepancy also for the low-gravity model "6Ia_lg'' is present again, by the same degree as we have found for the dwarfs.

	Figure 19: As Fig. 17, however for the supergiants from the grid by Lenorzer et al. (2004). Upper panel: CMFGEN models with "standard'' gravities. Note that the differences in the wings of the Balmer lines and in the H$\alpha$-emission almost vanish if our results are compared to the "low gravity'' CMFGEN models in the lower panel (see text).
Open with DEXTER

In summary, we find a very good agreement with the optical spectra from CMFGEN if the problem of different density stratifications is accounted for. The only disturbing fact is the strong difference in the He I singlets for dwarfs between 36 000 to 41 000 K and for supergiants between 31 000 to 35 000 K.

Although it is presently not clear which profiles are "correct'' or whether the truth lies in between, we like to point out the following. In our analyses of Galactic O-stars (Repolust et al. 2004), no problems were found in matching both the observed singlet and triplet lines in dwarfs. Concerning the supergiants, we had a problem for almost all stars cooler than O6, namely the well-known "generalized dilution effect'' (see the discussion and references in Repolust et al. 2004). Briefly, we could fit all He I lines (singlets and triplets) in parallel with the He II lines, except for He I 4471 (triplet) which was predicted to be too weak. One might argue that this is a symptom of generally incorrect He I lines, and speculate that this problem is related to the inconsistency seen here. Assuming that the He I-singlets produced by FASTWIND are erroneous it might then be possible to fit all He I singlets and the $\lambda$ 4471 triplet at cooler temperatures. In this case, however, we (and CMFGEN) would encounter the problem that the other triplet lines would be too strong and the He II lines too weak.

Presently, there is no way out of this dilemma other than to perform a number of detailed comparisons, with respect to both the models and the observations. Since the actual problem concerns the ratio of triplet to singlet lines and the problem is most pronounced for dwarfs, it should be possible to find a solution by comparing the theoretical predictions for this ratio (in terms of equivalent widths) as a function of $T_{\rm eff}$ vs. the observed ratio as a function of spectral type for a large sample of stars. Such work is in progress now.

Figure 20: Iso-contours of equivalent widths for He I 4471, as predicted by FASTWIND, using results from our model-grid with Helium abundance $Y_{\rm He}$ = 0.1 and solar metallicity for background elements. Left: negligible wind. Right: typical O-supergiant wind. Note the effect of wind-emission, shifting the iso-contours to the left. The locations of the corresponding models are indicated by asterisks.

Open with DEXTER

   
9 Model grids

As already outlined in Sect. 1, the parameter space to be investigated for the analysis of one object alone is large and almost prohibitive for the detailed analysis of very large samples of stars which have recently been collected (e.g., by means of the multi-object spectrograph FLAMES). Alternatively, a somewhat coarser analysis by means of the "traditional'' model-grid method is still applicable if an appropriate grid can be constructed. In this section, we will give some suggestions for this objective and report on first progress.

Although the presence of a wind introduces a large number of additional parameters to be considered in a fine fit ($\dot M$, $v_\infty$, $\beta$ and $R_\star$), there is a fortunate circumstance which allows for the construction of such model-grids with only one more parameter compared to grids from hydrostatic, plane-parallel models, at least if we do not aim at the analysis of specific (UV) resonance lines.

As has been shown by, e.g., Puls et al. (1996, see also de Koter et al. 1998; Schmutz et al. 1989, for diversifications), the wind-emission from recombination dominated transitions (so-called $\rho^2$-lines) remains rather unaffected by the specific choice of the individual values of $\dot M$, $v_\infty$ and $R_\star$ as long as the wind-strength parameter Q (also denoted as the "optical depth invariant''),

$$% Q=\frac{\dot M}{(v_\inftyR_\star)^{\frac{3}{2}}},$$ (93)

does not vary. In this case, most of the other lines also preserve their shape. An example is given in Fig. 18, where we have varied the mass-loss rate of a model of $\alpha$ Cam (cf. Repolust et al. 2004) by a factor of two (and accordingly the radius by 2^1.5) with almost no effect on the resulting H/He spectrum.

This behaviour (i.e., spectrum (and emergent fluxes) depend almost exclusively on Q and not on its individual constituents) follows from the fact that

• $\rho^2$-dependent line processes (e.g., recombination lines and resonance lines from ions one stage below the major one) scale with Q in the wind regime;
• the wind density scales with $\dot M/(v_\inftyR_\star^2)$ (continuity equation); and
• (resonance) lines from major ions scale with $\dot M/(v_\infty^2 R_\star)$, e.g. Hamann (1981).

Thus, the common power of "1.5'' with respect to $v_\infty$ and $R_\star$ used in Q, which refers to the scaling of $\rho^2$ lines, is also the best compromise to deal with the other physical parameters affecting a stellar model (most importantly, the line-blocking which depends both on density and line opacity).

Exploiting this knowledge, we have constructed a set of nine model-grids for the analysis of H/He profiles with three different helium abundances, $Y_{\rm He}$ = 0.1, 0.15 and 0.2, and three different background metallicities, z = 1.0, 0.5 and 0.2 (cf. Sect. 3), respectively. Each grid with given helium abundance and metallicity is three-dimensional with respect to the parameters  $T_{\rm eff}$, $\log g$ and $\log Q$, and the grid-spacing is roughly equidistant. The individual values for parameters incorporated into $\log Q$(which are actually needed to calculate a specific model) and additional ones have been assumed according to present knowledge:

• $R_\star$ from "empirical'' values, as a function of spectral type ( $T_{\rm eff}$) and luminosity class ($\log g$);
• $v_\infty$ as a function of photospheric escape velocity  $v_{\rm esc}$,
  

  $$% v_\infty= C(T_{\rm eff}) \cdot v_{\rm esc},$$ (94)

  
  
  in accordance with the results collected by Kudritzki & Puls (2000);
• velocity exponent $\beta$ from empirical values (see also Kudritzki & Puls 2000, and references therein), with $\beta = 0.9$ (as a compromise) for O-stars and increasing values towards later types;
• $\dot M$ from $\log Q$, $R_\star$ and $v_\infty$ as specified above;
• micro-turbulence $v_{\rm turb}$ = 15 km s^-1 throughout the grid as a compromise between O and B stars.

Our present grids comprise the range 20 000 K $\le$ $T_{\rm eff}$ $\le$ 50 000 K with $\Delta$  $T_{\rm eff}$ = 2500 K, $\log g$ between 2.2 $\le$ $\log g$ $\le$ 3.2 at $T_{\rm eff}$ = 20 000 K and 4.0 $\le$ $\log g$ $\le$ 4.5 at $T_{\rm eff}$ = 50 000 K.

The position of all models can be inferred from Fig. 20. With respect to $\log Q$ we have used values with -14.0 $\le$ $\log Q$ $\le$-11.4 ( $\Delta \log Q$ = 0.35 in most cases), where the lowest value corresponds to an almost negligible wind and the highest one to almost Wolf-Rayet conditions.

The denotation is such that we specify a letter for the wind density ("A'' to "H'', with densities $\log Q$ = -14.0, -13.5, -13.15, -12.8, -12.45, -12.1, -11.75, -11.4, respectively, if $\dot M$ is calculated in $M_\odot$/yr, $v_\infty$ in km s^-1 and $R_\star$ in $R_\odot$). Effective temperature and gravity are denoted by two numbers each. Thus, model "E2730'' refers to $\log Q$ = -12.45, $T_{\rm eff}$ = 27 500 K and $\log g$ = 3.0. Typical O-type supergiants correspond to series "E'', and typical B-type supergiants to series "D''.

For all these models we have calculated H/He profiles and equivalent widths in the optical and the IR. Thus, by simply over-plotting observed vs. simulated spectra one estimates the parameters  $T_{\rm eff}$, $\log g$, $Y_{\rm He}$ and wind-strength if the background metallicity is specified and the theoretical profiles have been convolved accounting for rotational broadening and resolution. In this way, the coarse analysis of one star is possible within few minutes and might be fine-tuned by calculating specific models (particularly with respect to $\beta$ if inferable from the emission line shapes).

In addition, a plot of various iso-contours of calculated equivalent widths gives deeper insight into certain dependencies. As an example, Fig. 20 shows the effect of wind emission on He I 4471. Further examples, particularly with respect to the spectral type classification criterium of O-star, $\log W' = \log ({\rm ew} 4471) - \log ({\rm ew} 4541)$, are given in Massey et al. (2005).

We intend to make these grids publicly available in the near future when the problem regarding the He I singlets has been solved.

   
10 Summary and outlook

In this paper we have described all updates applied to our previous version of FASTWIND (Paper I) regarding the approximative treatment of metal line-blocking/blanketing and the calculation of a consistent temperature structure.

The problem of line-blocking has been tackled in two steps. First, the occupation numbers of background elements are calculated by an approximative solution of the corresponding equations of statistical equilibrium with the option that the most abundant elements are treated almost "exactly'', i.e., by means of the Sobolev transport for line processes. Compared to alternative approaches (cf. Sect. 4) our method allows for the treatment of different spin systems, radially and frequency dependent radiation temperatures and a consistent ALI-iteration scheme. We have tested our solutions by comparing the approximative results with results from exact solutions and have not found any major discrepancies.

The occupation numbers derived in this way are subsequently used to calculate the line-blocked radiation field, again in an approximative way. To this end, we have formulated suitable means for the opacities (in analogy to Rosseland means but for frequency intervals not larger than 1000...1500 km s^-1) and emissivities (two-level-atom approach), and the resulting pseudo-continuum of overlapping lines is treated by means of a conventional continuum radiative transfer. Specific problems inherent in our approach (regarding a rigorous statistical description) have been pointed out and might lead to inaccurate solutions in a few cases. Investigations to improve our approach are presently under way in our group, as discussed in Sect. 5.

Our new version of FASTWIND allows for the calculation of a consistent temperature structure by applying a flux-correction method in the lower atmosphere and the electron thermal balance in the outer one. Regarding optical H/He lines, no major differences have been found compared to our previous NLTE Hopf-function method (cf. Paper I; and Repolust et al. 2004).

Due to the approximations applied and as intended, the performance of our code is very fast. The total computational time (starting all models from scratch) is of the order of 30 min on a PC with a 2 GHz processor if only H and He lines are considered as explicit ions, whereas the inclusion of other elements (e.g., Urbaneja 2004) into the "explicit'' treatment requires an additional 5 to 10 min each.

The new methods have been extensively tested by comparing with results from WM-Basic and CMFGEN, concerning temperature stratification, fluxes, number of ionizing photons and optical H/He profiles (comparison with CMFGEN only).

We have highlighted the importance of photospheric line-pressure, which is incorporated into the FASTWIND models and neglected in the standard version of CMFGEN, if not coupled to the plane-parallel code TLUSTY (see Sect. 1). Particularly, we have found indications that the use of the Sobolev approximation (within the force-multiplier concept) in WM-Basic can lead to an underestimate of this quantity, as already predicted by Pauldrach et al. (1986). On the other hand, the density/velocity stratification resulting from our approach (smoothly connecting the quasi-static photosphere and a $\beta$-law wind) agrees surprisingly well with the hydrodynamic structure as calculated from a consistent solution if $\beta$ is not too different from the "canonical'' value of 0.8... 1.0.

All three codes predict almost identical temperature structures and fluxes for $\lambda >$ 400 Å, whereas at lower wavelengths certain discrepancies are found. Compared to WM-Basic (using an identical line list for the background elements), our supergiant models differ only in the He II continua, where the FASTWIND-fluxes are somewhat larger, but still lower than the corresponding fluxes from CMFGEN. Since fluxes and corresponding numbers of ionizing photons can be extremely sensitive to subtle model differences in this wavelength regime, we consider any uncritical use of these quantities as unreliable.

Major discrepancies are also found in the range 300 Å < $\lambda$ < 400 Å, i.e., in the O II continuum and at the He II 304 resonance line. Compared to both WM-Basic and CMFGEN, our dwarf models produce less flux in this region (more blocking or less re-emission), whereas the supergiant models of FASTWIND and WM-Basic agree very well. The supergiant models of CMFGEN, on the other hand, show much less blocking which might point to some missing opacity. Again, the H I and He I continua agree very well in all three codes.

For the optical H/He lines, the coincidence between FASTWIND and CMFGEN is remarkable, except for the He I singlets in the temperature range between 36 000 to 41 000 K for dwarfs and between 31 000 to 35 000 K for supergiants, where CMFGEN predicts much weaker singlets. Up to now, the origin of this discrepancy could not be identified, but work is under way to solve this problem.

Although it is reassuring that the different codes agree well with respect to most of their predictions, this is only part of the story. One particularly disturbing fact concerns the present mismatch between the parameters obtained from an analysis in the optical and the UV, respectively. In the majority of cases, the UV gives lower effective temperatures, i.e., of the order of 2000 to 4000 K, if one compares the analyses of Galactic stars performed by Bianchi & Garcia (2002) and Garcia & Bianchi (2004) with results from Repolust et al. (2004) ( WM-Basic vs. FASTWIND) and the corresponding work for Magellanic Cloud stars by Hillier et al. (2003) and Bouret et al. (2003) (partly including also the optical range) with the results from Massey et al. (2004,2005) ( CMFGEN vs. FASTWIND). (Interestingly, the work by Crowther et al. 2002 ( CMFGEN) indicates higher temperatures for MC supergiants than derived by Massey et al. 2005.)

Part of this discrepancy (if combined UV/optical analyses are compared) might be related to the He I singlet vs. triplet problem as discussed above. Note, however, that this would account only for discrepancies in certain domains of the $T_{\rm eff}$ space and would typically result in maximum differences of the order of 2000 K, as has been found from a number of test calculations performed by one of us (J.P.) and F. Najarro (using CMFGEN), which will be reported on in a forthcoming publication. Moreover, the temperature scale for O-type dwarfs as derived by Martins et al. (2002) using CMFGEN and concentrating on the classification criterium He I 4471 (triplet) vs. He II 4541 is actually 1000 to 2000 K hotter than the calibration by Repolust et al. (2004).

In a recent paper, Martins et al. (2004) have discussed the uncertainties in $T_{\rm eff}$ obtained by relying on different diagnostic tools in the UV, analyzing four SMC-N81 dwarfs of spectral types O6.5 to O8.5. From the specific values derived from the UV-color index, the ionization balance of O IV/ V and Fe IV/ V and the N V1238/1242 and C III1426/1428 doublets, respectively, they quote a typical uncertainty of $\pm$3000 K in  $T_{\rm eff}$, which might easily account for part of the discrepancies with the optical.

Unfortunately, it is rather difficult to compare the differences obtained so far in a strict one-to-one case, simply because the corresponding samples barely overlap. In particular, a large fraction of the objects analyzed by means of CMFGEN are somewhat extreme, comprising either supergiants with (very) dense winds (Crowther et al. 2002) or dwarfs with very thin winds (Martins et al. 2004). The analysis of SMC stars by Bouret et al. (2003), on the other hand, covers only a sample of 6 dwarfs, in contrast to the larger sample by Massey et al. (2004,2005), and, therefore, it is not clear in how far selection effects do play a role. Finally, it is interesting to note that at least for one object in common, the O4I(f) star $\zeta$ Pup (HD 66811), the different analyses give almost identical results (Crowther et al. 2002; Repolust et al. 2004; and Pauldrach et al. in prep., analyzing the UV by means of WM-Basic).

Thus, we conclude that the present status of hot star parameters is not as clear as we would like it to be. We need to understand a number of additional physical processes and their influence on the derived parameters. Most important are the direct and indirect effects of the line-driven wind instability, i.e., the formation and interaction of clumps and shocks leading to X-ray emission and enhanced EUV-flux in the wind (e.g., Feldmeier et al. 1997; Pauldrach et al. 2001). Although incorporated to some extent into present codes, there are too many questions to be answered before we can consider these problems as solved. To give only two examples: We do not know the spatial distribution of the "clumping factor'', and also the X-ray emission is only on the verge of being understood (e.g., Oskinova et al. 2004; Kramer et al. 2003).

Before these effects can be treated in a realistic way, we need to primarily rely on diagnostic tools that are least "contaminated'', i.e., to concentrate on weak lines formed in the stellar photospheres (except, of course, the mass-loss indicators which will always be affected by clumping). Future investigations of O-type stars performed by FASTWIND will have to utilize not only H and He but also metal lines, as already incorporated in the analysis of B-stars (cf. Sect. 1). Particularly, one of the most important tools will be nitrogen with its strong sensitivity even at higher temperatures where He I begins to fail. Work in this direction is under way.

Acknowledgements

We like to thank a number of colleagues for their enduring willingness to discuss problems and provide assistance. Most important in this respect was and is Dr. Keith Butler, the living compendium in atomic physics. Particular thanks to Dr. Adi Pauldrach for providing his atomic data base (tailored for WM-Basic). Many thanks also to Drs. Paco Najarro and Tadziu Hoffmann for performing a number of test calculations with CMFGEN and WM-Basic, and their stimulating discussions concerning NLTE-effects. Finally, we would like to thank both our anonymous referee and Drs. Phil Massey, Alex de Koter and John Hillier for valuable comments on the manuscript.
J.P. appreciates support by NATO Collaborative Linkage Grant No. PST/CLG 980007. M.A.U. acknowledges financial support for this work by the Spanish MCyT under project PNAYA2001-0436, R.V. acknowledges support from the University of La Plata by a FOMEC grant (Pr. 724/98), and T.R. gratefully acknowledges financial support in form of a grant by the International Max-Planck Research School on Astrophysics (IMPRS), Garching.

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)
• Aufdenberg, J. P., Hauschildt, P. H., Baron, E., et al. 2002, ApJ, 570, 344 [NASA ADS] [CrossRef] (In the text)
• Baschek, B., v. Waldenfels, W., & Wehrse, R. 2001, A&A, 371, 1084 [EDP Sciences] [NASA ADS] [CrossRef] (In the text)
• Bianchi, L., & Garcia, M. 2002, ApJ, 581, 610 [NASA ADS] [CrossRef]
• Bouret, J.-C., Lanz, T., Hillier, D. J., et al. 2003, ApJ, 595, 1182 [NASA ADS] [CrossRef] (In the text)
• Butler K., & Giddings, J. R. 1985, Coll. Comp. Project No. 7 (CCP7), Newsletter 9, London, 7 (In the text)
• Castor, J., Abbott, D. C., & Klein, R. I. 1975, ApJ, 195, 157 [NASA ADS] [CrossRef]
• 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)
• de Koter, A., Heap, S. R., & Hubeny, I. 1998, ApJ, 509, 879 [NASA ADS] [CrossRef]
• Drew, J. 1985, MNRAS, 217, 867 [NASA ADS]
• Drew, J. 1985, ApJS, 71, 267 [NASA ADS] [CrossRef]
• Feldmeier, A., Puls, J., & Pauldrach, A. W. A. 1997, A&A, 322, 878 [NASA ADS]
• Fullerton, A. W., Crowther, P. A., De Marco, O., et al. 2000, ApJ, 538, L43 [NASA ADS] [CrossRef]
• Gabler, R., Gabler, A., Kudritzki, R. P., et al. 1989, A&A, 226, 162 [NASA ADS] (In the text)
• Garcia, M., & Bianchi, L. 2004, ApJ, 606, 497 [NASA ADS] [CrossRef]
• Gräfener, G., Koesterke, L., & Hamann, W.-R. 2002, A&A, 387, 244 [EDP Sciences] [NASA ADS] [CrossRef] (In the text)
• Grevesse, N., & Sauval, A. J. 1998, SSR, 85, 161 [NASA ADS]
• Hamann, W.-R. 1981, A&A, 100, 169 [NASA ADS] (In the text)
• Hauschildt, P. H., & Baron, E. 1999, J. Comp. Appl. Math., 109, 41 [CrossRef] (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. 2003, ASPC, 288, 199 [NASA ADS] (In the text)
• Hillier, D. J., & Miller, D. L. 1998, ApJ, 496, 407 [NASA ADS] [CrossRef] (In the text)
• Hillier, D. J., Lanz, T., Heap, S. R., et al. 2003, ApJ, 588, 1039 [NASA ADS] [CrossRef]
• Hubeny, I. 1992, in The Atmospheres of Early-Type Stars, ed. U. Heber, & C. S. Jeffery (Berlin: Springer), 377 (In the text)
• Hubeny, I., & Lanz, T. 1995, ApJ, 39, 875 [NASA ADS] (In the text)
• Hummer, D. G., & Rybicky, G. B. 1985, ApJ, 293, 258 [NASA ADS] [CrossRef]
• Jokuthy, A. 2002, Diploma Thesis, University Munich (In the text)
• Kramer, R. H., Cohen, D. H., & Owocki, S. P. 2003, ApJ, 592, 532 [NASA ADS] [CrossRef]
• Kubát, J., Puls, J., & Pauldrach, A. W. A. 1999, A&A, 341, 587 [NASA ADS] (In the text)
• Kudritzki, R. P., & Puls, J. 2000, ARA&A, 38, 613 [NASA ADS] (In the text)
• Lenorzer, A., Mokiem, M. R., de Koter, A., & Puls, J. 2004, A&A, 422, 275 [EDP Sciences] [NASA ADS] [CrossRef] (In the text)
• Lucy, L. B. 1999, A&A, 345, 211 [NASA ADS]
• Lucy, L. B., & Abbott, D. C. 1993, ApJ, 405, 738 [NASA ADS] [CrossRef]
• Martins, F., Schaerer, D., & Hillier, D. J. 2002, A&A, 382, 999 [EDP Sciences] [NASA ADS] [CrossRef] (In the text)
• Martins, F., Schaerer, D., Hillier, D. J., et al. 2004, A&A, 420, 1087 [EDP Sciences] [NASA ADS] [CrossRef]
• Massey, P., Kudritzki, R. P., Bresolin, F., et al. 2004, ApJ, 608, 1001 [NASA ADS] [CrossRef]
• Massey, P., Puls, J., Pauldrach, A. W. A., et al. 2005, ApJ, submitted
• Mazzali, P. A., & Lucy, L. B. 1993, A&A, 279, 447 [NASA ADS]
• Mazzali, P. A. 2000, A&A, 363, 705 [NASA ADS]
• Mihalas, D. 1975, Stellar atmospheres, 2nd edition (San Francisco: Freeman) (In the text)
• Mihalas, D., Kunasz, P. B., & Hummer, D. G. 1975, ApJ, 202, 465 [NASA ADS] [CrossRef] (In the text)
• Oskinova, L. M., Feldmeier, A., & Hamann, W.-R. 2004, A&A, 422, 675 [EDP Sciences] [NASA ADS] [CrossRef]
• Pauldrach, A. W. A., Puls, J., & Kudritzki, R. P. 1986, A&A, 164, 86 [NASA ADS]
• Pauldrach, A. W. A., Lennon, M., Hoffmann, T. L., et al. 1998, in Boulder-Munich Workshop II, ed. I. Howarth, PASP Conf. Ser., 131, 258
• Pauldrach, A. W. A., Hoffmann, T. L., & Lennon, M. 2001, A&A, 375, 161 [EDP Sciences] [NASA ADS] [CrossRef] (In the text)
• Przybilla, N., & Butler, K. 2004, ApJ, 609, 1181 [NASA ADS] [CrossRef]
• Puls, J., & Hummer, D. G. 1988, A&A, 191, 87 [NASA ADS]
• Puls, J., Kudritzki, R. P., Herrero, A., et al. 1996, A&A, 305, 171 [NASA ADS] (In the text)
• Puls, J., & Pauldrach, A. W. A. 1990, in Proc. Boulder Munich Workshop, ed. C. Garmany, PASPC, 7, 202 (In the text)
• Puls, J., Springmann, U., & Lennon, M. 2000, A&A, 141, 23 [NASA ADS] (In the text)
• Repolust, T., Puls, J., & Herrero, A. 2004, A&A, 415, 349 [EDP Sciences] [NASA ADS] [CrossRef]
• Rybicki, G. B., & Hummer, D. G. 1991, A&A, 245, 171 [NASA ADS] (In the text)
• Santolaya-Rey, A. E., Puls, J., & Herrero, A. 1997, A&A, 323, 488 [NASA ADS] (Paper I) (In the text)
• Schaerer, D., & Schmutz, W. 1994, A&A, 288, 231 [NASA ADS] (In the text)
• Schaerer, D., & de Koter, A. 1997, A&A, 322, 598 [NASA ADS]
• Schaerer, D., Meynet, G., Maeder, A., & Schaller, G. 1993, A&AS, 98, 523 [NASA ADS] (In the text)
• Schmutz, W. 1991, in Stellar Atmospheres: Beyond Classical Models, ed. L. Crivellari, I. Hubeny, & D. G. Hummer (Dordrecht: Kluwer), NATO ASI Ser. C, 341, 191 (In the text)
• Schmutz, W., Hamann, W.-R., & Wessolowski, U. 1989, A&A, 210, 236 [NASA ADS]
• Seaton, M. 1958, MNRAS, 118, 504 [NASA ADS] (In the text)
• Springmann, U. 1997, Thesis, University Munich (In the text)
• Taresch, G., Kudritzki, R. P., Hurwitz, M., et al. 1997, A&A, 321, 531 [NASA ADS] (In the text)
• Trundle, C., Lennon, D. J., Puls, J., et al. 2004, A&A, 417, 217 [EDP Sciences] [NASA ADS] [CrossRef]
• Urbaneja, M. A. 2004, Ph.D. Thesis, University of La Laguna
• Urbaneja, M. A., Herrero, A., Bresolin, F., et al. 2003, ApJ, 584, L73 [NASA ADS] [CrossRef]
• Wehrse, R., von Waldenfels, W., & Baschek, B. 1998, JQSRT, 60, 963 (In the text)