© ESO, 2012

1. Introduction

Since 1970, Kurucz with coworkers published extensive grids of model atmospheres assuming the local thermodynamic equilibrium (LTE), cf. Kurucz et al. (1974), Kurucz (1979, 1999) and Castelli & Kurucz (2003), and distributed computing codes (ATLAS5, ATLAS9, ATLAS12) for the common astronomical community, cf. Kurucz (1970, 1993) and websites of Kurucz¹ and Castelli². These fully line-blanketed LTE models of atmospheres are calculated using Kurucz’s atomic data for energy levels, transition probabilities and line broadening parameters for many atoms at several ionization stages and molecules. The data are progressively updated (cf. Kurucz 2005a, 2011) and more realistic opacities and equation of state are available now from Kurucz’s website. A huge number of spectral lines (cf. Kurucz 2011) can be included in calculations by using either opacity distribution functions (ODFs) or an opacity sampling (OS) approach. The first method is commonly used to construct models of atmospheres, whereas the second one is employed to obtain synthetic spectra for a given model of atmosphere. The new version of the ATLAS12 code (Kurucz 2005a) is dedicated to calculating atmospheric models within the OS approach. In this case opacities are calculated by the line-by-line method and no synthetic spectra must by recalculated because the radiation field at the final iteration should also give the correct emerging energy distribution, which can be compared with observations.

Non-LTE models of atmospheres and computing codes were elaborated and disseminated by Mihalas and Hubeny with coworkers (Mihalas 1972; Mihalas et al. 1975; Hubeny 1988; Hubeny & Lanz 1995, 2006; Lanz & Hubeny 2003, 2007). The Hubeny code is known as TLUSTY (Hubeny 1988) and the latest distributed version has the number 200, cf. TLUSTY website³. The state-of-the-art of line-blanketed non-LTE model atmospheres were recently calculated by Lanz & Hubeny (2003, 2007) for stars with T_eff ≥ 15   000 K. Although extensive model atoms of most abundant elements are used, non-LTE synthetic spectra were recalculated for a more complete list of line transitions with the best known broadening mechanisms of bound-bound transitions and were then distributed for general use through the TLUSTY website.

Since 1992, when the first results were published by Rogers & Iglesias (1992) and Iglesias & Rogers (1996) from the project based on the OPAL computing code, and by Seaton (1987), Seaton et al. (1994) and Badnell et al. (2005) from the international opacity project (OP), and by Hummer et al. (1993) from the international project of IRON devoted to compute ab-initio electron excitation cross-sections for astrophysical applications, new opportunities arise for modeling internal structure, evolution and pulsations of stars, as well as for the modeling of stellar atmospheres. Indeed, several cross-sections from the OP and IRON projects were implemented in the TLUSTY code, whereas Kurucz’s atomic data for Fe are included in OPAL and OP data. According to Seaton & Badnell (2004), the two sets of the Rosseland-mean opacities, OPAL and OP, agree between each other up to a few per cent at least up to lgT ≈ 6 despite differences in calculations of their atomic data and equation of states (EOS). In particular, the OPAL approach to the EOS is based on the many-body quantum statistical mechanics of partially ionized plasmas (cf. Rogers et al. 1996), whereas OP EOS involves the occupation probability formalism, which accounts for the destruction of loosely bound states of atoms by interactions with particles in the plasma (cf. Hummer & Mihalas 1988). The occupation probability formalism has already been used by Bergeron et al. (1991) and Tremblay & Bergeron (2009) to study hydrogen lines in atmospheres of white dwarf stars and has been implemented in the TLUSTY code by Hubeny & Lanz (1995), cf. also Hubeny et al. (1994).

In this paper we search for compatibility of LTE and non-LTE atmosphere models with stellar envelope models of early-type stars. The latter ones were calculated for the Rosseland-mean opacities and EOS taken from OPAL and OP projects. We also calculated LTE models of atmospheres using OP monochromatic opacities and OP EOS. The OP versus Kurucz treatment of autoionizing transitions has already been discussed by Castelli & Kurucz (2004). In the OP/IRON project autoionizing transitions are treated as resonances in photoionization cross-sections, cf. e.g. Bautista (1997). To make the problem tractable, each term is represented by only one level with its average energy. Kurucz treats autoionization as a line opacity with smooth continuum opacity, cf. Castelli & Kurucz (2004). In Kurucz’s approach, levels appropriate for autoionizing transitions with different quantum numbers J are included individually⁴.

The plan of the paper is as follows. In Sect. 2, we briefly compare published models of atmospheres with the internal structure of stars for the same effective temperature, T_eff and surface gravity, lgg. Significant differences are found in optically thin and optically thick layers. We therefore calculated LTE models of atmospheres with the OS approach using the ATLAS12 computer code implemented with the occupation probability formalism (Sect. 3.1), bound-free opacities of metals from the OP/IRON project (Sect. 3.2) and detailed profiles of line absorption coefficients of H I, He I and He II (Sect. 3.3). The new ATLAS/OP code is also developed for OP monochromatic opacities and OP EOS (Sect. 3.4). Examples of the new LTE models of atmospheres of B- and O-type stars are shown in Sect. 4.1 compared with OPAL and OP models of stellar envelopes. In Sect. 4.2, ATLAS/OP models are discussed for B stars. Remarks about LTE and non-LTE models of atmospheres calculated with the TLUSTY200 computer code are made in Sect. 4.3. Synthetic spectra of the analyzed models are also plotted in Sect. 4. Finally, conclusions are given in Sect. 5.

Fig. 1 Panel a) The Rosseland-mean opacities, κR    [cm2   g-1] , as a function of lgT for Castelli & Kurucz (2003) (marked as CK line) and Lanz & Hubeny (2007) (LH line) models of ${\mathit{T}}_{\mathrm{eff}}^{}\mathrm{=}\mathrm{25}\hspace{0.17em}\mathrm{000}$ K and lgg = 4.0 in comparison with OPAL (OPAL line) and OP (OP line) stellar envelope models. A non-LTE model calculated for a more complex model of O V ion (cf. Sect. 4.3) is also shown (T line). Panel b) The same as panel a), but on κR vs. lgR diagram. Panel c) The same as panel a), but for the radiation pressure gradients (in cgs units) as a function of lgT.

2. Comparison of published models of atmospheres with models of envelopes

In this section the published line-blanketed LTE models of atmospheres calculated by Castelli & Kurucz (2003) and non-LTE models of Lanz & Hubeny (2003, 2007) will be briefly reviewed and compared with models of stellar envelopes calculated for the same effective temperatures and surface gravities. The atmosphere models were taken from the Castelli and TLUSTY websites. The envelope models were calculated with computer codes (ENV/OPAL and ENV/OP) developed according to the Paczynski (1969) approach, but updated to OPAL and OP opacities and EOS. OPAL data were taken from the OPAL website⁵ and the OP data from the CDOP v. 3.3 base (Seaton 2005) described in Sect. 3.4. All OPAL models shown in this paper have the same normalized abundances of metals, but mass fractions of hydrogen and helium are specified by X and Y numbers, respectively, cf. Table 1, where details are given together with the “standard” chemical compositions ([m/H]  = 0) of the models of CK (Castelli & Kurucz 2003) and LH (Lanz & Hubeny 2007). Note that abundances of He, C, N, O, S, and Ar relative to H taken as the standard one by different authors are not the same. The abundance parameters of X, Y and Z are also shown in Table 1. As mentioned above, these numbers can be specified for OPAL data, while the chemical mixture for OP data can be fixed by users for all 17 most abundant elements (H, He, C, N, O, Ne, Na, Mg, Al, Si, S, Ar, Ca, Cr, Mn, Fe, and Ni) and therefore are not shown in Table 1. The reference chemical composition assumed in this paper is shown in the last column of Table 1.

As an example, Fig. 1 shows CK and LH models of atmospheres of T_eff = 25   000 K and lgg = 4.0 in comparison with stellar envelope models calculated for OPAL and OP data with X = 0.7347 and Z = 0.0170. Panels a–c on this figure show the diagrams κ_R vs. lgT, κ_R vs. lgR and lgg_rad vs. lgT, respectively. The well-known Rosseland-mean opacity, κ_R, is defined as $\begin{array}{ccc}{\mathit{\kappa }}_{\mathrm{R}}^{-1}& \mathrm{=}& \frac{\mathit{\pi }}{\mathrm{4}{\mathit{a}}_{\mathrm{S}\mathrm{-}\mathrm{B}}{\mathit{T}}^{\mathrm{3}}}\mathrm{\int }\mathrm{(}{\mathit{\kappa }}_{\mathit{\nu }}^{\mathrm{a}}\mathrm{+}{\mathit{\kappa }}_{\mathit{\nu }}^{\mathrm{s}}{\mathrm{)}}^{-1}\left(\frac{\mathrm{d}{\mathit{B}}_{\mathit{\nu }}}{\mathrm{d}\mathit{T}}\right)\mathrm{d}\mathit{\nu ,}\\ & \mathrm{=}& \mathrm{\int }{\left({\mathit{\kappa }}_{\mathit{\nu }}^{\mathrm{a}}\mathrm{+}{\mathit{\kappa }}_{\mathit{\nu }}^{\mathrm{s}}\right)}^{-1}\mathit{F}\mathrm{\left(}\mathit{u}\mathrm{\right)}\mathrm{d}\mathit{u,}\end{array}$where $\mathit{F}\mathrm{\left(}\mathit{u}\mathrm{\right)}\mathrm{=}\frac{\mathrm{15}}{\mathrm{4}{\mathit{\pi }}^{\mathrm{4}}}\frac{{\mathit{u}}^{\mathrm{4}}\exp \mathrm{(}\mathrm{-}\mathit{u}\mathrm{)}}{{\mathrm{[}\mathrm{1}\mathrm{-}\exp \mathrm{(}\mathrm{-}\mathit{u}\mathrm{\left)}\mathrm{\right]}}^{\mathrm{2}}}\mathit{,}$(3)and $\mathit{u}\mathrm{=}\frac{\mathit{h\nu }}{\mathit{kT}}\mathrm{·}$(4)Here, ${\mathit{\kappa }}_{\mathit{\nu }}^{\mathrm{a}}$ and ${\mathit{\kappa }}_{\mathit{\nu }}^{\mathrm{s}}$ are monochromatic absorption and scattering coefficients, respectively. B_ν(T) denotes the Planck function at temperature T and frequency ν and a_S−B is the Stefan-Boltzmann constant. The quantity R is defined by the local values of density ρ and T: $\lg \mathit{R}\mathrm{=}\lg \mathit{\rho }\mathrm{-}\mathrm{3.0}\hspace{0.17em}\lg \mathit{T}\mathrm{+}\mathrm{18.0.}$(5)As usual, the diffusion approximation is assumed to calculate the gradient of radiation pressure, g_rad, in stellar envelope codes $\frac{\mathrm{d}{\mathit{P}}_{\mathrm{rad}}}{\mathrm{d}\mathit{m}}\mathrm{=}\frac{\mathrm{d}}{\mathrm{d}\mathit{m}}\left(\frac{\mathrm{4}}{\mathrm{3}\mathit{c}}{\mathit{a}}_{\mathrm{S}\mathrm{-}\mathrm{B}}{\mathit{T}}^{\mathrm{4}}\right)\mathrm{=}\frac{{\mathit{a}}_{\mathrm{S}\mathrm{-}\mathrm{B}}}{\mathit{c}}{\mathit{\kappa }}_{\mathrm{R}}{\mathit{T}}_{\mathrm{eff}}^{\mathrm{4}}\mathrm{=}{\mathit{g}}_{\mathrm{rad}}\mathit{,}$(6)whereas $\frac{\mathrm{d}{\mathit{P}}_{\mathrm{rad}}}{\mathrm{d}\mathit{m}}\mathrm{=}\frac{\mathrm{4}\mathit{\pi }}{\mathit{c}}\mathrm{\int }\left({\mathit{\kappa }}_{\mathit{\nu }}^{\mathrm{a}}\mathrm{+}{\mathit{\kappa }}_{\mathit{\nu }}^{\mathrm{s}}\right){\mathit{H}}_{\mathit{\nu }}\mathrm{d}\mathit{\nu }\mathrm{=}\frac{{\mathit{a}}_{\mathrm{S}\mathrm{-}\mathrm{B}}}{\mathit{c}}{\mathit{\kappa }}_{\mathit{H}}{\mathit{T}}_{\mathrm{eff}}^{\mathrm{4}}\mathrm{=}{\mathit{g}}_{\mathrm{rad}}$(7)for modeling of stellar atmospheres, cf. Mihalas (1978). Here, m means the mass depth, H_ν the monochromatic Eddington flux of radiation, κ_H the flux-mean opacity and T_eff the effective temperature of a star.

