A&A 417, 317-324 (2004)
DOI: 10.1051/0004-6361:20034473

Improved discretization of the wavelength derivative term in CMF operator splitting numerical radiative transfer

P. H. Hauschildt^1,2 - E. Baron^3,4,5

1 - Hamburger Sternwarte, Gojenbergsweg 112, 21029 Hamburg, Germany
2 - Dept. of Physics and Astronomy & Center for Simulational Physics, University of Georgia, Athens, GA 30602-2451, USA
3 - Dept. of Physics and Astronomy, University of Oklahoma, 440 W. Brooks, Rm 131, Norman, OK 73019-0225, USA
4 - NERSC, Lawrence Berkeley National Laboratory, MS 50F-1650, 1 Cyclotron Rd, Berkeley, CA 94720-8139, USA
5 - Laboratoire de Physique Nucléaire et de Haute Énergies, CNRS-IN2P3, University of Paris VII, 75005 Paris, France

Received 8 October 2003 / Accepted 17 December 2003

Abstract
We describe two separate wavelength discretization schemes that can be used in the numerical solution of the comoving frame radiative transfer equation. We present an improved second order discretization scheme and show that it leads to significantly less numerical diffusion than the previous scheme. We also show that due to the nature of the second order term in some extreme cases it can become numerically unstable. We stabilize the scheme by introducing a mixed discretization scheme and present the results from several test calculations.

Key words: radiative transfer - methods: numerical

1 Introduction

The numerical solution of the radiative transfer equation plays a large role in our interpretation of spectroscopic data of astrophysical sources. New methods and faster computers have led to a resurgence of interest in solving the transfer equation (see for example Hubeny et al. 2003).

The numerical radiative transfer in the co-moving frame (CMF) method discussed in Hauschildt (1992) uses a discretization of the $\partial \lambda I/\partial \lambda$ terms in the RTE to obtain a formal solution. We show that a second order discretization scheme for the wavelength derivative leads to better numerical accuracy for a number of applications, such as stellar winds. In addition, the new discretization allows us to include the effects of the wavelength derivative in the construction of the approximate $\Lambda$ operator used in the operator splitting method. This improves the computational performance of the algorithm. It is possible to mix the two discretization scheme to tailor the performance of the algorithm to the problem being considered. In the following we describe the new discretization and the construction of the $\Lambda^*$ matrix for arbitrary bandwidths and discuss some test results. We have implemented this scheme into the PHOENIX code (see for example Hauschildt & Baron 1999).

2 Method

In the following discussion we use the same notation as in Hauschildt (1992) and reproduce the key equations for convenience. We use the spherically symmetric form of the special relativistic, time independent ( $\partial/\partial t \equiv 0$) RTE (hereafter, SSRTE), the restriction to plane parallel geometry is straightforward. The calculation of the characteristics is identical that of Hauschildt (1992) and we thus assume that the characteristics are known. First, we will describe the process for the formal solution, then we will describe how we construct the approximate $\Lambda$ operator, $\Lambda^*$.

2.1 Radiative transfer equation

In the CMF the SSRTE in the wavelength ($\lambda$) scale is given by (Mihalas & Weibel-Mihalas 1984; Hauschildt 1992)

 

$$a_r {\partial I \over \partial r} + a_\mu {\partial I \over \part... ...ambda{\partial \lambda I \over \partial \lambda} + 4a_\lambda I = \eta - \chi I$$ (1)

with

$$\gamma(\mu+\beta) ,$$                                    a_r = (2)

$$a_\mu \gamma(1-\mu^2) \left[{1+\beta\mu\over r} - \gamma^2\left(\mu+\beta\right){\partial \beta \over \partial r} \right],$$ = (3)

$$a_\lambda \gamma \left[ {\beta(1-\mu^2)\over r} +\gamma^2\mu\left(\mu+\beta\right){\partial \beta \over \partial r} \right]\cdot$$ = (4)

Along the characteristics, Eq. (1) has the form (Mihalas 1980)

$$\frac{{\rm d}I_l}{{\rm d}s} + a_l\frac{\partial \lambda I}{\partial \lambda} = \eta_l - (\chi_l+4a_l)I_l$$ (5)

where ds is a line element along a (curved) characteristic, I_l(s) is the specific intensity along the characteristic at point $s\ge 0$ (s=0denotes the beginning of the characteristic) and wavelength point $\lambda_l$. The coefficient a_l is defined by

$$a_l = \gamma \left[ \frac{\beta(1-\mu^2)}{r} +\gamma^2\mu\left(\mu+\beta\right){\partial \beta \over \partial r} \right]$$

where $\beta=v/c$, $\gamma=\sqrt{1-\beta^2}$ and r is the radius. $\eta_l$ and $\chi_l$ are the emission and extinction coefficients at wavelength $\lambda_l$, respectively.

2.1.1 First ${\partial I/\partial \lambda }$ discretization

As in Mihalas & Weibel-Mihalas (1984) and Hauschildt (1992), we discretize the wavelength derivative in the SSRTE with a fully implicit method in order to ensure stability:

$$\left.{\partial \lambda I \over \partial \lambda}\right\vert ... ...lambda_{l-1} I_{\lambda_{l-1}}\over \lambda_l - \lambda_{l-1}}$$

