A&A 401, 419-428 (2003)
D. Korcáková1,2 - J. Kubát2
1 - Ústav teoretické fyziky a astrofyziky PrF MU, Kotlárská 2, CZ-611 37 Brno, Czech Republic
2 - Astronomický ústav, Akademie ved Ceské republiky, CZ-251 65 Ondrejov, Czech Republic
Received 21 May 2002 / Accepted 17 January 2003
We describe the solution of the one-dimensional radiative transfer equation using the discontinuous finite element method. Tests of this method prove its reliability for static media. The static solution is then applied to the solution of the transfer equation in a moving medium using a method which we call here the Local Lorentz Transformation (LLT) method. This new method is applicable also to non-monotonic velocity fields. Direct solution of the radiative transfer equation in the moving medium using the discontinuous finite element method is also presented for cases both with and without inclusion of the effect of aberration. All three techniques are compared and they give comparable results.
Key words: radiative transfer - methods: numerical
The solution of the equation of radiative transfer belongs to the most important problems in astrophysics. There are many methods and computer codes for static media, but there are also many problems that do not allow ignoring of the motion of a medium, like studying atmospheres of hot stars with fast stellar winds in O stars or Wolf-Rayet stars, expanding nebulae, novae, supernovae, extragalactic radio sources, and fast changes in solar prominences. Consequently, the problem of the transfer of radiation is more complicated. The reasons for this are the coefficients of emissivity and opacity, which in the case of differential motion are angle dependent.
For continuum radiation we can simply use the static equation of radiative transfer, since the opacity and emissivity coefficients change only a little within the region of the Doppler frequency shift. However, this is not true for lines. Changes of the opacity caused by the Doppler shift may be large, and we have to take it into consideration. Generally, there are two basic classes of methods for its solution, namely the solution in the observer frame (usually a frame connected with the center of a star) and the solution in the comoving frame, which is the frame connected with (usually outflowing) material. A basic review and discussion of the applicability of both classes of methods can be found in Mihalas & Kunasz (1986).
Solution of the radiative transfer equation in the comoving frame is a commonly used possibility. It takes advantage of the fact that the opacity and emissivity coefficients are not affected by motion and may be treated as in the static case (i.e., if they are isotropic for the static case then they are isotropic also in a local reference frame in a moving atmosphere). However, the transfer equation becomes more complicated due to an additional term containing the frequency derivative of the intensity. This term describes the Doppler shift of radiation. The comoving-frame transfer equation is usually solved using Feautrier variables, which introduce additional restrictions to monotonic velocity fields (see Mihalas et al. 1975,1976a,1976b,1976c; Mihalas & Kunasz 1978). Feautrier variables have also been used by Hillier (1987). In addition, Feautrier variables are not usable for the solution of the relativistic transfer equation, where characteristic rays are not symmetric (Mihalas 1980) and use of short-characteristics is more appropriate (Hauschildt 1992).
Another effect, which is usually neglected, is the effect of aberration. The importance of aberration was discussed by Mihalas et al. (1976b); they concluded that for stellar winds (where ) the aberration and advection effects may be safely ignored. This was in principle confirmed by Hauschildt et al. (1995) who found the aberration and advection effects for novae to be negligible. On the other hand, Hauschildt et al. (1991) showed that for , omission of aberration leads to errors larger than 50% and Baron et al. (1996) found non-negligible effects of advection for velocities larger than . Another study of the effects of aberration and advection was performed by Peraiah (1987, 1991a), who found significant differences in the mean intensity and flux for maximum velocities between and both for planar and spherical geometry. He also obtained changes in spectral lines (Peraiah 1991a), but it is not clear what the real effect of aberration is.
In the first part of this paper we discuss briefly the one-dimensional equation of radiative transfer, which we will use in the subsequent parts. Here we use the slightly overlooked method of discontinuous finite elements, which has been introduced into astronomy by Castor et al. (1992), and Dykema et al. (1996). The solution of the static radiative transfer equation using this method is described in Sect. 3. In Sect. 4 we describe two new methods of solution of the one-dimensional equation of radiative transfer in moving media in a comoving frame based on the discontinuous finite elements method. The solution of the radiative transfer equation where the effect of aberration is not neglected is presented as well. In the last part we discuss the advantages and disadvantages of the method of discontinuous finite elements in these cases.
If the velocity of the medium is not very high with respect
to the observer
frame, we can neglect aberration, because the influence of the Doppler
effect is much stronger.
In this case, Eq. (3) reduces to (see Mihalas
1978, Eq. (14.99))
In the past, several classes of methods for the solution of the radiative transfer equation have been used. Possibly the most successful class of methods that have been used in a number of NLTE model atmosphere codes (e.g. Mihalas et al. 1975; Werner 1986; Hubeny 1988; Kubát 1994, 2002) is the solution of the transfer equation as a differential equation using the Feautrier (1964) method. Another class of methods is based on an integral method (see, e.g., Kalkofen 1974; Avrett & Loeser 1984).
In this paper, we deal with the in astrophysical radiative transfer slightly overlooked method of discontinuous finite elements (DFE). This method was originally developed for a study of neutron transfer by Goldin (1964). It is very similar to the finite element method, which has been used very often in various calculations in physics in the last few decades (see, e.g., Rektorys 1980). It was at first used for a solution of the radiative transfer equation by Castor et al. (1992). Further improvements and generalization of this method to more spatial dimensions were published by Dykema et al. (1996). Recently, a finite element method was used by Richling et al. (2001) for the solution of the monochromatic static radiative transfer equation in three dimensions.
The discontinuous finite element method is a generalization of the finite element method (see, e.g., Rektorys 1980). This method is based on the Galerkin method. To understand the discontinuous finite element method, we first describe briefly the Galerkin method.
Let us take a differential equation
Appropriate base functions for the finite element method may be
found, e.g., according to Rektorys (1980).
We section the region, where we seek the solution u, to cells.
We consider the base function to be nonzero only in a single cell.
We choose the base functions to be as simple as possible, e.g.,
polynomials or splines.
Let us assume that the subspace in a particular cell is L-dimensional,
i.e. we have L nonzero base functions wl (
L=K/number of cells).
Then we solve the following system of algebraic equations for every
Contrary to the above described Galerkin method, in the discontinuous
finite element method we do not neglect the residue
In this section we describe the solution of Eq. (2). The method basically follows Castor et al. (1992). We present this simple case here, because we shall use this method later in a more complicated situation of moving atmosphere.
First, we introduce the discretization of the geometrical
Note that the subscript d=1 corresponds to the deepest point.
We assume a linear dependence of the intensity on the depth within
cell (see Fig. 1).
|I1,D = 0.||(22)|
|Figure 2: Comparison of the Feautrier method (dashed line) and the method of discontinuous finite elements (solid line) for a static planar atmosphere.|
In this section we describe the application of the method of discontinuous finite elements to the solution of the transfer equation with velocity fields, Eqs. (3) and (4).
This method is based on the fact that when we solve a radiative transfer equation in a medium, which moves with a constant uniform velocity with respect to a static observer, the static solution of the radiative transfer equation is fully sufficient. This is obviously not true if there is a velocity gradient present in the medium. However, if we split the differentially moving medium into cells with constant uniform velocities moving with respect to each other, the transfer equation inside these cells may still be considered a static one thus preserving the simplicity of a static radiative transfer Eq. (2) locally. We proceed as follows.
We split the region into cells as described in Sect. 3.1.
In every cell we assume a constant velocity of the medium.
The velocity changes only at cell boundaries.
For quantities that do not change within the cell, but which
change on cell boundaries (grid points) we use "half indexes'', e.g.,
which means that this cell is between grid points d and d+1.
Since we consider individual cells as local inertial frames, we must
transform both intensity and frequency at the boundaries between cells.
Transformation laws for these quantities are
In practice we use the same frequency points as in the static case. We solve the Eq. (2) for each cell for these frequency points. At the cell boundary, we calculate the Doppler shift of frequency, Eq. (24), and interpolate the intensity I2,d (Sect. 3.1 and Fig. 1) to these new "rest'' frequencies. We tested both linear and parabolic polynomial interpolation and found the difference at the sixth order only. Because the equation of radiative transfer is Lorentz invariant, we solve the static Eq. (2) in each cell using the DFE method described in Sect. 3.1, i.e. we solve the Eq. (17) with and given by (18). The resulting intensity is given by Eq. (20). We chose the geometrical depth as an independent variable. Due to the strong dependence of the opacity on frequency in a spectral line the geometrical depth is more suitable as an independent variable than the optical depth.
We tested this method for several particular cases.
Here we show the result of calculations for a stellar wind with a
velocity expressed by the beta law (Cassinelli & Lamers 1999,
|Figure 3: Comparison of H line profiles calculated for the static case using the simple DFE method (dashed line) and for a stellar wind with and using the Local Lorentz Transformation together with DFE method (full line).|
This method is most effective in a case where the velocity gradient is not "too large''. We simply have to guarantee that the sets of new and old frequency points at a cell boundary significantly overlap. More exactly, the frequency shift must be smaller than the Doppler halfwidth. With increasing velocity gradient the number of cells has also to increase and, consequently, one needs more computing time to solve the equation. However, for "very'' large velocity gradients we may use the Sobolev method.
The local Lorentz transformation method is extremely useful for media with small and medium velocity gradients where it enables easy and fast implementation of the velocity field while preserving the basic simplicity of a static formulation of the radiative transfer equation.
|Figure 5: The scheme for the expression of the jump of the intensity at the boundary of the cell (see Eq. (34)). The thick square denotes the current cell.|
We again perform a scalar product with the functions wl, i.e.
we multiply the Eq. (32) by functions wl and
integrate it over the volume of the cell,
The solution is performed from deep layers to the stellar surface
to obtain the emergent intensity and in the opposite direction
for the incoming intensity.
The solution starts at the first (either uppermost or lowest) depth point
and the intensity for all frequencies is determined.
After the intensity is known for all frequencies,
we continue with the next
Lower boundary conditions can be obtained,
e.g., from the diffusion
approximation (21) and for the upper boundary condition we
can assume, for example, no incoming radiation.
The boundary conditions for the frequency follow directly from this solution
of transfer equation.
In the first cell
the change of opacity and emissivity caused
by the Doppler effect is very small, so it is sufficient to solve
the static equation there.
The intensity in a grid point
as in the previous case (20),
In this section we describe the solution of the radiative transfer
equation with the inclusion not only of the Doppler effect, but now also
We solve Eq. (3),
We introduce grids
For each cell we choose following space functions as base functions
where , , and . We express the intensity
As in the previous section we use following notation to enable easier description
and the expression for is similar to the case without aberration
After inserting this expression for the intensity into Eq. (38),
multiplying it by functions wl, and integrating
it over the volume of the cell we obtain
We obtained a system of eight linear algebraic equations for eight unknowns Ik that must be supplied with boundary conditions, e.g. the diffusion approximation (21) at the lower boundary and no incident radiation at the upper boundary. The solution is performed in the same manner as in the simpler case without aberration. First, the equation is solved for all frequencies, then for all angles, and finally, for all depth points. The resulting intensity in a grid point we obtain similarly as in the previous cases (20), (37),
|Figure 7: Comparison of the DFE methods of solution of the radiative transfer equation (described in Sects. 3.1 and 4) for the case of a zero velocity gradient. Left panel: the dependence of the Feautrier variable u on the optical depth for a continuum frequency. The solid line (denoted as static) indicates the static solution (2), the dashed line ( LLT) denotes solution using the local Lorentz transformation, the fine dashed line ( with-aberration) means solution including aberration, Eq. (38), and the dotted line ( without-aberration) means the solution without aberration, Eq. (26). Note that all curves coincide. Right panel: the dependence of the intensity in the H line on frequency. Types of lines have the same meaning as in the left panel.|
Test calculations for the simplest case of the static radiative transfer equation have been described in the Sect. 3.1. Here we would like to discuss the choice of an independent variable for the solution of this static equation.
We found that it is better to solve this equation with geometrical depth as an independent variable. If we choose the optical depth as an independent variable, we encounter problems, since the models converge to a wrong solution. The reason lies in the numerical calculation of the system of algebraic equation, which we solve using the LU decomposition. Although the geometrical depth does not change too much throughout the atmosphere, the optical depth changes by several orders of magnitude and its use as an independent variable results in numerical inaccuracy.
We should stress that contrary to what we have said here, Castor et al. (1992) used the optical depth as an independent variable. However, they report a reduced stability for the Galerkin prescription of the mass matrix. They overcome the stability problems with the help of replacing the mass matrix of the system of algebraic Eqs. (17) by a unit matrix. By this step one sets the term AuK-f in the Eq. (7) to be orthonormal with respect to the chosen set of the base functions , which is not generally true. We would like to avoid such a step since in doing it we may risk that the method will converge to a wrong solution, because it will be stable and we will not be able to identify errors.
In order to perform a basic test of our methods in moving media, we compared the solution of the static radiative transfer equation with the solution of the transfer equation in a moving atmosphere for the special case of a zero velocity gradient. All calculations have been done using the DFE method described in Sects. 3.1 and 4. We performed a test solution for a set frequencies from the H line assuming the Doppler profile. Opacities have been calculated for state parameters (electron density) and T (temperature) from an LTE model atmosphere calculated with a code of Kubát (2002). This model atmosphere is presented in the Table 1 (only in electronic form). Results of these calculations are shown in Fig. 7. In the left panel there is a dependence of the Feautrier variable u on an optical depth in four cases - static equation of radiative transfer, equation written in local frequency, equation including aberration, and calculation using local Lorentz transformation. As we can see, all curves coincide, which means that the basic algorithm of the solution of the transfer equation is correct. This result is supported also by a profile of the H line emerging from this model atmosphere.
In Fig. 8 we plotted H line profiles for a nonzero velocity gradient. The velocity field is described by the beta law (25) with the terminal velocity and the velocity in the photosphere . The H line profile is computed for the same condition as in the previous case. Three situations described in preceding sections are shown in Fig. 8 - the local Lorentz transformation, discontinuous finite element method with aberration, and without it.
The profiles shown are in a comoving frame.
For comparison, we plotted also the static solution.
We can see small deviation between the LLT method and methods using the
equation of the radiative transfer written in comoving frame.
The reason for it is in the assumption of a constant velocity in the
We can influence this deviation by
changing the density of the geometrical depth grid,
which is discussed later.
On the other hand, results obtained with the inclusion of aberration
and without it are in a good agreement.
|Figure 8: The dependence of the intensity on frequency in the case of a stellar wind of B stars. The solid line indicates the static solution, Eq. (2), the dashed line the solution using the local Lorentz transformation, the fine dashed line the solution including aberration, Eq. (38), and the dotted line the solution without aberration, Eq. (26).|
In order to see the maximum possible effect of aberration, we solved the transfer equation for a small value of . For this case, we obtained a small difference between the solution including aberration and without it (cf. Fig. 9) as expected. Bearing in mind that no effect of aberration was found in the previous section, we may state that in the region of validity of Eqs. (38) and (26), i.e., in non-relativistic velocity fields, we can safely neglect aberration and use the simpler Eq. (26) in accordance with previous results of Mihalas & Hummer (1976b).
The grid tests show a linear dependence of the computing time
on the number of depth points
(see Fig. 10, left panel).
the right panel
plotted the dependence of the accuracy on the
number of depth grid points.
As we can see, the optimum number of depth points is about one hundred.
For this choice the solution is accurate enough and the computing time
is not too high.
For the velocity field, introduced before, the dependence of computing
time on the number of grid points is linear, too
(see Fig. 11, left panel).
It is necessary to take enough frequency grid points in the
case, because we solve the differential equation with respect to the
now (see Fig. 11, right panel).
We performed the calculations with the same number of grid points in the
geometrical grid and in the frequency grid in the line (D=N=100).
|Figure 9: The dependence of the intensity on frequency for a large angle ( ). The solid line indicates the static solution, the dashed line the solution without aberration, and the dotted line the solution including aberration.|
|Figure 10: Tests of the depth grid. Left panel: the dependence of the computing time on the number of the depth grid points D. Linear fits of these dependences are also plotted. Right panel: the dependence of the accuracy of the solution on the number of the depth grid points D. Here, , where x= 50, 100, 200 and 600.|
|Figure 11: Tests of the frequency grid. Left panel: the dependence of the computing time on a number of frequency grid points within a line. Linear fits of these dependences are also plotted. Right panel: the comparison of the dependence of the line profile on the number of frequency grid points.|
In this paper we studied an application of the discontinuous finite element method for the solution of the radiative transfer equation both for static and differentially moving media.
First, we solved the static equation. We confirmed the fact that the discontinuous finite element method is comparable to the standard Feautrier method in the static case (Fig. 2).
Using the DFE method for the static case we developed a new method for the solution of the radiative transfer equation in moving media, which uses the Lorentz transformation of intensity and frequency between subsequent cells. Since this method uses local inertial frames we call it the local Lorentz transformation method. This method takes advantage of the possibility of a discontinuity at the cell boundary, which in this case comes from the Lorentz transformation.
We solved the transfer equation using the DFE method in the comoving frame. First, we neglected aberration and then we solved it more exactly including both the Doppler shift and aberration. Results of both methods are in good agreement, see Figs. 7 and 8.
In the region of validity of these equations (in non-relativistic velocity fields) it is not necessary to account for aberration, as we expected. This enables faster calculations, since one has to solve the set of only four algebraic equations instead of eight. On the other hand, if we aim at a more precise study of the radiation field, the possibility of taking aberration into account may become very useful.
In addition, we found that it is better to solve these equations with geometrical depth as an independent variable, since choosing the optical depth as an independent variable causes stability problems.
Another important point is a sufficiently fine frequency grid to enable good line profile resolution. If we do not resolve the line sufficiently, we obtain wrong line profiles.
The method of discontinuous finite elements may be applied to the solution of rather complex radiative transfer problems both in static and moving media. We tested it on a problem of expanding atmospheres of hot B stars. However, using this method in other types of moving media is possible as well. We plan to employ our results in future studies of Be stars.
The authors would like to thank the referee for valuable comments on the manuscript. This research has made use of NASA's Astrophysics Data System. This work was supported by a grant GA CR 205/01/0656 and by projects K2043105 and AV 0Z1 003909.