10.4: Nongray Radiative Transfer

While the elimination of the assumption of a gray opacity removes the easy incorporation of radiative equilibrium into the solution of the equation of radiative transfer, most methods described in Section 10.2 can be used for the nongray case. In spite of the diversity of methods available to the researcher for the solution of radiative transfer problems (there are more than are described here), most practical approaches can be divided into two categories: the solution of the integral equation for the source function and methods based on the solution of the differential equations for the radiation field. The solution of the integral equation for the source function is highly efficient, since no more information is generated than is necessary for the solution of the problem, and has also proved effective in dealing with problems of polarization, where complex redistribution functions are required (see Chapter 16). The differential equation approach is perhaps more widely used because a highly efficient algorithm has been developed which enables the investigator to utilize existing and proven mathematical packages for much of the numerical work. In addition, the differential equation approach has proved effective where geometries other than plane-parallel ones are required, and lends itself naturally to the incorporation of time-dependent and hydrodynamic terms when they may be needed. Of the myriads of specific applications, we will be concerned with only two.

a. Solutions of the Nongray Integral Equation for the Source Function

We derived the integral equation for the nongray source function in Section 10.1 [equation \ref{10.1.14}]. The approach we take is basically that described for the solution of the Schwarzschild-Milne equations in Section 10.2. Replacing the integral in equation \ref{10.1.14} with a suitable quadrature scheme, after removing the singularity of the first exponential integral as described in equation \ref{10.2.8}, we get \[\begin{aligned}
S_\nu\left(\tau_\nu\right)= & \epsilon_\nu B\left(\tau_\nu\right)+\frac{1}{2}\left(1-\epsilon_\nu\right)\left\{\sum_{j=1}^n\left[S_\nu\left(t_j\right)-S_\nu\left(\tau_\nu\right)\right] E_1\left|t_j-\tau_\nu\right| W_j\right. \\
& \left.+S_\nu\left(\tau_\nu\right)\left[2-E_2\left(\tau_\nu\right)\right]\right\}
\end{aligned}\label{10.3.1}\]

This functional equation, evaluated at the points of the quadrature, yields a set of linear algebraic equations for the source function at the quadrature points. These, in turn, can be put into standard form so that \[\begin{gathered}
\sum_{k=1}^n S_\nu\left(t_k\right)\left\{\frac{1}{2} \sum_{j=1}^n\left[1-\epsilon\left(t_i\right)\right]\left[\left(\delta_{k j}-\delta_{k i}\right) E_1\left|t_i-t_k\right| W_i+\delta_{i k} E_2\left(t_i\right)\right]+\epsilon\left(t_i\right) \delta_{i k}\right\} \\
=\epsilon_\nu\left(t_i\right) B_\nu\left(t_i\right) \quad i=1, \ldots, n
\end{gathered}\label{10.3.2}\]

These equations are strongly diagonal since the dominant contribution to the source function is always the local one. That contribution is measured by the last term in equation \ref{10.3.1}, and it represents the addition made to the equation to compensate for the removal of the local contribution within the integral. The strongly diagonal nature of the equations ensures that the solution is numerically stable. Indeed, when \(S_\nu=\mathrm{B}_\nu\) and \(\varepsilon_\nu=1\), the equations are formally diagonal. Thus, in practice they may be solved rapidly by means of the Gauss-Seidel iteration with \(S_\nu(t_\mathrm{i})=\mathrm{B}_\nu(t_\mathrm{i})\) as the initial guess. We remarked earlier that some care should be taken in choosing the quadrature scheme. It is a good practice to split the integral into two parts, with the first ranging from 0 to 1 and the second from 1 to 4. A 10-point Gauss-Legendre quadrature provides sufficient accuracy for the rapid change of the source function near the surface, while a 4-point Gauss-Laguerre quadrature scheme is adequate for the second as the source function approaches linearity.

A slightly different approach is taken by the Harvard group² in the widely used atmosphere program called ATLAS. They also solve the integral equation for the source function, but they deal with the singularity of the exponential integral in a somewhat different fashion. Instead of formally removing the singularity, they approximate the source function with cubic splines over a small interval. With an analytic form for the source function, it is possible to evaluate the integral, resulting in a multiplicative weight for the coefficients of the splines. This results in a series of weights which are numerically very similar to those present in equation \ref{10.3.2}. Again, a set of linear algebraic equations is produced for the source function at a discrete set of optical depths. The results of the two methods are nearly identical, with the gaussian quadrature scheme being somewhat more efficient.

b. Differential Equation Approach: The Feautrier Method

This method replaces the differential equations of radiative transfer with a set of finite difference equations for parameters related to the specific intensity at a discrete set of values for the angular variable \(\mathrm{m_i}\). However, the choice of values of mi is irrelevant to understanding the method, so we leave that choice arbitrary for the moment. Instead of solving the equation of transfer for the specific intensity, we write equations of transfer for combinations of inward- and outward-directed streams.