so that the wavelength-discretized SSRTE becomes

$${{\rm d}I_{\lambda_l}\over {\rm d}s} + a_\lambda {\lambda_l... ...\eta_{\lambda_l} - (\chi_{\lambda_l}+4a_\lambda)I_{\lambda_l}.$$

If we now define the optical depth scale along the ray as

$${\rm d}\tau \equiv \chi +a_\lambda \left(4+{\lambda_l\over \lambda_l-\lambda_{l-1}}\right) \equiv \hat\chi {\rm d}s,$$

and introduce the source function, $S=\eta/\chi$, we have

$${{\rm d} I\over {\rm d}\tau} = I - {\chi\over \hat\chi} ... ..._l-\lambda_{l-1}} I_{\lambda_{l-1}}\right) \equiv I-\hat S .$$

With this definition, the formal solution of the SSRTE along the characteristics can be written in the following way (cf. Olson & Kunasz 1987, for a derivation of the formulae):

$$I(\tau_i) I(\tau_{i-1}) \exp(\tau_{i-1}-\tau_i)$$ = (6)

$$+\int_{\tau_{i-1}}^{\tau_i} \hat S(\tau) \exp(\tau_{i-1}-\tau) ~ {\rm d}\tau$$

$$I(\tau_i) \equiv I_{i-1}\exp(-\Delta\tau_{i-1})+\Delta I_i$$ (7)

where we have suppressed the index labeling the ray; $\tau_i$ denotes the optical depth along the ray with  $\tau_1\equiv 0$ and  $\tau_{i-1} \le \tau_i$while $\tau$ is calculated using piecewise linear interpolation of $\hat\chi$ along the ray, viz.

 

$$\Delta\tau_{i-1} = (\hat\chi_{i-1}+\hat\chi_i)\vert s_{i-1}-s_i\vert/2.$$ (8)

The source function $\hat S(\tau)$ along a characteristic is interpolated by linear or parabolic polynomials so that

$$\Delta I_i = \alpha_i \hat S_{i-1} + \beta_i \hat S_i + \gamma_i \hat S_{i+1},$$ (9)

$\Delta\tau_i \equiv \tau_{i+1}-\tau_i$ is the optical depth along the ray from point i to point i+1. The coefficients $\alpha_i$, $\beta_i$, and $\gamma_i$ are given in Olson & Kunasz (1987).

For the tangent rays (for example ray 1 in Fig. 1 of Hauschildt 1992), the formal solution starts at point 2 with I₁ given as the outer boundary condition and proceeds along the ray. The formal solution for a core-intersecting ray is split into two parts: (i) integration from point 1 to point N, where I₁ is given as the outer boundary condition and (ii) integration from point N+2 to point 2N, where I_N+1 is given as the inner boundary condition.

2.1.2 Second ${\partial I/\partial \lambda }$ discretization

The discretization of the ${\partial I/\partial \lambda }$ derivative in the previous section can be deferred. We first rewrite Eq. (5) as

 

$$\frac{{\rm d}I_l}{{\rm d}\tau} - \frac{a_l}{\hat\chi_l}\frac{\partial \lambda I}{\partial \lambda} = I-\frac{\eta_l}{\hat\chi_l}$$ (10)

with

$$\hat\chi_l = \chi_l+4a_l$$

and the definition of the CMF optical depth

$${\rm d}\tau \equiv -\hat\chi_l {\rm d}s$$

along the characteristic. To obtain an expression for the formal solution, we rewrite Eq. (10) as