There are two main sources of discrepancy between the analyzed models, viz. the temperature distribution at optically thin regions and opacity effects at great optical depths. We first consider the temperature as a function of depth. The grey model of atmosphere with the optical scale corresponding to the Rosseland-mean opacities (diffusion approximation) is assumed for stellar envelope calculations. In the plane-parallel limit we have ${\mathit{T}}^{\mathrm{4}}\mathrm{\left(}{\mathit{\tau }}_{\mathrm{R}}\mathrm{\right)}\mathrm{=}\frac{\mathrm{3}}{\mathrm{4}}{\mathit{T}}_{\mathrm{eff}}^{\mathrm{4}}\mathrm{(}{\mathit{\tau }}_{\mathrm{R}}\mathrm{+}\mathit{q}\mathrm{(}{\mathit{\tau }}_{\mathrm{R}}\mathrm{\left)}\mathrm{\right)}\mathit{,}$(8)where Hopf function q(τ_R) is in the range from 0.577 at the surface to 0.710 at great optical depths, cf. Mihalas (1978). This temperature distribution should be compared with the exact relation for the plane-parallel LTE models in radiative equilibrium, which can be expressed as ${\mathit{T}}^{\mathrm{4}}\mathrm{\left(}{\mathit{\tau }}_{\mathit{H}}\mathrm{\right)}\mathrm{=}\frac{\mathrm{3}}{\mathrm{4}}{\mathit{T}}_{\mathrm{eff}}^{\mathrm{4}}\frac{{\mathit{\kappa }}_{\mathit{J}}}{{\mathit{\kappa }}_{\mathit{B}}}\left(\frac{{\mathit{\tau }}_{\mathit{H}}}{\mathrm{3}{\mathit{f}}_{\mathit{K}}}\mathrm{+}\frac{\mathrm{1}}{\mathrm{3}{\mathit{f}}_{\mathit{H}}\mathrm{\left(}\mathrm{0}\mathrm{\right)}}\right)\mathit{,}$(9)where τ_H is the optical depth corresponding to the flux-mean opacity, κ_J/κ_B is the ratio of the J- to B-mean opacities; f_K = K/J is the Eddington factor for J (zeroth moment, mean intensity integrated over frequency) and K (second moment, so called K-integral), and f_H(0) = H(0)/J(0). The CK models were calculated up to τ_R = 100, where T(τ_R = 100) ≈ 3T_eff for all models of early-type stars. At great optical depths Eqs. (8) and (9) lead to the same results due to the following relations κ_H ≈ κ_R, τ_H ≈ τ_R, J ≈ B, ${\mathit{f}}_{\mathit{K}}\mathrm{\approx }\frac{\mathrm{1}}{\mathrm{3}}$ and one can assume ${\mathit{f}}_{\mathit{H}}\mathrm{\left(}\mathrm{0}\mathrm{\right)}\mathrm{=}\frac{\mathrm{1}}{\mathrm{2}}$ (Eddington approximation). The same is true for the gradients of radiation pressure, cf. Eqs. (6) and (7). One can therefore expect that atmospheric and stellar envelope models of early-type stars should agree very well with each other at the bottom of the photosphere if the Rosseland-mean opacities and densities are the same. Figure 1 shows that this is not the case. Note, however, the good agreement between OPAL and OP envelope models.

At optically thin regions the diffusion approximation (Eq. (8)) breaks down and temperature distributions of our stellar envelope models are not correct. Equation (9) indicates that the true temperature distribution of LTE models in radiative equilibrium depends on the radiation field in a much more complicated manner than the diffusion approximation predicts. Again, the same is true for the gradient of radiation pressure. Furthermore, at optically thin regions non-LTE effects may also play an important role. However, the LH model (cf. Fig. 1) does not agree either with the envelope nor CK models at optically thick regions, where LTE conditions and diffusion approximations are fulfilled with high accuracy. Clearly, the non-LTE models of Lanz & Hubeny (2003, 2007) underestimate stellar opacities.

As mentioned above, Fig. 1b shows κ_R versus lgR defined by Eq. (5). The stellar envelope models are plotted up to lgT = 5.68, while the CK and LH models ended at lgT = 4.87 and 5.0, respectively. The two local maxima of stellar opacities at lgT = 4.6 and lgT = 5.3 (cf. Fig. 1a) result in two loops seen on the diagram plotted in Fig. 1b. The quantity of lgR can also be expressed as lgR = lgP_g − 4lgT + lgμ + 33.86, where P_g is the gas pressure and μ – mean atomic mass. The gradient of the total pressure, $\frac{\mathrm{d}\mathrm{(}{\mathit{P}}_{\mathrm{g}}\mathrm{+}{\mathit{P}}_{\mathrm{rad}}\mathrm{)}}{\mathrm{d}\mathit{m}}$, follows from the hydrostatic equilibrium equation. At optically thick regions the gradient of radiation pressure mimics the behavior of κ_R (cf. Eq. (6)) and therefore affects gas pressure, density, and stellar opacity in layers where κ_R has local maxima. This is the reason why two loops are seen in Fig. 1b.

The gradient of the radiation pressure plotted in Fig. 1c can also be compared with the gravity acceleration of lgg = 4.0, which means that the models shown in this figure are well inside the hydrostatic equilibrium. The role of the gradient of the radiation pressure increases for stars with lower lgg and higher T_eff values.

We also calculated OPAL and OP Rosseland-mean opacities for the input parameters of lgT and lgρ taken from CK models. This approach has the advantage that differences in gradients of radiation pressure, geometry of models, etc. do not influence the results. Discrepancies between the CK and OPAL, and OP data approach a factor of about 2 at lgT ≈ 5.06 for the model of T_eff = 45   000 K and lgg = 5.

More details about OPAL and OP data are shown in Fig. 2a−c, where the Rosseland-mean opacities, the ratio of OP to OPAL Rosseland-mean opacities and the ratio of gas pressures as functions of lgT are plotted. OP data were calculated for the input parameters (lgT and lgρ) taken from ENV/OPAL model of envelope shown in Fig. 2a. The abundance parameters equal to X = 0.7080 and Z = 0.0172 and OPAL mixture of metals listed in Table 1 were used in both cases. Evidently, OPAL and OP opacities differ from each other by a few percent at atmospheric layers, but differences can reach about 20 per cent at lgT ≈ 5.5. Although the OPAL and OP equations of state are based on different approaches, the resulting gas pressures are almost the same, cf. Fig. 2c. The agreement is better than 0.1 per cent for the chemical composition used here.

Fig. 2 Panel a) The Rosseland-mean opacities, κR    [cm2   g-1] , as a function of lgT for OPAL (thick line) and OP (thin line) models of Teff = 25   000 K and lgg = 4.0. The chemical abundance parameters X = 0.708 and Z = 0.0172 are assumed. Panel b) The ratio of OP to OPAL Rosseland-mean opacities are displayed for input parameters of (lgT and lgρ) taken from the OPAL model shown in panel a). Panel c) The same as panel b), but for gas pressures.

3. Modifications of the ATLAS12 code

3.1. ATLAS12/MHD-Q code

Most of the basic physical assumptions remain the same as those used in the ATLAS12 computing code distributed by Kurucz (cf. Sect. 1). Only significant changes from this code will be described here. Assuming LTE, we need to know the atomic internal partition functions for a given element k at different ionization states j${\mathit{U}}_{\mathit{j,k}}\mathrm{=}{\mathrm{\Sigma }}_{\mathit{i}} {\mathit{W}}_{\mathit{i,j,k}}{\mathit{g}}_{\mathit{i,j,k}}{\mathrm{e}}^{\mathrm{-}{\mathit{E}}_{\mathit{i,j,k}}\mathit{/}\mathrm{k}\mathit{T}}$(10)to calculate occupations of atomic levels (Boltzman equations) ${\mathit{N}}_{\mathit{i,j,k}}\mathrm{=}{\mathit{W}}_{\mathit{i,j,k}}{\mathit{g}}_{\mathit{i,j,k}}{\mathrm{e}}^{\mathrm{-}{\mathit{E}}_{\mathit{i,j,k}}\mathit{/}\mathrm{k}\mathit{T}}{\mathit{N}}_{\mathit{j,k}}\mathit{/}{\mathit{U}}_{\mathit{j,k}}$(11)and ionization stages (Saha equations) $\frac{{\mathit{N}}_{\mathit{j,k}}\mathit{/}{\mathit{U}}_{\mathit{j,k}}}{{\mathit{N}}_{\mathit{j}\mathrm{+}\mathrm{1}\mathit{,k}}\mathit{/}{\mathit{U}}_{\mathit{j}\mathrm{+}\mathrm{1}\mathit{,k}}}\mathrm{=}{\left(\frac{{\mathit{h}}^{\mathrm{2}}}{\mathrm{2}\mathit{\pi }\mathit{mk}\mathit{T}}\right)}^{\mathrm{3}\mathit{/}\mathrm{2}}{\mathrm{e}}^{\mathrm{-}{\mathit{E}}_{\mathit{j,k}}\mathit{/}\mathit{kT}}{\mathit{N}}_{\mathrm{e}}\mathit{/}\mathrm{2}$(12)of species, cf. Mihalas (1978). Formally, the partition functions diverges if the sum extends over an infinite number of bound states of an isolated atom. In reality, the highest levels cannot remain bound because they are strongly perturbed by neighboring particles. Several theories were developed, cf. Hummer & Mihalas (1988). In particular, Kurucz (1970) assumed that levels with very high quantum number have a definite cutoff (W_i,j,k = 0) and that bound levels (W_i,j,k = 1) maintain their original values. This effect was described in terms of lowering of the ionization energy by the amount of $\mathrm{\Delta }\mathit{E}\mathrm{=}\mathrm{1.44}\mathrm{\times }{\mathrm{10}}^{-7} {\mathrm{Z}}_{\mathrm{a}}\mathit{/}\mathit{D,}$(13)where ΔE is in eV, Z_a is 1 for neutral particles, 2 for those of once ionized, etc., and D is the Debye length, $\mathit{D}\mathrm{=}{\left(\frac{\mathit{k}\mathrm{T}}{\mathrm{4}\mathit{\pi }{\mathit{e}}^{\mathrm{2}}}\right)}^{\mathrm{1}\mathit{/}\mathrm{2}}\mathrm{[}{\mathit{N}}_{\mathrm{e}}\mathrm{+}{\mathrm{\Sigma }}_{\mathit{i}} {\mathit{Z}}_{\mathit{i}}^{\mathrm{2}}{\mathit{N}}_{\mathit{i}}{\mathrm{]}}^{\mathrm{-}\mathrm{1}\mathit{/}\mathrm{2}}\mathit{,}$(14)cf. also Mihalas (1978). At low temperatures the upper levels are not populated, therefore the value of the cutoff may have a small effect on the sum in Eq. (10), but at high temperatures, the cutoff determines the value of the partition function. A more detailed discussion of this point is given in Sect. 4.1.

