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. Hauschildt1,2 - E. Baron3,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 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  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  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 ( ) 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  operator, .

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

 (1)

with
 ar = (2) = (3) = (4)

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

 (5)

where ds is a line element along a (curved) characteristic, Il(s) is the specific intensity along the characteristic at point  (s=0denotes the beginning of the characteristic) and wavelength point . The coefficient al is defined by

where , and r is the radius. and  are the emission and extinction coefficients at wavelength , respectively.

##### 2.1.1 First 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:

so that the wavelength-discretized SSRTE becomes

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

and introduce the source function, , we have

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):
 = (6) (7)

where we have suppressed the index labeling the ray; denotes the optical depth along the ray with  and  while  is calculated using piecewise linear interpolation of  along the ray, viz.

 (8)

The source function along a characteristic is interpolated by linear or parabolic polynomials so that
 (9)

is the optical depth along the ray from point i to point i+1. The coefficients , , and  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 I1 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 I1 is given as the outer boundary condition and (ii) integration from point N+2 to point 2N, where IN+1 is given as the inner boundary condition.

##### 2.1.2 Second discretization

The discretization of the derivative in the previous section can be deferred. We first rewrite Eq. (5) as

 (10)

with

and the definition of the CMF optical depth

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

with
 = (11) = (12)

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

 (13)

with the definitions

and

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

Note that the integration of 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

 = (14) = (15)

Here, and in the following, we use a sorted wavelength grid with 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 in the form

where

The formal solution then assumes the form

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

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

We describe the construction of 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  . To compute row j of the discrete -operator (or -matrix), , we sequentially label the intersection points of the characteristic kwith the shell i and define auxiliary quantities  and  as follows:
for i<j-1

for i=j-1

for i=j

for i=j+1

for

For the calculation of , we obtain:
for i = k+2

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

for i = k+j+2

for i = k+j+3

for i = k+j+4

for

Using the and , we can now write the -Matrix as

where wki,j are the angular quadrature weights, is the set and  is the set  . This expression gives the full -matrix, it can easily be specialized to compute only certain bands of the -matrix. In that case, not all of the and  have to be computed, reducing the CPU time from that required for the computation of the full -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  so that for the first scheme in Hauschildt (1992) we have and for the second discretization we have  . With this the optical depth scale along the ray d and  become

and is given by

whereas becomes

With this we can arbitrarily combine the two discretizations by varying  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 scheme handles the  and  parts differently since the (known)  is treated like a source term whereas the  term leads to a changed definition of d in the CMF along the ray. Therefore, the  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  and  will differ most).

This is not the case for the discretization which assumes that the 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 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 . 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  with . The line of the 2-level atom is parameterized by its strength relative to the continuum  and its thermalization fraction . For the intrinsic line profile we use a Gaussian profile with a of  . For the radial structure, we set and use a radial grid with  and  (so that the extension of the test atmosphere is 100) and a logarithmic optical depth grid from  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  . 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 ( ). Top part: CMF H in arbitrary units (full line: , dotted line: ), 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  and  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 ( ) and scattering dominated ( ) 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 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 discretization produces the expected flat continuum (once the initial conditions in wavelength are "forgotten''), the  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 ( ). Top part: CMF H in arbitrary units (full line: , dotted line: ), 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 ( ). Top part: CMF H in arbitrary units (full line: , dotted line: ), 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 ( ). Top part: CMF H in arbitrary units (full line: , dotted line: ), 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  . The discretization produces a perfectly flat continuum after the effect of the (grey) initial condition dies away. In contrast, the discretization produces considerable numerical diffusion in this test, on the average the solution deviates by more than 10% from the 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 ( ). Top part: CMF H in arbitrary units (full line: , dotted line: ), 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  was used. Open with DEXTER

It is clear that the discretization gives significantly better results than the 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 with 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 discretization scheme produces an emission feature that is about 2-3% stronger than the 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 scattering dominated line with an absorptive background continuum ( ). Top part: CMF H in arbitrary units (full line: , dotted line: ), bottom part: relative differences between CMF H (full line) and CMF J (dotted line) for the two discretizations. Open with DEXTER

 Figure 7: Strong scattering dominated line with a scattering dominated continuum ( ). Top part: CMF H in arbitrary units (full line: , dotted line: , dashed line:  ), bottom part: relative differences between CMF H for the two discretizations with and . In addition, the dotted line shows the relative differences between calculations with and . Open with DEXTER

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

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

#### 3.3 Numerical instability for

In certain conditions, it is possible for the 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 ( ) with a maximum speed of  , the inner boundary is at  . 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 , 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 ( ) and reaches a scattering optical depth . 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 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 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  as shown in the Fig. 8. The differences between the spectra for and 1.0 are very small, the case shows a few percent difference compared to . It is clear that the instability of the case can easily be avoided by using due to the stabilizing features of the discretization. In our test case even is sufficient to prevent the instability.

 Figure 8: Example to illustrate the possibility of numerical instabilities in the 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 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 , which in all our tests leads to recovered stability and is more accurate than a higher value of . However, until we determine the exact conditions that trigger the instability the cautious user should vary  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.