$$\frac{{\rm d}I}{{\rm d}\tau} = I - \hat S - \tilde S$$

with

$$\hat S \frac{\chi}{\hat\chi_l}S = \frac{\eta_l}{\hat\chi_l}$$ = (11)

$$\tilde S -\frac{a_l}{\hat\chi}\frac{\partial \lambda I}{\partial \lambda}\cdot$$ = (12)

With this we obtain the following expression for the formal solution (see also Eq. (20) in Hauschildt 1992)

$$I_{i,l}= I_{i-1,l}\exp(-\Delta\tau_{i-1}) + \delta\hat I_{i,l}+ \delta\tilde I_{i,l}$$ (13)

with the definitions

$$\delta\hat I_{i,l}= \alpha_{i,l}\hat S_{i-1,l}+ \beta_{i,l}\hat S_{i,l}+ \gamma_{i,l}\hat S_{i+1,l}$$

and

$$\delta\tilde I_{i,l}= \tilde\alpha_{i,l}\tilde S_{i-1,l}+ \tilde\beta_{i,l}\tilde S_{i,l}.$$

The index i labels the (spatial) points along a characteristic, the index l denotes the wavelength point. The coefficients  $\alpha_{i,l}$, $\beta_{i,l}$, and  $\gamma_{i,l}$ are given in Hauschildt (1992) and Olson & Kunasz (1987), here they are calculated for a fixed wavelengths for all points along a characteristic. $\hat S$ is a vector of known quantities (the old mean intensities and thermal sources). The $\tilde S$ contain the effects of the velocity field on the formal solution and are given by

$$\begin{eqnarray*}\tilde S_{i-1,l}& = & -\frac{a_{i-1,l}}{\hat\chi_{i-1,l}} \lef... ...{\partial \lambda I}{\partial \lambda}\right\vert _{_{i,l}}\cdot \end{eqnarray*}$$

Note that the integration of $\tilde S$ is performed using linear elements in order to allow for a recursive formal solution. Parabolic elements can be used but require the solution of matrix equations to obtain the formal solution.

If the velocity field is monotonically increasing or decreasing we use a stable upwind discretization of the wavelength derivative. In both cases, the problem becomes an initial value problem and can be solved for each wavelength once the results of the previous (smaller or longer) wavelength points are known. For a monotonically increasing velocity field this gives

$$\tilde S_{i-1,l} -\frac{a_{i-1,l}}{\hat\chi_{i-1,l}} \left[ \frac{\lambda_{l}}{\la... ...1}} I_{i-1,l}- \frac{\lambda_{l-1}}{\lambda_l-\lambda_{l-1}} I_{i-1,l-1}\right]$$ = (14)

$$\tilde S_{i,l} -\frac{a_{i,l}}{\hat\chi_{i,l}} \left[ \frac{\lambda_{l}}{\lambda... ...l-1}} I_{i,l} - \frac{\lambda_{l-1}}{\lambda_l-\lambda_{l-1}} I_{i,l-1}\right].$$ = (15)

Here, and in the following, we use a sorted wavelength grid with $\lambda_{l-1}<\lambda_l<\lambda_{l+1}$ for monotonically increasing velocities (the order of the grid would be reversed for monotonically decreasing velocity fields). For non-monotonic velocity fields this is no longer the case and the formal solution needs to explicitly account for the wavelength couplings (Baron & Hauschildt 2003).

With these formulae we can write $\delta\tilde I_{i,l}$ in the form