In the ATLAS codes partition functions of all species from H (Z = 1) up to Es (Z = 99) for neutral atoms and a few first ionization stages are calculated. The data are continuously updated by Kurucz and the latest versions of the ATLAS codes include partition functions of sufficiently high ionization stages to obtain the ion populations of elements up to T ≈ 200   000 K, cf. Kurucz (2005b, 2011). In this paper models calculated with the code based on ATLAS12 with the Kurucz partition functions are denoted as ATLAS12/KPF.

As mentioned above, we implemented the occupation probability formalism (cf. Hummer & Mihalas 1988; Mihalas et al. 1988, 1990; Däppen et al. 1988; Nayfonov et al. 1999) into the ATLAS12 code. This theory is known as the MHD (Mihalas, Hummer & Däppen) formalism. Hubeny et al. (1994) extended this approach to non-LTE plasma. The partition functions are again given by Eq. (10), but now W_i,j,k means the probability that level i remains bounded despite of perturbations by neighboring particles. Levels with low values of W_i,j,k are said to be dissolved, viz. (1 − W_i,j,k) is the probability that the level  i is dissolved. Hummer & Mihalas (1988) argued that Stark splitting is the most efficient mechanism leading to dissolution of levels. The occupation probability of level i is then expressed in terms of the electric microfield distribution function W(β) as ${\mathit{W}}_{\mathit{i,j,k}}^{\mathrm{s}}\mathrm{=}{\mathrm{\int }}_{\mathrm{0}}^{{\mathit{\beta }}_{\mathit{i,j,k}}}\mathit{W}\mathrm{\left(}\mathit{\beta }\mathrm{\right)}\mathrm{d}\mathit{\beta }\mathrm{=}\mathit{Q}\mathrm{\left(}{\mathit{\beta }}_{\mathit{i,j,k}}\mathrm{\right)}\mathit{.}$(15)The Q(β) = Q_H(β) functions are tabulated e.g. by Hummer (1986) for the Holtsmark theory of the distribution function of W(β) = W_H(β). Here, β measures the electric microfield at the distance r from the core with a charge of Z_a of radiated ions in units of the mean field strength for perturber ions with a charge Z_p and density N_p. For the purpose of the estimations given below we note that for low β${\mathit{W}}_{\mathrm{H}}\mathrm{\left(}\mathit{\beta }\mathrm{\right)}\mathrm{=}\frac{\mathrm{3}\mathit{\pi }}{\mathrm{3}}{\mathrm{\Sigma }}_{\mathit{l}}\mathrm{(}\mathrm{-}\mathrm{1}{\mathrm{)}}^{\mathit{l}}\mathrm{\Gamma }\left(\frac{\mathrm{4}\mathit{l}\mathrm{+}\mathrm{6}}{\mathrm{3}}\right)\frac{{\mathit{\beta }}^{\mathrm{2}\mathit{l}\mathrm{+}\mathrm{2}}}{\mathrm{(}\mathrm{2}\mathit{l}\mathrm{+}\mathrm{1}\mathrm{)}\mathrm{!}}\mathit{,}$(16)whereas for high β${\mathit{W}}_{\mathrm{H}}\mathrm{\left(}\mathit{\beta }\mathrm{\right)}\mathrm{=}\mathrm{1.496}{\mathit{\beta }}^{\mathrm{-}\mathrm{5}\mathit{/}\mathrm{2}}\mathrm{(}\mathrm{1}\mathrm{+}\mathrm{5.107}{\mathit{\beta }}^{\mathrm{-}\mathrm{3}\mathit{/}\mathrm{2}}\mathrm{+}\mathrm{14.43}{\mathit{\beta }}^{-3}\mathrm{+}{\mathit{...}}^{\mathrm{)}}\mathit{,}$(17)cf. Mihalas (1978). These series expansions of W_H(β) lead to ${\mathit{Q}}_{\mathrm{H}}\mathrm{\left(}{\mathit{\beta }}_{\mathit{i,j,k}}\mathrm{\right)}\mathrm{\approx }\mathrm{0.1415}{\mathit{\beta }}_{\mathit{i,j,k}}^{\mathrm{3}}\mathrm{-}\mathrm{0.0392}{\mathit{\beta }}_{\mathit{i,j,k}}^{\mathrm{5}}\mathrm{+}\mathrm{0.0077}{\mathit{\beta }}_{\mathit{i,j,k}}^{\mathrm{7}}\mathrm{+}\mathit{...}$(18)and ${\mathit{Q}}_{\mathrm{H}}\mathrm{\left(}{\mathit{\beta }}_{\mathit{i,j,k}}\mathrm{\right)}\mathrm{\approx }\mathrm{1.0}\mathrm{-}\mathrm{0.9973}{\mathit{\beta }}_{\mathit{i,j,k}}^{\mathrm{-}\mathrm{3}\mathit{/}\mathrm{2}}\mathrm{-}\mathrm{2.5467}{\mathit{\beta }}_{\mathit{i,j,k}}^{-3}\mathrm{-}\mathrm{4.7972}{\mathit{\beta }}_{\mathit{i,j,k}}^{\mathrm{-}\mathrm{9}\mathit{/}\mathrm{2}}\mathrm{+}\mathit{...}$(19)for β_i,j,k ≪ 1 (relevant to highly excited levels) and β_i,j,k ≫ 1 (appropriate for low-lying levels).

When plasma correlation effects are important, the post-Holtsmark microfield distribution is used to upgrade the original formalism, cf. Hubeny et al. (1994) and Nayfonov et al. (1999). The electric microfield distribution, W(β;Z_r,a), depends now on two additional parameters, the radiator charge Z_r and the correlation parameter $\mathit{a}\mathrm{=}\mathrm{0.09}{\mathit{N}}_{\mathrm{e}}^{\mathrm{1}\mathit{/}\mathrm{6}}{\mathit{T}}^{\mathrm{-}\mathrm{1}\mathit{/}\mathrm{2}}\mathit{.}$(20)Results of numerical integrations of Eq. (15) with W(β) = W(β;Z_r,a) can be approximated by ${\mathit{W}}_{\mathit{i,j,k}}^{\mathrm{pH}}\mathrm{=}\frac{\mathit{f}}{\mathrm{1}\mathrm{+}\mathit{f}}\mathit{,}$(21)where $\begin{array}{ccc}\mathit{f}& & \mathrm{=}\frac{\mathrm{0.1402}\mathrm{(}\mathit{x}\mathrm{+}\mathrm{4}{\mathit{Z}}_{\mathrm{r}}{\mathit{a}}^{\mathrm{3}}\mathrm{)}{\mathit{\beta }}_{\mathit{i,j,k}}^{\mathrm{3}}}{\mathrm{1}\mathrm{+}\mathrm{0.1285}\mathit{x}{\mathit{\beta }}_{\mathit{i,j,k}}^{\mathrm{3}\mathit{/}\mathrm{2}}}\mathit{,}\\ \mathit{x}& & \mathrm{=}\mathrm{(}\mathrm{1}\mathrm{+}\mathit{a}{\mathrm{)}}^{\mathrm{3.15}}\end{array}$and ${\mathit{Z}}_{\mathrm{r}}\mathrm{=}{\mathit{Z}}_{\mathit{a}}\mathrm{-}\mathrm{1.}$(24)The readers are kindly referred to the papers of Hubeny et al. (1994) and Nayfonov et al. (1999), where more details are given. For a Holtsmark distribution (a = Z_r = 0) this formula fits the series expansions given by Eqs. (18) and (19) to 1 per cent or better for β_i,j,k < 1 and β_i,j,k > 3, respectively.

