Method 1

In our first approach, Compton scattering is taken into account in the radiation transfer equation using the Kompaneets operator (Kompaneets, 1957):

$$\frac{\partial^2 f_{\nu} J_{\nu}}{\partial \tau_{\nu}^2} = \frac{... ...- \frac{\sigma_{\rm e}}{k_{\nu}+\sigma_{\rm e}} \frac{kT}{m_{\rm e} c^2} \times$$ (3)

$$x \frac{\partial}{\partial x} \left(x \frac{\partial J_{\nu}}{\pa... ...frac{T_{\rm eff}}{T} x J_{\nu} \left( 1 + \frac{CJ_{\nu}}{x^3} \right) \right),$$

where $x=h \nu /kT_{\rm eff}$ is the dimensionless frequency, $f_{\nu}(\tau_{\nu}) \approx 1/3$ is the variable Eddington factor, $J_{\nu}$ is the mean intensity of radiation, $B_{\nu}$ is the black body (Planck) intensity, $T$ is the local electron temperature, $T_{\rm eff}$ is the effective temperature of WD, and $C=c^2 h^2~/~2(kT_{\rm eff})^3$. The optical depth $\tau_{\nu}$ is defined as

$$d \tau_{\nu} = (k_{\nu}+\sigma_{\rm e}) \, dm.$$ (4)

These equations have to be completed by the energy balance equation

$$\int_0^{\infty} k_{\nu}\left(J_{\nu} - B_{\nu}\right) d\nu - \sigma_{\rm e} \frac{kT}{m_{\rm e} c^2} \times$$ (5)

$$\left( 4 \int_0^{\infty} J_{\nu} \, d\nu - \frac{T_{\rm eff}}{T} \int_0^{\infty} x J_{\nu} (1+\frac{CJ_{\nu}}{x^3}) \, d\nu \right)=0,$$

the ideal gas law

$$P_{\rm g} = N_{\rm tot} kT,$$ (6)

where $N_{\rm tot}$ is the number density of all particles, and also by the particle and charge conservation equations. We assume local thermodynamical equilibrium (LTE) in our calculations, so the number densities of all ionisation and exitation states of all elements have been calculated using Boltzmann and Saha equations.

For solving the above equations and computing the model atmosphere we used a version of the computer code ATLAS (Kurucz, 1993; Kurucz, 1970), modified to deal with high temperatures; see Ibragimov et al. (2003) and Swartz et al. (2002) for further details. This code was also modified to account for Compton scattering.

The scheme of calculations is as follows. First of all, the input parameters of the WD are defined: the effective temperature $T_{\rm eff}$ and surface gravity $g=GM_{\rm wd}/R^2_{\rm wd}$. Then a starting model using a grey temperature distribution is calculated. The calculations are performed with a set of 98 depth points $m_{\rm i}$ distributed logarithmically in equal steps from $m\approx 10^{-6}$ g cm$^{-2}$ to $m_{\rm max}$. The appropriate value of $m_{\rm max}$ is found from the condition $\sqrt{\tau_{\nu,\rm b-f,f-f}(m_{\rm max})\tau_{\nu}(m_{\rm max})} >$ 1 at all frequencies. Satisfying this equation is necessary for the inner boundary condition of the radiation transfer.

For the starting model, all number densities and opacities at all depth points and all frequencies (we use 300 logarithmically equidistant frequency points) are calculated. The radiation transfer equation (3) is non-linear and is solved iteratively by the Feautrier method (Mihalas, 1978, see also Zavlin & Shibanov 1991; Pavlov et al. 1991; Grebenev & Sunyaev 2002). We use the last term of the equation (3) in the form $xJ_{\nu}^i(1+CJ_{\nu}^{i-1}/x^3)$, where $J_{\nu}^{i-1}$ is the mean intensity from the previous iteration. During the first iteration we take $J_{\nu}^{i-1}=0$. Between iterations we calculate the variable Eddington factors $f_{\nu}$ and $h_{\nu}$, using the formal solution of the radiation transfer equation in three angles at each frequency. Usually 2-3 iterations are sufficient to achieve convergence.

We used the usual condition at the outer boundary

$$\frac{\partial J_{\nu}}{\partial \tau_{\nu}} = h_{\nu} J_{\nu},$$ (7)

where $h_{\nu} \approx 1/2$ is surface variable Eddington factor. The inner boundary condition is

$$\frac{\partial J_{\nu}}{\partial \tau_{\nu}} = \frac{\partial B_{\nu}}{\partial \tau_{\nu}}.$$ (8)

The outer boundary condition is found from the lack of incoming radiation at the WD surface, and the inner boundary condition is obtained from the diffusion approximation $J_{\nu} \approx B_{\nu}$ and $H_{\nu} \approx 1/3 \times \partial B_{\nu}/\partial \tau_{\nu}$.

The boundary conditions along the frequency axis are

$$J_{\nu} = B_{\nu}$$ (9)

at the lower frequency boundary ( $\nu_{\rm min}=10^{12}$ Hz, $h\nu_{min}$ $\ll kT_{\rm eff}$) and

$$x \frac{\partial J_{\nu}}{\partial x} - 3J_{\nu} + \frac{T_{\rm eff}}{T} x J_{\nu} \left( 1 + \frac{CJ_{\nu}}{x^3} \right)=0$$ (10)

at the higher frequency boundary ( $\nu_{\rm max}\approx 3 \cdot 10^{17}$ Hz, $h\nu_{\rm max} \gg kT_{\rm eff}$). Condition (9) means that at the lowest energies the true opacity dominates the scattering $k_{\nu} \gg \sigma_{\rm e}$, and therefore $J_{\nu} \approx B_{\nu}$. Condition (10) means that there is no photon flux along frequency axis at the highest energy.

The solution of the radiative transfer equation (3) was checked for the energy balance equation (5) together with the surface flux condition

$$4 \pi \int_0^{\infty} H_{\nu} (m=0) d\nu = \sigma T_{\rm eff}^4 = 4 \pi H_0.$$ (11)

The relative flux error as a function of depth,

$$\varepsilon_{H}(m) = 1 - \frac{H_0}{\int_0^{\infty} H_{\nu} (m) d\nu},$$ (12)

was calculated, where $H_{\nu} (m)$ is radiation flux at given depth. This latter quantity was found from the first moment of the radiation transfer equation:

$$\frac{\partial f_{\nu} J_{\nu}}{\partial \tau_{\nu}} = H_{\nu}.$$ (13)

Temperature corrections were then evaluated using three different procedures. The first is the integral $\Lambda$-iteration method, modified for Compton scattering, based on the energy balance equation (5). It is valid in the upper atmospheric layers. The second one is the Avrett-Krook flux correction, which uses the relative flux error, and is valid in the deep layers. And the third one is the surface correction, which is based on the emergent flux error. See Kurucz (1970) for a detailed description of the methods.

The iteration procedure is repeated until the relative flux error is smaller than 1%, and the relative flux derivative error is smaller than 0.01%. As a result of these calculations, we obtain the self-consistent WD model atmosphere together with the emergent spectrum of radiation.

Our method of calculation was tested on a model of bursting neutron star atmospheres (Pavlov et al., 1991; Madej, 1991), and a model DA white dwarf atmosphere with $T_{\rm eff}=2\cdot 10^5$ K, $\log~g = 6.3$ (Madej, 1994). Agreement with the earlier calculations is extremely good. We illustrate the emergent spectrum from the latter calculation in Fig.1.

Figure 1: Spectra (a) and temperature structures (b) of a hot DA white dwarf model atmosphere with $T_{\rm eff}= 200\, 000$ K and $\log~g = 6.3$ with (solid line) and without (dashed line) Compton scattering taken into account.