$$\begin{eqnarray*}\delta\tilde I_{i,l}&=& \tilde\alpha_{i,l}\left(p_{i-1,l}I_{i-1... ...\tilde\beta_{i,l}\left(p_{i,l}I_{i,l}- p_{i,l-1}I_{i,l-1}\right) \end{eqnarray*}$$

where

$$\begin{eqnarray*}p_{i-1,l}& = & -\frac{\lambda_{l}}{\lambda_l-\lambda_{l-1}}\fra... ...}{\lambda_l-\lambda_{l-1}}\frac{a_{i,l}}{\hat\chi_{i,l}}\cdot\\ \end{eqnarray*}$$

The formal solution then assumes the form

$$\begin{eqnarray*}\left(1-\tilde\beta_{i,l}p_{i,l}\right)I_{i,l}&=& \left(\tilde... ...\\ &&- \tilde\beta_{i,l}p_{i,l-1}I_{i,l-1} +\delta\hat I_{i,l}. \end{eqnarray*}$$

This equation can be solved recursively along a characteristic to calculate the I_i,l for all i and fixed l. The form of the wavelength discretization is now second order accurate in  $\Delta\lambda$.

The new form of the formal solution requires only small changes to the construction of the $\Lambda^*$ operator discussed in Hauschildt (1992) and Hauschildt et al. (1994).

We describe the construction of $\Lambda^*$ for arbitrary bandwidth using the example of a tangential characteristic. The intersection points (including the point of tangency) are labeled from left to right, the direction in which the formal solution proceeds. For convenience, we label the characteristic tangent to shell k+1 as k. Therefore, the characteristic k has 2k+1 points of intersection with discrete shells  $1\ldots k+1$. To compute row j of the discrete $\Lambda$-operator (or $\Lambda$-matrix), $\Lambda_{ij}$, we sequentially label the intersection points of the characteristic kwith the shell i and define auxiliary quantities  $\lambda_{ij}^k$and  $\hat\lambda_{ij}^k$ as follows:
for i<j-1

$$\lambda^k_{i,j} = 0$$

for i=j-1

$$\lambda^k_{j-1,j} = \left(1-\tilde\beta^k_{j-1} p^k_{j-1,l}\right)^{-1}\gamma^k_{j-1}$$

for i=j

$$\begin{eqnarray*}\lambda^k_{j,j} & = & \left(1-\tilde\beta^k_{j} p^k_{j,l}\right... ...{j-1,l}+\exp(-\Delta\tau^k_{j-1})\right] + \beta^k_{j} \right) \end{eqnarray*}$$

for i=j+1

$$\begin{eqnarray*}\lambda^k_{j+1,j} & = & \left(1-\tilde\beta^k_{j+1} p^k_{j+1,l... ...^k_{j,l}+\exp(-\Delta\tau^k_{j})\right] + \alpha^k_{j+1} \right) \end{eqnarray*}$$

for $j+1 < i \le k+1$

$$\begin{eqnarray*}\lambda^k_{i,j} & = & \left(1-\tilde\beta^k_{i} p^k_{i,l}\right... ...ha^k_{i} p^k_{i-1,l}+\exp(-\Delta\tau^k_{i-1})\right] \right). \end{eqnarray*}$$

For the calculation of $\hat\lambda^k_{i,j}$, we obtain:
for i = k+2

$$\begin{eqnarray*}\hat\lambda^k_{i,j} & = & \left(1-\tilde\beta^k_{i} p^k_{i,l}\r... ...pha^k_{i} p^k_{i-1,l}+\exp(-\Delta\tau^k_{i-1})\right] \right) \end{eqnarray*}$$

for k+2 <i < k+j+2

$$\begin{eqnarray*}\hat\lambda^k_{i,j} & = & \left(1-\tilde\beta^k_{i} p^k_{i,l}\r... ...ha^k_{i} p^k_{i-1,l}+\exp(-\Delta\tau^k_{i-1})\right] \right) \end{eqnarray*}$$

for i = k+j+2

$$\begin{eqnarray*}\hat\lambda^k_{i,j} & = & \left(1-\tilde\beta^k_{i} p^k_{i,l}\r... ...i-1,l}+\exp(-\Delta\tau^k_{i-1}) \right] + \alpha^k_i \right) \end{eqnarray*}$$

for i = k+j+3

$$\begin{eqnarray*}\hat\lambda^k_{i,j} & = & \left(1-\tilde\beta^k_{i} p^k_{i,l}\r... ...{i-1,l}+ \exp(-\Delta\tau^k_{i-1}) \right] + \beta_i \right) \end{eqnarray*}$$

for i = k+j+4

$$\begin{eqnarray*}\hat\lambda^k_{i,j} & = & \left(1-\tilde\beta^k_{i} p^k_{i,l}\r... ...i-1,l}+\exp(-\Delta\tau^k_{i-1}) \right] + \gamma^k_i \right) \end{eqnarray*}$$

for $k+j+5 \le i \le 2k+1$

$$\begin{eqnarray*}\hat\lambda^k_{i,j} & = & \left(1-\tilde\beta^k_{i} p^k_{i,l}\r... ...ha^k_{i} p^k_{i-1,l}+\exp(-\Delta\tau^k_{i-1})\right] \right). \end{eqnarray*}$$