For normal stellar compositions of elements one can neglect the effect of perturber ions other than protons, setting Z_p = 1 and N_ion = N_e. Another choice discussed by Hummer & Mihalas (1988) is to distribute the positive charge equally over all ions, so that each would have a charge of $\begin{array}{}\mathrm{˜}\\ {\mathit{Z}}_{\mathrm{p}}\end{array}\mathrm{=}{\mathrm{\Sigma }}_{\mathit{Z}}\mathit{Z}{\mathit{N}}_{\mathit{Z}}\mathit{/}{\mathrm{\Sigma }}_{\mathit{Z}}{\mathit{N}}_{\mathit{Z}}\mathrm{=}{\mathit{N}}_{\mathrm{e}}\mathit{/}{\mathit{N}}_{\mathrm{ion}}\mathit{.}$(25)Under the second circumstances the critical value of β can be estimated as ${\mathit{\beta }}_{\mathit{i,j,k}}\mathrm{=}\frac{{\left(\frac{{\mathit{Z}}_{\mathit{a}}^{\mathrm{2}}}{\mathrm{2}{\mathit{n}}^{\mathrm{2}}}\right)}^{\mathrm{2}}}{\mathrm{4}{\mathit{Z}}_{\mathit{a}}}{\left(\frac{\mathrm{4}\mathit{\pi }{\mathit{a}}_{\mathrm{0}}^{\mathrm{3}}}{\mathrm{3}}\right)}^{\mathrm{-}\mathrm{2}\mathit{/}\mathrm{3}}{\mathit{N}}_{\mathrm{e}}^{-1}{\mathit{N}}_{\mathrm{ion}}^{\mathrm{1}\mathit{/}\mathrm{3}} {\mathit{K}}_{\mathit{n}}\mathit{,}$(26)where N_e = Σ_ZZN_Z, and N_ion = Σ_ZN_Z are the total electron and ion densities, respectively, a₀ is the Bohr radius, n means the principal quantum number of the state of i, and the quantal Stark-ionization correction factor is given by ${\mathit{K}}_{\mathit{n}}\mathrm{=}\mathrm{1} \mathrm{for} \mathit{n}\mathrm{\le }\mathrm{3}$(27)and ${\mathit{K}}_{\mathit{n}}\mathrm{=}\frac{\mathrm{16}}{\mathrm{3}}\frac{\mathit{n}}{\mathrm{(}\mathit{n}\mathrm{+}\mathrm{1}{\mathrm{)}}^{\mathrm{2}}} \mathrm{for} \mathit{n}\mathit{>}\mathrm{3.}$(28)We now estimate the value of the principal quantum number n′ of a level for which the occupation probability decreases to 0.9. Setting N_ion = 0.9N_e into Eq. (26) and using Eq. (19) one can find β_n′ = 5.3 and ${\mathit{n}}^{\mathrm{\prime }}\mathrm{=}\mathrm{964} {\mathit{N}}_{\mathrm{e}}^{\mathrm{-}\mathrm{2}\mathit{/}\mathrm{15}}{\mathit{Z}}_{\mathit{a}}^{\mathrm{3}\mathit{/}\mathrm{5}}\mathit{.}$(29)For example, at the bottom of an atmosphere model of T_eff = 25   000 K and lgg = 4.0, where N_e ≈ 10¹⁶   cm^-3, one can obtain n′ ≈ 7 for Z_a = 1, whereas n′ = 61 at optically thin regions, where N_e ≈ 10⁹   cm^-3, respectively. On the other hand, one can expect that the largest atomic orbit cannot extend beyond the nearest ion, which leads to $\mathrm{⟨}\mathit{r}{\mathrm{⟩}}_{\mathrm{nl}}\mathrm{\approx }\frac{{\mathit{a}}_{\mathrm{0}}}{\mathrm{2}{\mathit{Z}}_{\mathit{a}}}\mathrm{[}\mathrm{3}{\mathit{n}}^{\mathrm{2}}\mathrm{-}\mathit{l}\mathrm{(}\mathit{l}\mathrm{+}\mathrm{1}\mathrm{\left)}\mathrm{\right]}\mathit{<}{\mathit{r}}_{\mathrm{0}}\mathrm{\approx }\mathrm{(}\frac{\mathrm{3}}{\mathrm{4}\mathit{\pi }{\mathit{N}}_{\mathrm{ion}}}{\mathrm{)}}^{\mathrm{1}\mathit{/}\mathrm{3}}$ and ${\mathit{n}}_{\max }\mathrm{\approx }\mathrm{9000} {\mathit{N}}_{\mathrm{e}}^{\mathrm{-}\mathrm{1}\mathit{/}\mathrm{6}}{\mathit{Z}}_{\mathit{a}}^{\mathrm{1}\mathit{/}\mathrm{2}}\mathit{.}$(30)Here, the mean radius of an atom  ⟨ r ⟩ _nl was estimated for the orbital quantum number l = 0. The occupation probability quickly decreases with n for a given value of the electron density N_e. From Eq. (18) one can obtain ${\mathit{W}}_{\mathit{i,j,k}}^{\mathrm{s}}\mathrm{\approx }\mathrm{1.22}\mathrm{\times }{\mathrm{10}}^{\mathrm{46}} {\mathit{N}}_{\mathrm{e}}^{-2} {{\mathit{Z}}_{\mathit{a}}}^{\mathrm{9}} {\mathit{n}}^{-15}\mathit{,}$(31)as far as β_i,j,k < 1. This indicates that the occupation probability decreases by several orders of magnitude for levels between n′ to n_max, viz. ${\mathit{W}}_{{\mathit{n}}_{\max }}^{\mathrm{s}}\mathit{/}{\mathit{W}}_{{\mathit{n}}^{\mathrm{\prime }}}^{\mathrm{s}}\mathrm{=}\mathrm{(}{\mathit{n}}^{\mathrm{\prime }}\mathit{/}{\mathit{n}}_{\max }{\mathrm{)}}^{\mathrm{15}}$ and effective truncation of the partition functions takes place due to decrease of the occupation probability of levels before the largest atomic orbit with n_max has the tendency to extend beyond the nearest ion. For ${\mathit{W}}_{\mathit{i,j,k}}^{\mathrm{s}}\mathrm{=}{\mathrm{10}}^{-4}$, one can find $\mathit{n}″\mathrm{\approx }\mathrm{2.2}\mathrm{\times }{\mathrm{10}}^{\mathrm{3}} {\mathit{N}}_{\mathrm{e}}^{\mathrm{-}\mathrm{2}\mathit{/}\mathrm{15}} {\mathit{Z}}_{\mathit{a}}^{\mathrm{3}\mathit{/}\mathrm{5}}\mathit{,}$(32)which is near to Kurucz’s empirical formula for cutoff energy levels in opacity calculations.

In addition to the static perturbations of radiated atoms (ions) described above, inelastic collisions of radiated atoms (ions) with electrons should be considered for energy levels with high values of n. In the Hummer & Mihalas (1988) model of such collisions the probability that the level is bound can be expressed as ${\mathit{W}}_{\mathit{n}}^{\mathrm{c}}\mathrm{=}\exp \left[\mathrm{-}\mathrm{4.84}\mathrm{\times }{\mathrm{10}}^{-17} {\mathit{N}}_{\mathrm{e}}{\mathit{n}}^{\mathrm{6}}{{\mathit{Z}}_{\mathit{a}}}^{-4}{\mathit{f}}_{\mathit{n}}\mathrm{\left(}{\mathit{Z}}_{\mathit{a}}\mathit{,T}\mathrm{\right)}/\left(\mathrm{1}\mathrm{-}{\left(\frac{\mathit{n}}{{\mathit{n}}^{\mathrm{c}}\mathrm{+}\mathrm{1}}\right)}^{\mathrm{2}}\right)\right]\mathit{.}$(33)We interpolate f_n(Z_a,T) and n^c functions tabulated by Hummer & Mihalas (1988) to the input parameters N_e, T and Z_a assuming that the values of f_n(Z_a = 8,T) and n^c(Z_a = 8) can be used whenever Z_a ≥ 8. As pointed out by Hummer & Mihalas (1988), the model is not well defined for very high values of n and we simply put $\frac{\mathit{n}}{{\mathit{n}}^{\mathrm{c}}\mathrm{+}\mathrm{1}}\mathrm{=}\frac{{\mathit{n}}^{\mathrm{c}}}{{\mathit{n}}^{\mathrm{c}}\mathrm{+}\mathrm{1}}$ in Eq. (33) if n ≥ n^c.

For cool and high-density plasma collisions between neutral atoms may be important. Again, following Hummer & Mihalas (1988), the occupation probability can be approximated as ${\mathit{W}}_{\mathit{i,j}\mathrm{=}\mathrm{0}\mathit{,k}}^{\mathrm{n}}\mathrm{=}\exp \left[\mathrm{-}\frac{\mathrm{4}\mathit{\pi }}{\mathrm{3}}{\mathrm{\Sigma }}_{\mathit{j}\mathrm{=}\mathrm{0}\mathit{,}{\mathit{k}}^{\mathrm{\prime }}}{\mathit{N}}_{\mathit{j}\mathrm{=}\mathrm{0}\mathit{,}{\mathit{k}}^{\mathrm{\prime }}}\mathrm{(}{\mathit{r}}_{\mathit{i,j}\mathrm{=}\mathrm{0}\mathit{,k}}\mathrm{+}{\mathit{r}}_{\mathrm{1}\mathit{,j}\mathrm{=}\mathrm{0}\mathit{,}{\mathit{k}}^{\mathrm{\prime }}}{\mathrm{)}}^{\mathrm{3}}\right]\mathit{.}$(34)Here, j = 0 means neutral atoms, r_i,j = 0,k the radius of the radiated atom k at state i, while r_{1,j = 0,k′} corresponds to the radii of perturbed atoms k′ at the ground state with densities N_j = 0,k′.

Assuming that the static and dynamic perturbations are essentially independent, we use the product of ${\mathit{W}}_{\mathit{i,j,k}}\mathrm{=}{\mathit{W}}_{\mathit{i,j,k}}^{\mathrm{s}}{\mathit{W}}_{\mathit{n}}^{\mathrm{c}}{\mathit{W}}_{\mathit{i,j}\mathrm{=}\mathrm{0}\mathit{,k}}^{\mathrm{n}}$ to calculate the probability that the level is occupied. The last factor is equal to 1.0 for ionized species (j ≠ 0). In stellar atmospheres with the solar abundance of elements discussed in this paper, ${\mathit{W}}_{\mathit{n}}^{\mathrm{c}}$ and ${\mathit{W}}_{\mathit{i,j}\mathrm{=}\mathrm{0}\mathit{,k}}^{\mathrm{n}}$ are close to unity for low and moderate values of n, and it is sufficient to use only the function of ${\mathit{W}}_{\mathit{i,j,k}}^{\mathrm{s}}$. For lower densities, however, both n′ and n_max are high and the model of inelastic collisions could be considered. As supposed by Hummer & Mihalas (1988), it can influence the internal partition functions of elements. Collisions between neutral atoms may be important for cool plasma with a high density of neutral species.

We used the occupation probability formalism described above to calculate partition functions from neutral atoms to hydrogenic ions of the following elements: H–Si, S, Ar, and Ca. The necessary atomic data were taken from the OP TOPbase website⁶, where models of atoms with electron configurations up to n = 10 are available. The quantum defect description of states was applied for levels with higher values of n, up to n_max. In the OP TOPbase project extensive atomic data are also available for the ions of Fe III–Fe XXV, while one can find Fe II, Ti XX, Cr XXI, and Ni XXVI in the OP TIPbase project. Moreover, some ions of other elements (Sc–Ni) one can also find in the OP TIPbase website⁷. The atomic data at OP TIPbase are not as rich as for species from the OP TOPbase. Therefore we used the partition functions computed by Kurucz in 1988 that were recently updated to the ATLAS codes for the iron group elements: Sc I–IX, Ti I–IX, V I–IX, Cr I–IX, Mn I–IX, Fe I–IX, Co I–IX, Ni I–IX, and Cu I–IX. As mentioned above, these data include sufficiently high ionization stages to obtain ion populations up to T ≈ 200   000 K.