Feautrier Equations Consider the variables \[u \equiv \frac{1}{2}\left[I\left(+\mu, \tau_\nu\right)+I\left(-\mu, \tau_\nu\right)\right] \quad v \equiv \frac{1}{2}\left[I\left(+\mu, \tau_\nu\right)-I\left(-\mu, \tau_\nu\right)\right]\label{10.3.3}\]

Here we have paired the outward directed stream \(\mathrm{I(+\mu,\tau)}\) with its inward -directed counterpart \(\mathrm{I(-\mu,\tau)}\) into quantities that resemble a "mean" intensity u and a "flux" v. One of the benefits of the linearity of the equation of transfer is that we can add or subtract such equations and still get a linear equation. Thus, by adding an equation for a \(+\mu\) stream to one for a \(-\mu\) stream we get \[\mu \frac{d v}{d \tau_\nu}=u-S \quad \mu>0\label{10.3.4}\]

Similarly, by subtracting one from the other, we get \[\mu \frac{d u}{d \tau_\nu}=v \quad \mu>0\label{10.3.5}\]

Using this result to eliminate v from equation \ref{10.3.4}, we have \[\mu^2 \frac{d^2 u}{d \tau_\nu^2}=u-S \quad \mu>0\label{10.3.6}\]

This is a second-order linear differential equation, so we will need two constraints or constants of integration. At the surface \(\mathrm{I(-\mu,0)}=0\), so \(\mathrm{v(0)=u(0)}\), and from equation \ref{10.3.5} we have \[\left.\mu \frac{d u}{d \tau_\nu}\right|_{\tau_\nu=0}=u(0)\label{10.3.7}\]

The other constraint on the differential equation comes from invoking the diffusion approximation at large depths. Under this assumption \[I_\nu \simeq J_\nu \simeq B_\nu \simeq S_\nu\label{10.3.8}\]

We may now use the equation of transfer itself to generate a perturbation expression for \(\mathrm{I}(\mu,\tau_\nu)\) at large depths: \[\begin{aligned}
I_\nu\left(\mu, \tau_\nu\right) & =\mu \frac{d I_\nu\left(\mu, \tau_\nu\right)}{d \tau_\nu}+S_\nu\left(\tau_\nu\right) \\
& \approx \mu \frac{d B_\nu\left(\tau_\nu\right)}{d \tau_\nu}+B_\nu\left(\tau_\nu\right) \quad \tau_\nu \gg 1
\end{aligned}\label{10.3.9}\]

Substituting the left-hand side into the definition for v we get \[v=\mu \frac{d B_\nu\left(\tau_\nu\right)}{d \tau_\nu}=\mu \frac{d u}{d \tau_\nu} \quad \tau \gg 1\label{10.3.10}\]

Equations \ref{10.3.7} and \ref{10.3.10} are the two constraints needed to specify the solution. Now consider the finite difference approximations required to solve equation \ref{10.3.6} subject to these constraints.

Solution of the Feautrier Equations We saw earlier, in Chapter 4, how the method of solution used to solve the Schwarzschild equations of stellar structure was supplanted by the Henyey method utilizing finite differences. Many of the reasons that lead to the superiority of the Henyey method are applicable to the Feautrier method for solution of the equations of radiative transfer. It is for that reason that we describe the numerical method in some detail.

First, we must pick a set of \({\tau_\mathrm{k}}’\mathrm{s}\) for which we desire the solution. We must be certain that the largest \(\tau_N\) is deep enough in the atmosphere to ensure that the assumptions resulting in the boundary condition given in equation \ref{10.3.9} are met. In addition, it is useful if the density of points near the surface is large enough that the solution will be accurately described. This is particularly important when we are dealing with the transport of radiation within a spectral line. Now we define the following finite difference operators: \[\begin{aligned}
\Delta f_{k+1 / 2} & \equiv f\left(\tau_{k+1}\right)-f\left(\tau_k\right) \\
\Delta \tau_{k+1 / 2} & \equiv \tau_{k+1}-\tau_k \\
\Delta \tau_k & \equiv \frac{1}{2}\left(\Delta \tau_{k+1 / 2}+\Delta \tau_{k-1 / 2}\right)
\end{aligned}\label{10.3.11}\]