Using the $\lambda^k_{ij}$ and $\hat\lambda^k_{ij}$, we can now write the $\Lambda$-Matrix as

$$\Lambda_{ij} = \sum_k \left( \sum_{\{l\}} w^k_{l,j}\lambda^k_... ...\{l'\}} w^k_{2(k+1)-l',j}\hat\lambda^k_{2(k+1)-l',j} \right)$$

where w^k_i,j are the angular quadrature weights, $\{l\}$ is the set $\{i \le k+1\}$ and $\{l'\}$ is the set  $\{i > k+1\}$. This expression gives the full $\Lambda$-matrix, it can easily be specialized to compute only certain bands of the $\Lambda$-matrix. In that case, not all of the $\lambda^k_{i,j}$ and  $\hat\lambda^k_{i,j}$ have to be computed, reducing the CPU time from that required for the computation of the full $\Lambda$-matrix.

2.1.3 Stability

In the second discretization scheme, the wavelength derivative contains an explicit term, which is required to derive a recursive method with second order accuracy. We can stabilize the discretion via a method similar to the standard Crank-Nicholson scheme (Abramowitz & Stegun 1972). To combine the two discretization schemes described above we introduce a factor  $\zeta \in [0,1]$ so that for the first scheme in Hauschildt (1992) we have $\zeta \equiv 1$ and for the second discretization we have  $\zeta \equiv 0$. With this the optical depth scale along the ray d$\tau$ and $\hat\chi$ become

$${\rm d}\tau = \chi +a_\lambda \left(4+\zeta{\lambda_l\over \lambda_l-\lambda_{l-1}}\right) \equiv \hat\chi {\rm d}s,$$

and $\hat S$ is given by

$$\hat S = {\chi\over \hat\chi} \left(S + \zeta{a_\lambda\... ...{l-1}\over \lambda_l-\lambda_{l-1}} I_{\lambda_{l-1}}\right).$$

whereas $\tilde S$ becomes

$$\tilde S = -(1-\zeta)\frac{a_l}{\hat\chi}\frac{\partial \lambda I}{\partial \lambda}\cdot$$

With this we can arbitrarily combine the two ${\partial I/\partial \lambda }$discretizations by varying $\zeta$ over the interval [0,1]. This allows us to combine them to optimize accuracy of the overall solution of the SSRTE.

2.1.4 Discussion

It should be expected that the two discretization schemes introduced above will show different numerical behavior. The $\zeta = 1$ scheme handles the  $I_{\lambda_l}$ and  $I_{\lambda_{l-1}}$ parts differently since the (known)  $I_{\lambda_{l-1}}$ is treated like a source term whereas the $I_{\lambda_l}$ term leads to a changed definition of d$\tau$ in the CMF along the ray. Therefore, the  $I_{\lambda_l}$ is basically assumed to vary exponentially across a characteristic. This will introduce more diffusive behavior for large optical depths and strong scattering (where the intrinsic assumptions for  $I_{\lambda_l}$ and  $I_{\lambda_{l-1}}$ will differ most).

This is not the case for the $\zeta = 0$ discretization which assumes that the $\lambda\partial I/\partial \lambda$ itself can be fitted to a first order polynomial and treated as a source term. The discretization step is delayed until the integration of the source term. Although the discretization is implicit by design, the behavior of ${\partial I/\partial \lambda }$ is a priori unknown and the second order nature of the discretization leads to an explicit term. Under extreme conditions, e.g., a strong line in an optically thin rapidly expanding gas, it is possible that the numerical integration is dominated by the "known'' term due to large changes in ${\partial I/\partial \lambda }$. This leads to numerical instability. We therefore expect that the two discretizations will show different behavior in numerical tests.

3 Tests

In order to compare and test the two different discretization schemes, we have devised a simple test model. We consider the case of a 2-level atom with background continuum on a radial grid (assuming spherical symmetry). The background continuum is assumed to be grey with adjustable thermalization fraction  $\epsilon_{\rm c} = \kappa_{\rm c}/\chi_{\rm c}$ with $\chi_{\rm c}=\kappa_{\rm c}+\sigma_{\rm c}$. The line of the 2-level atom is parameterized by its strength relative to the continuum  $\chi_l/\chi_{\rm c}$ and its thermalization fraction $\epsilon_l = \kappa_l/\chi_l$. For the intrinsic line profile we use a Gaussian profile with a $\Delta \lambda_D$ of  $30\hbox{$\,$ km$\,$ s$^{-1}$ }$. For the radial structure, we set $\chi_{\rm c} \propto r^{-2}$ and use a radial grid with  $R_{\rm out}=101$ and  $R_{\rm in} = 1$ (so that the extension of the test atmosphere is 100) and a logarithmic optical depth grid from  $\tau = 10^4$ to 10^-6 in the continuum to fix the scale of the opacities. The velocity law is linear (homologous expansion) with a prescribed maximum velocity of  $1000\hbox{$\,$ km$\,$ s$^{-1}$ }$. The computations are performed on a wavelength grid that uses a stepsize of 0.1 line width inside the line and 5 times the line width outside the line if not specified in detail for some of the test calculations.

3.1 Continuum tests

	Figure 1: Absorptive continuum ( $\epsilon _{\rm c}=1$). Top part: CMF H in arbitrary units (full line: $\zeta = 0$, dotted line: $\zeta = 1$), bottom part: relative differences between CMF H (full line) and CMF J (dotted line) for the two discretizations. The + signs give the location of the actual wavelength points used in the computation.
Open with DEXTER

As a first test we verified that there are no differences between the solutions for $\zeta = 0$ and $\zeta = 1$ for zero expansion velocity. The next test is for a pure continuum case (i.e., line strength of zero). This is actually a difficult test for the method as the flat continuum should be reproduced by the algorithm. The results for purely absorptive ( $\epsilon _{\rm c}=1$) and scattering dominated ( $\epsilon _{\rm c}=10^{-2}$) continua are shown in Figs. 1 and 2, respectively. In this case the CMF H (the first Eddington moment of the intensity in the CMF) should be constant. From Fig. 1 we see that both discretization schemes deliver nearly identical results with differences below 0.3% in the case of an absorption dominated continuum. For the case of a scattering dominated continuum the differences are significantly larger, nearly 6%, cf. Fig. 2. From the top part of the figure it is obvious that the differences start to increase just around the region of resolution change in wavelength at the (zero strength) line. That scattering increases the diffusive effect of the $\zeta = 1$wavelength discretization as shown by comparing Figs. 1 and 2. This effect becomes larger and more pronounced if the wavelength resolution is lower, which is shown in Figs. 3 and 4. In this test, we have decreased the wavelength resolution by roughly a factor of 10. The case of the absorption dominated continuum shows similar behavior as the test with higher wavelength resolution with the relative differences increasing slightly. However; the scattering dominated case produces much larger relative differences of up to 13% between the two discretizations. Whereas the $\zeta = 0$ discretization produces the expected flat continuum (once the initial conditions in wavelength are "forgotten''), the $\zeta = 1$ discretization produces a pseudo line feature just due to the change in resolution. This illustrates the higher accuracy of the second order scheme, which since it is properly centered produces more accurate results on an irregular grid.

	Figure 2: Scattering dominated continuum ( $\epsilon _{\rm c}=10^{-2}$). Top part: CMF H in arbitrary units (full line: $\zeta = 0$, dotted line: $\zeta = 1$), bottom part: relative differences between CMF H (full line) and CMF J (dotted line) for the two discretizations. The + signs give the location of the actual wavelength points used in the computation.
Open with DEXTER

	Figure 3: Absorptive continuum ( $\epsilon _{\rm c}=1$). Top part: CMF H in arbitrary units (full line: $\zeta = 0$, dotted line: $\zeta = 1$), bottom part: relative differences between CMF H (full line) and CMF J (dotted line) for the two discretizations. The + signs give the location of the actual wavelength points used in the computation.
Open with DEXTER

	Figure 4: Scattering dominated continuum ( $\epsilon _{\rm c}=10^{-2}$). Top part: CMF H in arbitrary units (full line: $\zeta = 0$, dotted line: $\zeta = 1$), bottom part: relative differences between CMF H (full line) and CMF J (dotted line) for the two discretizations. The + signs give the location of the actual wavelength points used in the computation.
Open with DEXTER

The case of a constant step size wavelength grid is shown in Fig. 5. The step size in wavelength is  $300\hbox {$\,$ km$\,$ s$^{-1}$ }$. The $\zeta = 0$ discretization produces a perfectly flat continuum after the effect of the (grey) initial condition dies away. In contrast, the $\zeta = 1$ discretization produces considerable numerical diffusion in this test, on the average the solution deviates by more than 10% from the $\zeta = 0$ solution. In addition, the continuum is not flat but shows a noticeable increase in CMF H (and J) toward the red.

	Figure 5: Scattering dominated continuum ( $\epsilon _{\rm c}=10^{-2}$). Top part: CMF H in arbitrary units (full line: $\zeta = 0$, dotted line: $\zeta = 1$), bottom part: relative differences between CMF H (full line) and CMF J (dotted line) for the two discretizations. In this calculation equidistant wavelength grid with a spacing of  $300\hbox {$\,$ km$\,$ s$^{-1}$ }$ was used.
Open with DEXTER

It is clear that the $\zeta = 0$ discretization gives significantly better results than the $\zeta = 1$ approach for the continuum tests.

3.2 Line tests

In contrast to the previous set of tests, we now introduce a strong line. We set the line strength $\chi_l/\chi_{\rm c}=100$ with $\epsilon_l=10^{-4}$ to simulate a strong, scattering dominated line of a 2-level atom against a background continuum. Figures 6 and 7 show the results for absorption and scattering dominated background continua. In both cases the relative differences in the continua are comparable to the previous tests. The differences in the lines are actually comparatively small in both cases, significantly smaller than might be expected from the continuum test results considering the fact that the line is strongly scattering dominated. This is most likely caused by the smoothing effect of the rapidly varying opacity across the line. The effect of continuum scattering is to considerably widen the line due to line photons being scattered by the continuum. The more diffusive $\zeta = 1$ discretization scheme produces an emission feature that is about 2-3% stronger than the $\zeta = 0$ discretization relative to the continuum. It is interesting to note the strong effect of even a weak (compared to the line itself) background continuum on the shape of the line and the differences between the two discretizations.

	Figure 6: Strong $(\chi _l/\chi _{\rm c}=100)$ scattering dominated $(\epsilon _l=10^{-4})$ line with an absorptive background continuum ( $\epsilon _{\rm c}=1$). Top part: CMF H in arbitrary units (full line: $\zeta = 0$, dotted line: $\zeta = 1$), bottom part: relative differences between CMF H (full line) and CMF J (dotted line) for the two discretizations.
Open with DEXTER

Figure 7: Strong $(\chi _l/\chi _{\rm c}=100)$ scattering dominated $(\epsilon _l=10^{-4})$ line with a scattering dominated continuum ( $\epsilon _c=10^{-2}$). Top part: CMF H in arbitrary units (full line: $\zeta = 0$, dotted line: $\zeta = 1$, dashed line:  $\zeta = 0.5$), bottom part: relative differences between CMF H for the two discretizations with $\zeta = 0$ and $\zeta = 1$. In addition, the dotted line shows the relative differences between calculations with $\zeta = 0.5$ and $\zeta = 0$.

Open with DEXTER

For the final test we show in Fig. 7 he results of a calculation for $\zeta = 0.5$. Compared to $\zeta = 1$ this calculation shows much smaller differences from the $\zeta = 0$ case in the continuum. From the top portion of the figure it is clear that in the line itself the $\zeta = 0.5$ model falls roughly half way in between the $\zeta = 1$ and $\zeta = 0$ cases.

In all cases the $\zeta = 0$ discretization gives numerically better results. However, there are situations where the $\zeta = 0$discretization is numerically unstable, which we discuss in the following section.

3.3 Numerical instability for $\zeta = 0$

In certain conditions, it is possible for the $\zeta = 0$ discretization to become numerically unstable. This is illustrated in the test case shown in Fig. 8. The two test lines are in complete LTE (S=B), the background continuum is scattering dominated. The structure is taken from typical SN Ia structure, but we have completely ignored the strong line blanketing inherent in SNe Ia by including line opacity from only Ca II. This procedure is actually quite useful for line identification in detailed model calculations. The velocity field is homologous ( $v\propto r$) with a maximum speed of  $30~000\hbox{$\,$ km$\,$ s$^{-1}$ }$, the inner boundary is at  $v=2000\hbox{$\,$ km$\,$ s$^{-1}$ }$. We use this model rather than the simpler more academic test cases for two reasons: 1) this is the model where we actually discovered the instability for $\zeta = 0$, and 2) we have been unable to find a simple test case that so strongly reproduces the instabilities. The continuum is optically thin in absorption ( $\tau_{\rm abs} \approx 4\times 10^{-4}$) and reaches a scattering optical depth $\tau_{\rm scatt} \approx 1.5$. We use a Planck-function as the inner boundary condition for the intensities for simplicity, physically it is better to employ a nebular [symmetry] boundary condition which is shown in Fig. 9. We see that the instability is evident (for the $\zeta = 0$ case) even though the line has gone into emission (note that the different scales between Figs. 8 and 9 give the impression that the instability is somewhat suppressed in the nebular case, but it is simply due to the larger range required to plot the strong emission features). The unstable wiggles are not associated with the grid and varying the wavelength resolution does not significantly affect the output spectrum. While we believe from our tests that the instability is produced predominantly in optically thin atmospheres, we have not been able to ascertain the precise conditions that trigger the instability.

The instability in the $\zeta = 0$ discretization is obvious starting at the rest wavelength of the line. Even in the line trough the oscillations created by the instability are apparent. The overall line shapes are distorted, even in the red emission feature. The instability vanishes for  $\zeta\ne 0$ as shown in the Fig. 8. The differences between the spectra for $\zeta = 0.5$ and 1.0 are very small, the $\zeta=0.1$ case shows a few percent difference compared to $\zeta = 1$. It is clear that the instability of the $\zeta = 0$ case can easily be avoided by using $\zeta>0$ due to the stabilizing features of the $\zeta = 1$ discretization. In our test case even $\zeta=0.1$ is sufficient to prevent the instability.

	Figure 8: Example to illustrate the possibility of numerical instabilities in the $\zeta = 0$ discretization. The structure is taken from an type Ia supernova calculation, only Ca II lines are included in the opacity and the lines are assumed to be in complete LTE.
Open with DEXTER

	Figure 9: The same structure as was used in Fig. 8, but with nebular (symmetry) boundary conditions.
Open with DEXTER

4 Conclusions

We have studied the discretization of the $\frac{\partial}{\partial\lambda}$ term in the co-moving frame radiative transfer equation. We have found that we can write a second order discretization scheme, which is in general more accurate than our previous method, but it is possible for this scheme to become unstable. We have constructed a hybrid Crank-Nicholson like scheme, which is unconditionally stable. We have been unable to identify the exact conditions which lead to the unstable behavior, but we have shown that even with a small admixture of the unconditionally stable scheme, stability is recovered. We recommend using a small value of $\zeta=0.1$, which in all our tests leads to recovered stability and is more accurate than a higher value of $\zeta$. However, until we determine the exact conditions that trigger the instability the cautious user should vary $\zeta$ and determine the sensitivity of the results to its value.

Acknowledgements

We thank the referee for a very helpful report which significantly improved the presentation of this paper. This work was supported in part by NASA grants NAG 5-8425 and NAG 5-3619 to the University of Georgia and by NASA grant NAG5-3505, NSF grants AST-0204771 and AST-0307323, and an IBM SUR grant to the University of Oklahoma. PHH was supported in part by the Pôle Scientifique de Modélisation Numérique at ENS-Lyon. Some of the calculations presented here were performed at the San Diego Supercomputer Center (SDSC), supported by the NSF, at the National Energy Research Supercomputer Center (NERSC), supported by the U.S. DOE, and at the Höchstleistungs Rechenzentrum Nord (HLRN). We thank all these institutions for a generous allocation of computer time.

References

• Abramowitz, M., & Stegun, I. 1972, Handbook of Mathematical Functions (New York: Dover) (In the text)
• Baron, E., & Hauschildt, P. H. 2003, A&A, submitted (In the text)
• Hauschildt, P. H. 1992, JQSRT, 47, 433 (In the text)
• Hauschildt, P. H., & Baron, E. 1999, JCAM, 102, 41 (In the text)
• Hauschildt, P. H., Störzer, H., & Baron, E. 1994, JQSRT, 51, 875 (In the text)
• Hubeny, I., Mihalas, D., & Werner, K. 2003, Stellar Atmosphere Modeling (San Francisco: ASP), ASP Conf. Ser., 288 (In the text)
• Mihalas, D. 1980, ApJ, 237, 574 [CrossRef] (In the text)
• Mihalas, D., & Weibel-Mihalas, B. 1984, Foundations of Radiation Hydrodynamics (Oxford: Oxford University Press)
• Olson, G. L., & Kunasz, P. B. 1987, JQSRT, 38, 325 (In the text)