The occupation probability formalism also affects the spectroscopic properties of the plasma. According to Hubeny et al. (1994), the absorption coefficient at a frequency ν of an atom in state i is given by ${\mathit{\chi }}_{\mathit{i}}\mathrm{\left(}\mathit{\nu }\mathrm{\right)}\mathrm{=}{\mathit{N}}_{\mathit{i}}\mathrm{\left[}{\mathrm{\Sigma }}_{{\mathit{i}}^{\mathrm{\prime }}\mathit{>}\mathit{i}}{\mathit{W}}_{{\mathit{i}}^{\mathrm{\prime }}\mathit{,j,k}}{\mathit{\sigma }}_{\mathit{i}{\mathit{i}}^{\mathrm{\prime }}}\mathrm{\right(}\mathit{\nu }\mathrm{)}\mathrm{+}{\mathit{\sigma }}_{\mathit{ic}}^{\mathrm{tot}}\mathrm{(}\mathit{\nu }{\mathrm{\left)}}^{\mathrm{\right]}}\mathrm{(}\mathrm{1}\mathrm{-}{\mathrm{e}}^{\mathrm{-}\mathit{h\nu }\mathit{/}\mathit{kT}}\mathrm{)}\mathit{.}$(35)To a good approximation ${\mathit{\sigma }}_{\mathit{i}{\mathit{i}}^{\mathrm{\prime }}}\mathrm{=}\frac{\mathit{\pi }{\mathit{e}}^{\mathrm{2}}}{\mathit{mc}}{\mathit{f}}_{\mathit{i}{\mathit{i}}^{\mathrm{\prime }}}{\mathit{\phi }}_{\mathit{i}{\mathit{i}}^{\mathrm{\prime }}}\mathrm{\left(}\mathit{\nu }\mathrm{\right)}$(36)for bound-bound transitions, where f_ii′ is the oscillator strength and φ_ii′(ν) the normalized profile of the absorption coefficient. For bound-free transitions one can use ${\mathit{\sigma }}_{\mathit{ic}}^{\mathrm{tot}}\mathrm{\left(}\mathit{\nu }\mathrm{\right)}\mathrm{=}{\mathit{\sigma }}_{\mathit{ic}}\mathrm{\left(}\mathit{\nu }\mathrm{\right)} \mathrm{if} \mathit{\nu }\mathrm{\ge }{\mathit{\nu }}_{\mathit{ic}}$(37)and ${\mathit{\sigma }}_{\mathit{ic}}^{\mathrm{tot}}\mathrm{\left(}\mathit{\nu }\mathrm{\right)}\mathrm{=}\mathrm{[}\mathrm{1}\mathrm{-}{\mathit{W}}^{\mathrm{\ast }}\mathrm{(}{\mathit{n}}^{\mathrm{\ast }}\mathrm{\left)}\mathrm{\right]}{\mathit{\sigma }}_{\mathit{ic}}^{\mathrm{\ast }}\mathrm{\left(}\mathit{\nu }\mathrm{\right)} \mathrm{if} \mathit{\nu }\mathit{<}{\mathit{\nu }}_{\mathit{ic}}\mathit{.}$(38)Here, ν_ic is the ionization frequency for the unperturbed atom and W ^∗ (n ^∗ ) corresponding to the effective quantum number ${\mathit{n}}^{\mathrm{\ast }}\mathrm{=}{\left(\frac{\mathrm{1}}{{\mathit{n}}^{\mathrm{2}}}\mathrm{-}\frac{\mathit{h\nu }}{{\mathit{E}}_{\mathit{j,k}}}\right)}^{\mathrm{-}\mathrm{1}\mathit{/}\mathrm{2}}$(39)designates the highest level that can be reached from the state of i with the principal quantum number n by absorption of a photon of energy hν. We calculate W ^∗ (n ^∗ ) by interpolating the occupation probability W_i,j,k of levels i to the effective quantum number n ^∗ . ${\mathit{\sigma }}_{\mathit{ic}}^{\mathrm{\ast }}\mathrm{\left(}\mathit{\nu }\mathrm{\right)}$ is the extrapolated photoionization cross-sections σ_ic(ν) to the lower frequency side of the ionization limit ν_ic.

The absorption coefficient χ_i(ν) given by Eq. (35) is corrected for the stimulated emission and can be used under LTE conditions. The first term in Eq. (35) corresponds to the bound-bound transitions at frequency ν from the level of i to the undissolved (bound) part of the state i′, which occurs with the probability W_i′,j,k. The second term represents the bound-free transitions. Equation (36) refers to the “true” transitions between undissolved parts of the levels of i and i′, while the transitions between the undissolved to the dissolved components of the levels i and i′ are included in Eq. (35) with ${\mathit{\sigma }}_{\mathit{ic}}^{\mathrm{tot}}\mathrm{\left(}\mathit{\nu }\mathrm{\right)}$ given by Eq. (38). The transitions from the dissolved part of i to higher levels and the true free-free transitions are already included in the free-free opacities. We note that the formulae (35)–(38) differ from those used by Däppen et al. (1987) and Bergeron et al. (1991). How this influences the opacities for line transitions between low and high excited levels is discussed in Sect. 4.

There are several slightly different MHD-type occupation probability functions used in the literature. In addition to the post-Holtsmark distribution function used by Hubeny et al. (1994), Nayfonov et al. (1999), and Tremblay & Bergeron (2009), the true Holtsmark distribution, which gives ${\mathit{W}}_{\mathit{i,j,k}}^{\mathrm{s}}\mathrm{=}{\mathit{W}}_{\mathit{n}}^{\mathrm{H}}\mathrm{=}{\mathit{Q}}_{\mathrm{H}}\mathrm{\left(}{\mathit{\beta }}_{\mathit{i,j,k}}\mathrm{\right)}$, was assumed for studying white dwarfs by Bergeron et al. (1991). A simplified approximate form of the Holtsmark microfield for randomly distributed ions was introduced by Hummer & Mihalas (1988) to reduce the computational effort for the free energy of the gas. The resulting occupation probability function, denoted here as ${\mathit{W}}_{\mathit{n}}^{\mathrm{HM}}$, corresponds to the nearest-neighbor theory (cf. Unsöld 1948) but includes an additional empirical factor, which results in ${\mathit{W}}_{\mathit{n}}^{\mathrm{HM}}\mathrm{\approx }{\mathit{W}}_{\mathit{n}}^{\mathrm{H}}$, whenever ${\mathit{W}}_{\mathit{n}}^{\mathrm{H}}\mathrm{\ge }\mathrm{0.1}$. This approximate function was used by Seaton et al. (1994) to calculate OP EOS and OP opacities.

There is, however, a problem of consistency of line profiles and the occupation probability function ${\mathit{W}}_{\mathit{i,j,k}}^{\mathrm{s}}$ used for hydrogen as first pointed out by Seaton (1990), cf. also Bergeron (1993). The commonly accepted theory describing these line profiles is the unified theory of Stark broadening elaborated by Vidal et al. (1970, 1971, 1973), cf. also Lemke (1997), where a grid of hydrogen profiles are available. This theory is valid for the entire line profile from the impact regime in the core to the quasi-static far wings. To improve the consistency of the line profiles with the true Holtsmark or post-Holtsmark distributions used to calculate the occupation probability functions, one can use the function of ${\mathit{W}}_{\mathit{i,j,k}}^{\mathrm{HL}}\mathrm{=}{\mathit{W}}_{\mathit{i,j,k}}^{\mathrm{pH}}\mathrm{(}\mathrm{2} {\mathit{\beta }}_{\mathit{i,j,k}}\mathrm{;}{\mathit{Z}}_{\mathit{r}}\mathit{,a}\mathrm{)}$ with the additional factor of 2 introduced to the argument of the post-Holtsmark distribution, cf. Bergeron (1993) and Hubeny et al. (1994). This approximation applies only to hydrogen, assuming that both protons and electrons contribute to the profile equally. Another approach, originally proposed by Seaton (1990), uses ${\mathit{W}}_{\mathit{i,j,k}}^{\mathrm{s}}\mathrm{=}{\mathit{W}}_{\mathit{i,j,k}}^{\mathrm{HM}}\mathrm{\approx }{\mathit{W}}_{\mathit{i,j,k}}^{\mathrm{H}}$, as in Hummer & Mihalas (1988) or ${\mathit{W}}_{\mathit{i,j,k}}^{\mathrm{s}}\mathrm{=}{\mathit{W}}_{\mathit{i,j,k}}^{\mathrm{pH}}$ as in Nayfonov et al. (1999), with line profiles renormalized in such a way that only electric microfields β with an amplitude inferior to the critical field β_i′,j,k for the upper level i′ of the transition contribute to the broadening of the lines by the protons. Another renormalization of line profiles was proposed by Tremblay & Bergeron (2009) for atmospheres of white dwarf stars.

Fig. 3 Panel a) The occupation probability function Wn (in logarithmic units) for the bound states of H I with the principal quantum number n for a typical physical parameters at τR ≈ 0.2 of a model of Teff = 25   000 K and lgg = 4.0. Explanation of the functions displayed is given in the text. Panel b) The same as panel a), but for He II.

Figure 3 shows MHD-type probability functions of ${\mathit{W}}_{\mathit{n}}^{\mathrm{H}}$ (marked as H lines), ${\mathit{W}}_{\mathit{n}}^{\mathrm{pH}}$ (pH lines), ${\mathit{W}}_{\mathit{i,j,k}}^{\mathrm{HL}}$ (HL line), ${\mathit{W}}_{\mathit{n}}^{\mathrm{HM}}$ (HM lines), together with ${\mathit{W}}_{\mathit{n}}^{\mathrm{c}}$ (Coll lines), respectively, for a model of T_eff = 25   000 K and lgg = 4.0 at τ_R ≈ 0.2, where lgN_e ≈ 14.4. Clearly, the two Holtsmark-type distribution functions, ${\mathit{W}}_{\mathit{n}}^{\mathrm{pH}}$ and ${\mathit{W}}_{\mathit{n}}^{\mathrm{H}}$, are almost the same, but their simplified approximate form, ${\mathit{W}}_{\mathit{n}}^{\mathrm{HM}}$, decreases with n much more sharply. Similarly, the two types of collisions described above result in sharp cutoff level populations, but for higher values of n. In all cases considered above, an effective truncation of the partition functions takes place due to decrease of the occupation probability of levels before the largest atomic orbit with n_max has the tendency to extend beyond the nearest ion, and before the above-mentioned theory of lowering of the ionization energy (Eqs. (13) and (14)) relaxes levels as unbounded ones. In Fig. 3 we also plot the Kurucz cutoff energy levels.

3.2. Bound-free transitions for metals

In ATLAS12/MHD-Q, bound-free transitions for metals (Z > 2) can be calculated either as in the original ATLAS12 code distributed by Kurucz or for large-scale data given by the OP/IRON project. The latter, available from the TOPbase at CDS (Cunto et al. 1993) or the XSTAR (Bautista & Kallman 2001) atomic databases, include detailed photoionization cross-sections with extensive resonance structures for all energy levels of bound states with n ≤ 10 and l ≤ 7. We used resonance-averaged photoionization cross-sections (cf. Bautista et al. 1998) for more than 1000 bound states of atoms and ions of C through Si (Z = 6−14), S (Z = 16) and Fe (Z = 26).

3.3. Bound-bound transitions

All bound-bound transitions for metals were calculated assuming the Voigt function with the natural, Stark, van der Waals, and thermal broadening of the line absorption coefficients. The broadening parameters are taken from Kurucz’s atomic base.

Detailed cross-sections tabulated by Lemke (1997) or by Shoening & Butler (after Hubeny & Lanz 2000) are implemented in the ATLAS12/MHD-Q code for H I bound-bound transitions. These profiles include the convolution of the Stark and thermal Doppler profiles. The approximate Stark broadening given by Peterson (after Kurucz 1970) or by Hubeny et al. (1994) as well as the procedure of HPROF4 developed by Peterson for the Stark effect of hydrogen lines included in the original ATLAS12 code can also be used.

Similarly, detailed profiles of selected He I lines given by Barnard et al. (1974, 1975) and Griem (1972), as well as profiles for the Balmer, Pashen and Pickering series of He II calculated by Schoening & Butler (1989a,b) are implemented. Helium lines can be also calculated under the assumption of the Voigt function or profiles with the approximate Stark broadening given by Hubeny et al. (1994).

3.4. ATLAS/OP code