The subscript k+½ simply means that this is an estimate of the parameter appropriate for the value of τ midway between k and k + 1. Unlike the Henyey scheme, where this information was obtained from an earlier model structure, the Feautrier method obtains the information by linear interpolation from the existing solution. Now we replace the derivatives with the following finite difference operators: \[\begin{aligned}
\left.\frac{d f(\tau)}{d \tau}\right|_{\tau=\tau_k} & \simeq \frac{\Delta f_{k+1 / 2}}{\Delta \tau_{k+1 / 2}} \\
\left.\frac{d^2 f(\tau)}{d \tau^2}\right|_{\tau=\tau_k} & \simeq \frac{\Delta f_{k+1 / 2} / \Delta \tau_{k+1 / 2}-\Delta f_{k-1 / 2} / \Delta \tau_{k-1 / 2}}{\Delta \tau_k}
\end{aligned}\label{10.3.12}\]

The second derivative in equation \ref{10.3.6} can now be replaced by these operators operating on \(\mathrm{u}(\mu)\) to yield the following linear algebraic equations for \(\mathrm{u}(\mu)\) at the chosen optical depth points \(\tau_\mathrm{k}\): \[\begin{aligned}
&\mu^2 \frac{u_{k-1}(\mu)}{\Delta \tau_k \Delta \tau_{k-1 / 2}}-\frac{\mu^2}{\Delta \tau_k}\left(\frac{1}{\Delta \tau_{k-1 / 2}}+\frac{1}{\Delta \tau_{k+1 / 2}}\right) u_k(\mu)+\frac{\mu^2 u_{k+1}(\mu)}{\Delta \tau_k \Delta \tau_{k+1 / 2}} \\
&=u_k(\mu)-S_k
\end{aligned}\label{10.3.13}\]

Now it is time to pick those values of µ for which we desire the solution. Let u be considered a vector whose elements are the values of u at the particular values of \(\mu_\mathrm{i}\), so that \[\vec{u}=\left[u\left(\mu_1\right), u\left(\mu_2\right), u\left(\mu_3\right), \ldots, u\left(\mu_n\right)\right]\label{10.3.14}\]

The linear equations \ref{10.3.13} can now be written as a system of matrix-vector equations of the form \[\mathbf{A}_k \cdot \vec{u}_{k-1}-\mathbf{B}_k \cdot \vec{u}_k+\mathbf{C}_k \cdot \vec{u}_{k+1}=-\vec{S}_k\label{10.3.15}\]

The elements of matrices A, B, and C involve only the values of \(\mu_\mathrm{i}\) and \(\tau_\mathrm{k}\) that were chosen to describe the solution. The constraints given by equations \ref{10.3.7} and \ref{10.3.9} require that \[\mathbf{A}_\boldsymbol{1}=\mathbf{0} \quad \mathbf{C}_\boldsymbol{N}=\mathbf{0}\label{10.3.16}\]

Thus we have set the conditions required to solve the equations for \(\mathrm{u_k(\mu_i)}\) from which the specific intensity can be recovered and all the moments that depend on it. Equations \ref{10.3.16} happen to be tridiagonal, which ensures that they can be solved efficiently and accurately. We have glossed over the source function \(S_\mathrm{k}\) in our discussion by assuming that it is known everywhere and depends only on \(\tau_\mathrm{k}\). However, the property of the source function that caused so many problems for earlier methods (and, indeed, resulted in the integral equations in the first place) is that the source function usually depends on the intensity itself. However, for scattering, the source function does depend on the intensity in a linear manner. Therefore, it is possible to represent the source function in terms of the unknowns \(\mathrm{u_k}\) and \(\mathrm{v_k}\) and include them in the equations, still preserving their tridiagonal form.

There is one caveat to this. The Feautrier method imposes a certain symmetry on the solution to the radiative transfer problem by combining inward- and outward-directed streams. If the redistribution function does not share this symmetry, it will not be possible to represent the scattering in terms of the functions \(\mathrm{u(\mu)}\) and \(\mathrm{v(\mu)}\). Thus, for some problems involving anisotropic scattering, the Feautrier method may not be applicable. In addition, when the redistribution function involves redistribution in frequency, the optical depth points must be chosen so that the deepest point will satisfy the assumptions required for the approximation given in equation \ref{10.3.9} for all frequencies. If this is not done, errors incurred at those frequencies for which the assumptions fail can propagate in an insidious manner throughout the entire solution.

The Feautrier method does not suffer from the exponential instabilities described for the discrete ordinate method, because it invokes the diffusion approximation at large depths (specifically the inner boundary). The diffusion approximation basically contains the information that the radiation field has been randomized in direction and thereby stabilizes the solution in the same manner as it stabilizes the Eddington solution. As we see in Chapter 11, knowledge of the mean intensity, the radiative flux, and occasionally the radiation pressure is usually sufficient to calculate the structure of the atmosphere. The Feautrier method finds more information than that and therefore is not as efficient as it might be. However, the numerical methods for solving the resulting linear equations are so fast that the overall efficiency of the method is quite good, and it provides an excellent method of solution for most problems of radiative transfer in stellar atmospheres. Remember that, like any numerical method, the Feautrier method should be used with great care and only on those problems for which it is suited.