ATLAS/OP is designed for computing LTE model atmospheres using monochromatic opacity and EOS data taken from the OP project distributed as the OPCD v. 3.3 base, cf. Seaton (2005). This OP base includes monochromatic cross-sections of the 17 most abundant elements (cf. Sect. 2) calculated for the temperature range of lgT = 3.5−8.0 with the step of 0.05 dex. The second parameter, lgN_e, is in the range specified for each value of lgT with the step of 0.5 dex. The OPCD base also contains several computer codes for calculating the Rosseland (and Planck) mean opacities as well as mass-densities for chemical compositions specified by the user. The accurate calculations of the Rosseland-mean opacities are made for equally-spaced points in the variable v, viz.: $\frac{\mathrm{1}}{{\mathit{\sigma }}_{\mathrm{R}}}\mathrm{=}{\mathrm{\int }}_{\mathrm{0}}^{\mathrm{\infty }}\frac{\mathrm{1}}{\mathit{\sigma }}\mathrm{d}\mathit{v,}$(40)where $\mathit{v}\mathrm{\left(}\mathit{u}\mathrm{\right)}\mathrm{=}{\mathrm{\int }}_{\mathrm{0}}^{\mathit{u}}\frac{\mathit{F}\mathrm{\left(}\mathit{u}\mathrm{\right)}}{\mathrm{[}\mathrm{1}\mathrm{-}\exp \mathrm{(}\mathrm{-}\mathit{u}\mathrm{\left)}\mathrm{\right]}}\mathrm{d}\mathit{u}$(41)and σ_R is the Rosseland-mean cross-section in atomic units (${\mathit{a}}_{\mathrm{0}}^{\mathrm{2}}$). The function of F(u) is given by Eq. (3) and the total cross-section, σ, for a mixture of elements with fractional abundances f_k is $\mathit{\sigma }\mathrm{=}{\mathrm{\Sigma }}_{\mathit{k}}{\mathit{\sigma }}_{\mathit{k}}{\mathit{f}}_{\mathit{k}}\mathit{,}$(42)where cross-sections σ_k are given in atomic units and include the true absorption (not corrected for stimulated emission) and scattering cross-sections (divided by the stimulated emission factor). The Rosseland-mean opacity per unit mass (κ_R) is equal to ${\mathit{\kappa }}_{\mathrm{R}}\mathrm{=}\frac{{\mathit{\sigma }}_{\mathrm{R}}}{\mathit{\mu }}\mathit{,}$(43)where μ is the mean atomic weight. Finally, the resulting lgκ_R and its derivatives are interpolated to the input values of lgT and lgρ using bi-cubic interpolations.

With the exception of H and He, monochromatic opacities σ_k are stored in the OPCD v. 3.3 base at 10 000 variably spaced points of the variable u from 0.073157 to 20. (For H and He another u-mesh points are used.) Of course, data for different values of lgT do not correspond to the same wavelength points. This does not matter for calculations of the Rosseland (and Planck) mean opacities for the fixed set of (lgT, lgN_e)-values and then interpolation to the input parameters of temperature and density. However, we need the radiation field at wavelength points simultaneously covering cool (optically thin) and hot (optically thick) regions, where T is changing by a factor of about 4. In the ATLAS/OP code we therefore added ATLAS opacities for the continuum at wavelength points outside the range of OPCD v. 3.3 data.

In the ATLAS/OP code we made bi-cubic interpolations of the total OP cross-sections for a given mixture of elements (corrected for scattering and stimulated emission) to the input parameters of lgT and lgρ. The interpolations of monochromatic opacities was performed in exactly the same manner as the original OP codes do for the Rosseland- and Planck-mean opacities.

As mentioned above, OP monochromatic data must also be interpolated to fixed wavelength points appropriate for the assumed T_eff. This can be performed either by interpolating the total opacities σ or individual opacities of each element σ_k (cf. Eq. (42)) to fixed wavelength points. The results are not the same due to the intermediate spectral resolution of the OP data used. Interpolations of the individual monochromatic opacities of each element, σ_k, to the needed wavelength points and then summing the total cross-section σ gives more accurate results.

Fig. 4 The internal partition functions (in logarithmic units) vs. lgT for ATLAS12/MHD-Q model of Teff = 15   000 K and lgg = 4.0 (thick lines) in comparison with Kurucz’s data shown as thin lines. The panels show H I, He I and II, and the first three ionization stages of C, N, O, Ne and Si, respectively. The lowest right-hand panel shows the same, but for Fe III–V. The ionization states of I–V are plotted by means of lines as explained in the right-hand side of the figure.

Fig. 5 The ionization fractions (in logarithmic units) of species for the same models as in Fig. 4. Symbols used are the same as in Fig. 4.

Fig. 6 Panel a) The Rosseland-mean opacities, κR    [cm2   g-1] , as a function of lgT for ATLAS12/MHD-Q (dotted line) in comparison with OP (continuous line) and OPAL (dashed line) envelope models of Teff = 25   000 K and lgg = 4.0. Panel c) The same as panel a), but OP/IRON bound-free absorptions of metals are included in the ATLAS12/MHD-Q model (dotted line). Panel e) The same as panelc), but OP opacities are used for hydrogen and helium. Panel g) The same as panel a), but the ATLAS12/OP model is plotted by means of an A/OP line. The right-hand side panels show the same as the left-hand side ones, but for the diagrams of κR vs. lgR.

Fig. 7 Same as Fig. 6, but for ${\mathit{T}}_{\mathrm{eff}}^{}\mathrm{=}\mathrm{45}\hspace{0.17em}\mathrm{000}$ K and lgg = 5.0. The CK model (plotted as A/CK line) is also shown in panel g).

Fig. 8 Panel a) Emerging flux of radiation, H  [erg   cm-2   s-1   nm-1] , near the Lyman series of H I for the ATLAS12/MHD-Q model of Teff = 25   000 K and lgg = 4.0. Panel b) The ratio of the emerging fluxes for models calculated with different description of the occupation probability formalism (see text). Panel c) The same as panel a), but for the Balmer series of H I. Panel d) The same as panel b), but for the Balmer series of H I.

Fig. 9 Emerging flux of radiation, H  [erg   cm-2   s-1   nm-1] , near the Lyman series of H I (panel a)) of the model of ATLAS12/MHD-Q shown in Figs. 6a and b is compared with models calculated with the ATLAS/PFK code with different profiles for H I, He I and He II line absorption coefficients (see text).

Fig. 10 Same as Fig. 9, but now for the emerging flux of radiation near the Balmer series of H I.

Fig. 11 OP monochromatic cross-sections (in atomic units) of H, He, C, N, O and Fe elements for lgT = 4.35 K and lgNe = 15.0 corrected for abundance parameters of the elements, σkfk. Panel a) shows the cross-sections for H (thick line) and He (dashed line). Panel b) shows data for C, N and O elements in comparison with the dominant absorber of H. Panel c) shows the same as panel b), but for Fe.

Fig. 12 Emerging flux of radiation, H  [erg   cm-2   s-1   nm-1] , near the Lyman series of H I for ATLAS/OP models shown in panels a) − d) for (Teff,lgg) = (20 000, 4.5), (28 000, 4.5), (250 000, 4.0), (25 000, 5.0), respectively.

Fig. 13 Same as Fig. 12, but near the Balmer limit of H I.

4. Results

4.1. LTE models of stellar atmospheres with the occupation probability formalism

Figure 4 shows examples of the internal partition functions, U_j,k, calculated in terms of the occupation probability formalism taking into account both the static and dynamic perturbations of species in comparison with the Kurucz data. The static perturbations were calculated for the post-Holsmark microfield of randomly distributed ions, while the dynamic perturbations caused by inelastic collisions of radiated atoms (ions) with electrons and collisions between neutral atoms were calculated as described in Sect. 3.1. Clearly, large differences are present for H I (U_0,1). In this case, however, hydrogen is markedly ionized and only a weak effect is seen on EOS. For other species, significant differences exist at optically thick layers. Note, however, the very good agreement of our results with partition functions calculated by Kurucz (cf. Sect. 3.1) for Fe III–V, as Fig. 4 illustrates. The same is true for Fe VI–IX. We found a similar behavior of the partition functions for all models of B stars. For H I, discrepancies decrease with lower values of T_eff and become small for T_eff ≈ 10   000 K, where the exponent factor in Eq. (10) is able to effectively cut off occupations of high levels, cf. Sect. 3.1. The partition functions shown in Fig. 4 have indeed only weak effect on the ionization fractions of He and metals for B stars, as Fig. 5 shows for the same model of T_eff = 15   000 and lgg = 4.0. The truncation effect in the partition functions remains significant only for the ionization fraction, N_j,k/Σ_jN_j,k, of H I and cancels exactly when the ratio of N_j,k/U_j,k (cf. Eq. (11)) is computed. For He and metals, differences arise for higher values of temperature and density. Quantitatively, this also depends on the MHD-type occupation probability function used. The ionization fractions generated by the ATLAS12/MHD-Q code with the ${\mathit{W}}_{\mathit{i,j,k}}^{\mathrm{s}}\mathrm{=}{\mathit{W}}_{\mathit{i,j,k}}^{\mathrm{HM}}$ function (cf. Sect. 3.1) agree well with Mihalas et al. (1988, 1990) results for very broad physical conditions. This means that both the Saha-type equation of state (cf. Kurucz 1970; Mihalas 1978) implemented by the occupation probability functions and the MHD-type equation of state, which is based on the minimization of the free-energy of gas, give the same results for much broader physical parameters than studied in this paper. More details will be given in a future paper, where properties of stellar envelope models will be discussed.

Figures 6 and 7 show the diagrams κ_R vs. lgT and κ_R vs. lgR for models calculated with the ATLAS12/MHD-Q, ENV/OPAL, and ENV/OP codes for T_eff = 25   000 and 45   000 K, respectively. Chemical mixtures corresponding to X = 0.708 and Z = 0.0172 were used, cf. Table 1. The ATLAS12/MHD-Q models shown in these figures were calculated with an opacity sampling (OS) approach implemented by MHD-Q formalism assuming ${\mathit{W}}_{\mathit{i,j,k}}^{\mathrm{s}}\mathrm{=}{\mathit{W}}_{\mathit{i,j,k}}^{\mathrm{HL}}$ for the occupation probability of levels. Atomic data of LOWLINES.DAT, HIGHLINES.DAT and GFALLNLTE.DAT taken from the Kurucz website were used for metals. In panels a and b, H I, He I and He II profiles tabulated by Schoening & Butler (after Hubeny & Lanz 2000), Griem (1972) and Schoening & Butler (1989a,b), when available, were used, otherwise approximate profiles given by Hubeny et al (1994) were assumed. Photoionization cross-sections from the original ATLAS12 code were used. Panels c and d show the same, but OP/IRON photoionization cross-sections (cf. Sect. 3.2) were used. ATLAS12/MHD-Q models calculated with total OP cross-sections of H and He and OP/IRON photoionization cross-sections for metals are shown in panels e and f. Finally, models calculated with the ATLAS12/OP code (cf. Sect. 3.4) are plotted in panels g and h. In Fig. 7g we also plot the Castelli & Kurucz (2003) model. Evidently, this model predicts a high opacity excess near lgT ≈ 5.06, which reflects differences in the atomic bases for bound-bound transitions for metals and possibly the effect of splitting the LS terms into individual total angular momentum J-levels for autoionizing transitions, cf. Sect. 1, instead of treating each LS term as a whole as in the OP cross-sections.

As mentioned above, the ATLAS12/MHD-Q emerging fluxes of radiation presented in this paper were taken from the last iteration during construction of the models. (In the last few iterations, high spectral resolution was used, taking into account 1.2 million wavelength points.) Figures 8–10 show examples of emerging flux of radiation near the limits of the H I Lyman and H I Balmer series. The behavior of the theoretical spectra at these regions is dominated by hydrogen line opacities and, in principle, it offers the possibility of studying the occupation probability formalism discussed in Sect. 3.1 in astrophysical conditions. In particular, Fig. 8a shows the emerging flux of radiation near the Lyman series of H I for the model of T_eff = 25   000 K and lgg = 4.0 assuming that the hydrogen absorption coefficients are given by Eqs. (35)–(38) with ${\mathit{W}}_{\mathit{i,j,k}}^{\mathrm{s}}\mathrm{=}{\mathit{W}}_{\mathit{i,j,k}}^{\mathrm{HL}}$ for occupation probability functions. Figure 8b illustrates how models are affected by the function of ${\mathit{W}}_{\mathit{i,j,k}}^{\mathrm{s}}\mathrm{=}{\mathit{W}}_{\mathit{i,j,k}}^{\mathrm{pH}}$ and Däppen’s et al. (1987) approach to line absorption coefficients mentioned in Sect. 3.1. Panels c and d of Fig. 8 show the same, but near the Balmer limit. This clearly influences the emerging spectra only near the Lyman and Balmer limits.

Now, we illustrate how the approximate profiles of H I, He I and He II affect atmospheric models of B stars. Figures 9 and 10 show the comparison of emerging fluxes of radiation for models constructed using the ATLAS12/MHD-Q and ATLAS12/KPF (cf. Sect. 3.1) codes. The former flux denoted as H is shown in panels a; the others are denoted as H(PFKx). We extended the partition functions of the ATLAS12/KPF code to all ions of species up to Z = 28 to obtain more consistent models at great optical depths. This happens only for highest ionization states of species, which are not present in the original code. We assumed that these missing partition functions can be approximated by the values of statistical weight of the ground states of the considered ions. The flux ratio of H(KPFf)/H is shown in panels b for the case when the same hydrogen and helium cross-sections, taken from ATLAS12/MHD-Q, were used in both codes. Panels c correspond to H(KPFd) with the Stark broadening of H I lines following the Kurucz (1970) and Voigt function for helium lines. Panels d shows the same as panel c, but the emerging flux H(PFKe) corresponds to the model calculated with the HPROF4 procedure for H I (cf. Sect. 3.3). As one can see, the results based on ATLAS12/KPF code with approximate cross-sections for H I, He I and He II differ markedly from those of corresponding to ATLAS12/MHD-Q model. It significantly affects the emerging flux distribution and the energy balance of the atmosphere. The former effect is of considerable importance if the models are to be employed in a subsequent analysis of stars. The latter effect changes the structure of the atmospheric model itself and hence the line profiles of H I, He I and He II, as well as metals, especially in UV range. It can to some degree influence the derived abundances and turbulent velocities from observed spectra of stars.

4.2. Examples of LTE models with OP monochromatic opacities and OP EOS for B stars

Figure 11 shows examples of OP monochromatic cross-sections, σ_k, corrected for abundance parameters f_k of the elements taken from the last column of Table 1 (cf. Eq. (42)) for lgT = 4.35 K and lgN_e = 15.0. The panels a−c display He, C–O and Fe in comparison with the dominant absorber of H (thick line), respectively. We note that some narrow lines can be missing for a given value of T due to moderate resolution in u. Other lines can appear when T is changing by 0.05 dex. as in the OPCD v 3.3 data. Better sampling over T and u is necessary to improve the quality of the data. In practical applications, however, this mainly influences synthetic spectra, but structures of atmospheric models are not influenced very much because very high spectral resolution is not necessary to calculate the total flux of radiation integrated over frequency with a good accuracy. Indeed, Figs. 6 and 7 show that the ATLAS/OP models fit the internal structure of stars calculated with the ENV/OP code well. This does not mean, however, that emerging spectra are able to fit observed high-resolution spectra of stars with a good accuracy. In addition to moderate sampling over u, T and N_e of OP data used, there is also a problem connected with the fact that wavelengths of bound-bound transitions are known only approximately because of the model atoms used in the OP/IRON project. An example is well illustrated in Fig. 12, where the C III line is located at λ ≈ 96 nm instead of λ ≈ 97.7 nm.

Moreover, the OP project (and the OPCD v. 3.3 base) includes spectral lines only for transitions in which the final state has an effective quantum number of ${\mathit{n}}_{\mathrm{f}}^{\mathrm{\ast }}\mathrm{\le }{\mathit{n}}_{\mathrm{0}}^{\mathrm{\ast }}$, where ${\mathit{n}}_{\mathrm{0}}^{\mathrm{\ast }}\mathrm{=}\mathrm{10}$ or $\mathrm{2}{\mathit{n}}_{\mathit{i}}^{\mathrm{\ast }}$ whenever the effective quantum number of an initial state is ${\mathit{n}}_{\mathit{i}}^{\mathrm{\ast }}\mathrm{\le }\mathrm{5}$ or  >5, respectively, viz., all levels with ${\mathit{n}}^{\mathrm{\ast }}\mathit{>}{\mathit{n}}_{\mathrm{0}}^{\mathrm{\ast }}$ are assumed to be fully dissolved, cf. Seaton et al. (1994). This effect is clearly seen in Figs. 12 and 13, where the Lyman and Balmer series of H I are plotted. A more careful examination of Fig. 12 reveals that the pseudocontinuum of He II extends from the ionization limit to about 97.8 nm.

Despite these limitations, the continuum level of energy fluxes of radiation agrees well with the ATLAS12/MHD-Q models, cf. Figs. 8, 12 and 13.

Fig. 14 Panel a) Emerging flux of radiation, lgH, near the Lyman series of H I for the ATLAS12/MHD-Q and BG25000g400v2.uv.7 models, respectively. Panel b) The same as panel a), but near the Balmer limit of H I.

Fig. 15 Same as Fig. 14, but near the line series of He I, He II and the Paschen limit of H I (panel b)).

Fig. 16 Rosseland-mean opacities, κR [cm2 g-1], vs. lgT for LTE (dotted lines) and non-LTE (dashed lines) models calculated with the TLUSTY200 computing code for Teff = 25   000 K and lgg = 4.0 with different mixtures of elements compared with the ATLAS12/MHD-Q models (thick lines). The chemical compositions of elements are the following: H+He (panels a) and b)), H–Ca (panels c) and d)), H+He+Fe (panels e) and f)) and all elements considered in the codes (panels g) and h)). The right-hand side panels display the same models as the left-hand side ones, but only the outermost layers are shown.

Fig. 17 Panel a) Emerging fluxes of radiation, lgH, near the Lyman limit of H I from LTE (dashed line) and non-LTE (continuous line) models calculated with the TLUSTY200 computing code including H–Fe species. Panels b) The ratio of the spectra shown in panel a).

Fig. 18 Same as Fig. 17, but near the limit of the Balmer series of H I.

4.3. Remarks about non-LTE models

In this section some additional remarks about non-LTE models of atmospheres are made in the context of the aim of this paper. The LH models of B-type stars (Lanz & Hubeny 2007) explicitly include 46 ions of H, He, C, N, O, Ne, Mg, Al, Si, S, and Fe that are represented by about 53 000 individual levels grouped into 1127 superlevels. Line opacity includes about 39 000 lines of elements from H to S and many more lines of Fe dynamically selected from a list of about 5.65 million lines taken mainly from the Kurucz atomic base. Despite this, the LH models underestimate stellar opacities, as discussed in Sect. 2. How this influences the emerging flux of radiation is illustrated in Figs. 14 and 15, where the ATLAS12/MHD-Q data are compared with synthetic spectra calculated by Lanz & Hubeny that are distributed for general use through the TLUSTY website (cf. Sect. 1).

Clearly, there is a quantitatively similar behavior of H I lines near the Lyman limit (cf. Fig. 14a) in agreement with the fact that both models use the same occupation probability formalism for H I lines. Note that of the markedly lower continuum level of the LH spectrum in comparison with the ATLAS12/MHD-Q, however. The same is true for the spectral region in which the Balmer lines are present, cf. Fig. 14b.

Figure 15a displays the wavelength regions near the ionization limit of the He I line series. This ionization limit is clearly seen in the BG25000g400v2.vis.7 spectrum (lower curve in panel a) calculated by Lanz & Hubeny. Again, the upper spectrum corresponds to the ATLAS12/MHD-Q model. As mentioned above, the occupation probability formalism was implemented in ATLAS12/MHD-Q code not only for hydrogen but also for other species and then used in calculations of the stellar monochromatic opacities and the equation of state. Finally, Fig. 15b shows the wavelength region near the Paschen limit of H I populated by He I and He II line series.

It is also interesting to compare models of atmospheres calculated with the ATLAS12/MHD-Q and TLUSTY200 codes for individual chemical species. Examples of such a comparison are given in Fig. 16, where models included only H−He, H−Ca, H+He+Fe and all absorbers considered in the codes are shown. All species up to Z = 30 are always included in the evaluation of the equation of state. As one can see in Fig. 16, the Rosseland-mean opacities obtained from ATLAS12/MHD-Q code are higher than those obtained from the TLUSTY200 code. Generally, differences increase when metallic atoms are added to the calculations, but note the fairly good agreement between models that include absorption caused by exactly the same elements, cf. panel c of Fig. 16, where H, He and Fe only are included in the two codes. We conclude that more attention should be given to model atoms of C–Ca in calculations of non-LTE models. These elements contribute significantly to stellar opacity and play a dominant role in cooling the outermost layers of atmospheres (cf. panels b and d of Fig. 16).

A deeper insight into the problem follows from a comparison of LTE and non-LTE models calculated with the TLUSTY200 code with the same input parameters. The results are shown in Figs. 17 and 18. In both cases the spectra were taken from the last iteration during the constructions of the models. Evidently, the non-LTE effects are much weaker than Figs. 14 and 15 suggest. The ratio of non-LTE to LTE spectra plotted in panels b indicates that non-LTE effects are important for spectral lines (especially in the UV range) while only weak effects are present in the Balmer, Paschen, etc. continuum levels.

As pointed out by Lanz & Hubeny (2007), there are still several points that are only crudely approximated. Some of them result from the present lack of sufficient atomic data, but some of them might require constructing more detailed model atoms, or possibly also splitting some superlevels into individual levels. On the other hand, there are many parameters in the TLUSTY200 code with default values recommended by Hubeny & Lanz. Many of them are the result of a compromise between computer time, computer memory, etc. and the physics of the problem. These parameters can be verified by users when new models are calculated.

We also calculated LTE and non-LTE models using the TLUSTY200 computing code with an extended list of model atoms. In particular, Fig. 1 shows the effect of the increasing number of O V levels from 6 (as in LH models of B stars, cf. Lanz & Hubeny 2007) to 40 levels (as in LH models of O stars, cf. Lanz & Hubeny 2003). Adding Mg I, P IV–P V, S VI, Fe VI, Ni V, and Ni VI ions and changing the default parameters, that control the width of lines to improve coverage of the frequency interval corresponding to some bound-bound transitions improves the agreement with the ATLAS12/MHD-Q models, but it does not entirely solve the problem pointed out in Sect. 2. The lack of less abundant elements in the LH-type models possibly plays an important role. We plan to investigate non-LTE models in a more detailed manner in a future paper.

5. Conclusions

We quantitatively compared LTE and non-LTE models of atmospheres with stellar envelope structures calculated for OPAL and OP data. An inspection of the published models of atmospheres (Castelli & Kurucz 2003; Lanz & Hubeny 2003, 2007) revealed systematic differences between them and with envelope models for early-type stars (cf. Sect. 2). The Castelli & Kurucz (2003) models show a high excess of opacity at the bottom of the photospheres of hot stars. If the Castelli & Kurucz models are correct, the new opacity bump at lgT ≈ 5.06 should be considered in addition to the opacity bumps at lgT ≈ 4.6 (dominated by hydrogen and helium lines) and lgT ≈ 5.3 (the main metallic opacity bump, which plays a crucial role in stellar oscillations). It should be noted, however, that the opacity bump at lgT ≈ 5.06 does not occur either in the Lanz & Hubeny models or in the newly calculated models presented in this paper. It is clear that non-LTE models calculated by Lanz & Hubeny underestimate stellar opacities.

We calculated new LTE models of atmospheres (denoted as ATLAS12/MHD-Q) with the modified ATLAS12 computing code distributed in a basic version by Kurucz. The major modifications are the implementation of the occupation probability formalism (cf. Sect. 3.1) and the OP/IRON photoionization cross-sections (cf. Sect. 3.2). The occupation probability formalism is important for two aspects: the stellar monochromatic opacities near ionization limits of line series and the equation of state. All atomic data needed to calculate the equation of state were taken from the OP TOPbase project. Partition functions of all ions were calculated for the following elements: H−Si, S, Ar, and Ca. For iron group elements, Sc−Cu, we used the data recently updated by Kurucz, which include ionization stages from I to IX. Partition functions for higher ions of Fe were calculated using OP TOPbase atomic data. We found that both the Saha-type equation of state with the partition functions implemented by the occupation probability functions and the MHD-type equation of state give the same results for the ionizations fractions and thermodynamic variables in a broad range of physical parameters. In models studied in this paper, the ATLAS/KPF (cf. Sect. 3.1) gives almost the same results. Implementing the OP/IRON resonance-averaged photoionization cross-sections into the ATLAS12/MHD-Q code improves the compatibility of the atmospheric models with the internal structure of stars, but differences remain.

More complete model atoms (cf. Kurucz 2011) than used in this paper (cf. Sect. 4.1) are necessary to obtain better agreement with the CK models. Probably, the Castelli & Kurucz (2004) approach to autoionizing transitions also plays important role, as our test calculations for Fe VIII suggest.

Models based on the ATLAS/OP code fit the internal structures very well because the same opacity and EOS data were used. However, these models are unable to fit observed spectra at high spectral resolution mode. Better sampling of OP data over $\mathit{u}\mathrm{=}\frac{\mathit{h\nu }}{\mathit{kT}}$, T and N_e is necessary to improve the accuracy of emerging spectra of ATLAS/OP models. Of course, synthetic spectra can be recalculated for ATLAS/OP models using a more complete list of line transitions with the best known broadening mechanisms of bound-bound transitions, as mentioned in Sect. 1. In principle, the ATLAS/OP models have an important diagnostic value for fully line-blanketed models of atmospheres and for OP data.

Concerning the above-mentioned new opacity bump at lgT ≈ 5.06, there is a suggestion from the Kurucz atomic data that the current versions of OPAL and OP opacities may need to be upgraded at least by a more complete list of bound-bound transitions and possibly by a more accurate treatment of autoionizing transitions. In terms of the OP project, this means that the splitting of spectroscopic terms of atoms into individual levels may be necessary for calculations of photoionization cross-sections with resonance structures. The new opacity bump changes markedly stellar envelope models. The well-known Fe-bump is located at lgT ≈ 5.30 in the OPAL and OP data, cf. Sect. 2. The new additional opacity bump, extracted in this paper from the Castelli & Kurucz (2003) models, is located on the outer side of the Fe-bump and therefore affects the driving zone for stellar pulsations. This markedly changes the predicted spectra of unstable oscillation modes, which will be a subject of a separate paper.

As mentioned in Sect. 4.3, we also plan to investigate non-LTE models of atmospheres in a more detailed manner in a future paper, but some conclusions and suggestions can be given here. A direct comparison of results obtained from non-LTE (e.g. Lanz & Hubeny 2003, 2007) and LTE models based on ATLAS codes (like Castelli & Kurucz 2003, or ATLAS12/MHD-Q models) can not be attributed entirely to non-LTE effects. Differences between current LTE and non-LTE models of atmospheres for B stars reflect mostly opacity differences with the exception of the outermost layers where non-LTE effects play an important role. The best way to quantitatively account for non-LTE effects is to calculate both LTE and non-LTE models with the same computing code for the same model atoms. As one can see from panels b of Figs. 17 and 18,

spectral lines are markedly affected by non-LTE effects. To study the observed energy flux distributions of stars measured in absolute units, non-LTE results should be corrected to models with a more complete treatment of the line-blanketing effect, otherwise lower values of T_eff and lgg as well as higher angular diameters derived from observations in the Paschen continuum will be obtained for B stars.

Acknowledgments

The author thanks the referee, Dr. Robert Kurucz, for helpful comments and for a critical reading of the manuscript.

References

All Tables

All Figures

Fig. 1 Panel a) The Rosseland-mean opacities, κR    [cm2   g-1] , as a function of lgT for Castelli & Kurucz (2003) (marked as CK line) and Lanz & Hubeny (2007) (LH line) models of ${\mathit{T}}_{\mathrm{eff}}^{}\mathrm{=}\mathrm{25}\hspace{0.17em}\mathrm{000}$ K and lgg = 4.0 in comparison with OPAL (OPAL line) and OP (OP line) stellar envelope models. A non-LTE model calculated for a more complex model of O V ion (cf. Sect. 4.3) is also shown (T line). Panel b) The same as panel a), but on κR vs. lgR diagram. Panel c) The same as panel a), but for the radiation pressure gradients (in cgs units) as a function of lgT.

In the text

Fig. 2 Panel a) The Rosseland-mean opacities, κR    [cm2   g-1] , as a function of lgT for OPAL (thick line) and OP (thin line) models of Teff = 25   000 K and lgg = 4.0. The chemical abundance parameters X = 0.708 and Z = 0.0172 are assumed. Panel b) The ratio of OP to OPAL Rosseland-mean opacities are displayed for input parameters of (lgT and lgρ) taken from the OPAL model shown in panel a). Panel c) The same as panel b), but for gas pressures.

In the text

	Fig. 3 Panel a) The occupation probability function Wn (in logarithmic units) for the bound states of H I with the principal quantum number n for a typical physical parameters at τR ≈ 0.2 of a model of Teff = 25   000 K and lgg = 4.0. Explanation of the functions displayed is given in the text. Panel b) The same as panel a), but for He II.
In the text

Fig. 4 The internal partition functions (in logarithmic units) vs. lgT for ATLAS12/MHD-Q model of Teff = 15   000 K and lgg = 4.0 (thick lines) in comparison with Kurucz’s data shown as thin lines. The panels show H I, He I and II, and the first three ionization stages of C, N, O, Ne and Si, respectively. The lowest right-hand panel shows the same, but for Fe III–V. The ionization states of I–V are plotted by means of lines as explained in the right-hand side of the figure.

In the text

	Fig. 5 The ionization fractions (in logarithmic units) of species for the same models as in Fig. 4. Symbols used are the same as in Fig. 4.
In the text

Fig. 6 Panel a) The Rosseland-mean opacities, κR    [cm2   g-1] , as a function of lgT for ATLAS12/MHD-Q (dotted line) in comparison with OP (continuous line) and OPAL (dashed line) envelope models of Teff = 25   000 K and lgg = 4.0. Panel c) The same as panel a), but OP/IRON bound-free absorptions of metals are included in the ATLAS12/MHD-Q model (dotted line). Panel e) The same as panelc), but OP opacities are used for hydrogen and helium. Panel g) The same as panel a), but the ATLAS12/OP model is plotted by means of an A/OP line. The right-hand side panels show the same as the left-hand side ones, but for the diagrams of κR vs. lgR.

In the text

	Fig. 7 Same as Fig. 6, but for ${\mathit{T}}_{\mathrm{eff}}^{}\mathrm{=}\mathrm{45}\hspace{0.17em}\mathrm{000}$ K and lgg = 5.0. The CK model (plotted as A/CK line) is also shown in panel g).
In the text

Fig. 8 Panel a) Emerging flux of radiation, H  [erg   cm-2   s-1   nm-1] , near the Lyman series of H I for the ATLAS12/MHD-Q model of Teff = 25   000 K and lgg = 4.0. Panel b) The ratio of the emerging fluxes for models calculated with different description of the occupation probability formalism (see text). Panel c) The same as panel a), but for the Balmer series of H I. Panel d) The same as panel b), but for the Balmer series of H I.

In the text

	Fig. 9 Emerging flux of radiation, H  [erg   cm-2   s-1   nm-1] , near the Lyman series of H I (panel a)) of the model of ATLAS12/MHD-Q shown in Figs. 6a and b is compared with models calculated with the ATLAS/PFK code with different profiles for H I, He I and He II line absorption coefficients (see text).
In the text

	Fig. 10 Same as Fig. 9, but now for the emerging flux of radiation near the Balmer series of H I.
In the text

	Fig. 11 OP monochromatic cross-sections (in atomic units) of H, He, C, N, O and Fe elements for lgT = 4.35 K and lgNe = 15.0 corrected for abundance parameters of the elements, σkfk. Panel a) shows the cross-sections for H (thick line) and He (dashed line). Panel b) shows data for C, N and O elements in comparison with the dominant absorber of H. Panel c) shows the same as panel b), but for Fe.
In the text

	Fig. 12 Emerging flux of radiation, H  [erg   cm-2   s-1   nm-1] , near the Lyman series of H I for ATLAS/OP models shown in panels a) − d) for (Teff,lgg) = (20 000, 4.5), (28 000, 4.5), (250 000, 4.0), (25 000, 5.0), respectively.
In the text

	Fig. 13 Same as Fig. 12, but near the Balmer limit of H I.
In the text

	Fig. 14 Panel a) Emerging flux of radiation, lgH, near the Lyman series of H I for the ATLAS12/MHD-Q and BG25000g400v2.uv.7 models, respectively. Panel b) The same as panel a), but near the Balmer limit of H I.
In the text

	Fig. 15 Same as Fig. 14, but near the line series of He I, He II and the Paschen limit of H I (panel b)).
In the text

Fig. 16 Rosseland-mean opacities, κR [cm2 g-1], vs. lgT for LTE (dotted lines) and non-LTE (dashed lines) models calculated with the TLUSTY200 computing code for Teff = 25   000 K and lgg = 4.0 with different mixtures of elements compared with the ATLAS12/MHD-Q models (thick lines). The chemical compositions of elements are the following: H+He (panels a) and b)), H–Ca (panels c) and d)), H+He+Fe (panels e) and f)) and all elements considered in the codes (panels g) and h)). The right-hand side panels display the same models as the left-hand side ones, but only the outermost layers are shown.

In the text

	Fig. 17 Panel a) Emerging fluxes of radiation, lgH, near the Lyman limit of H I from LTE (dashed line) and non-LTE (continuous line) models calculated with the TLUSTY200 computing code including H–Fe species. Panels b) The ratio of the spectra shown in panel a).
In the text

	Fig. 18 Same as Fig. 17, but near the limit of the Balmer series of H I.
In the text