10.5: Radiative Transport in a Spherical Atmosphere

Any discussion of the solution of radiative transfer problems would be incomplete without some mention of the problem introduced by a departure from the simplifying assumption of plane-parallel geometry. In addition, there are stars for which the plane-parallel approximation is inappropriate, and we would like to model these stars as well as the main sequence stars for which the plane-parallel approximation is generally adequate. The density in the outer regions of red supergiants is so low that the atmosphere will occupy the outer 30 percent to 40 percent of what we would like to call the radius of the star. Here, the plane-parallel assumption is clearly inappropriate for describing the star. We must include the curvature of the star in any description of its atmosphere. In doing so, we will require a parameter that was removed by the plane-parallel assumption - the stellar radius. This parameter can be operationally defined as the distance from the center to some point where the radial optical depth to the surface is some specified number (say unity). In doing so, we must remember that the radius may now become a wavelength-dependent number and so some mean value from which the majority of the energy escapes to the surrounding space may be appropriate for describing the star when a single value for the radius is required. However, for the calculation of the stellar interior, we need to know only the surface structure at a given distance from the center in order to specify the interior structure. Whether the distance corresponds to our idea of a stellar radius is irrelevant. In addition, we assume that the star is spherically symmetric.

a. Equation of Radiative Transport in Spherical Coordinates

In Chapter 9 we developed a very general equation of radiative transfer which was coordinate-independent [equation \ref{9.2.11}]. Writing the time-independent form for which the gravity gradient does not significantly affect the photon energy, we get \[\hat{n} \cdot \nabla I_\nu=\rho\left(\kappa_\nu+\sigma_\nu\right)\left(S_\nu-I_\nu\right)\label{10.4.1}\]

Writing the ∇ operator in spherical coordinates and making the usual definition for m (see Figure 10.2), we get \[\mu \frac{\partial I_\nu(r, \mu)}{\partial r}+\frac{1-\mu^2}{r} \frac{\partial I_\nu(r, \mu)}{\partial \mu}=\rho\left(\kappa_\nu+\sigma_\nu\right)\left[S_\nu(r, \mu)-I_\nu(r, \mu)\right]\label{10.4.2}\]

where we take the source function to be that of a nongray atmosphere with coherent isotropic scattering, so that \[S_\nu=\epsilon_\nu B_\nu[T(r)]+\left(1-\epsilon_\nu\right) J_\nu(r)\label{10.4.3}\]

Our approach to the solution of the equation of transfer will be to obtain and solve some equations for the important moments of the radiation field.

Radiative Equilibrium and Moments of the Radiation Field For a steady state atmosphere, our condition for radiative equilibrium [equation \ref{9.4.4}] becomes \[\nabla \cdot \int_0^{\infty} \vec{F}_\nu d\nu=4 \int_0^{\infty} \kappa_\nu \rho\left(B_\nu-J_\nu\right) d\nu=0\label{10.4.4}\]

However, in spherical coordinates, the divergence of the total flux yields the same condition that we obtained for stellar interiors [equation \ref{4.2.1}]: \[\pi F \equiv \pi \int_0^{\infty} F_\nu d\nu=\frac{L}{4 \pi r^2}\label{10.4.5}\]

Now the condition of radiative equilibrium is obtained from the zeroth moment of the equation of transfer [equation \ref{9.4.3}], while the first moment of the equation of transfer [equation \ref{9.4.6}] yields an expression for the radiation pressure tensor. For an atmosphere with no time-dependent processes, these moment equations become \[\nabla \cdot \vec{F}_\nu=4 \kappa_\nu \rho\left(B_\nu-J_\nu\right) \quad \nabla \cdot \mathbf{K}_\nu=-\frac{\rho\left(\kappa_\nu+\sigma_\nu\right) \vec{F}_\nu}{4}\label{10.4.6}\]

Noting that there is no net flow of radiation in either the \(\theta\) or \(\phi\) coordinates for a spherically symmetric atmosphere, we see that the divergence of the flux in spherical coordinates becomes \[\frac{\partial\left(r^2 F_\nu\right)}{\partial r}=4 r^2 \kappa_\nu \rho\left(B_\nu-J_\nu\right)\label{10.4.7}\]

If we make the assumption that the radiation field is nearly isotropic, then \(\nabla \cdot \mathbf{K}_\nu\) becomes \(\nabla \mathrm{K}_\nu\) where \(\mathrm{K}_\nu\) is the scalar moment that we have identified with the radiation pressure [see equations \ref{9.3.14} through \ref{9.3.16}]. Perhaps the easiest way to find the representation of equation \ref{10.4.6} in spherical coordinates is to multiply equation \ref{10.4.2} by m and integrate over all m. This yields the second of the required moment equations, \[\frac{\partial K_\nu}{\partial r}+\frac{3 K_\nu-J_\nu}{r}=-\frac{\rho\left(\kappa_\nu+\sigma_\nu\right) F_\nu}{4}\label{10.4.8}\]

Figure 10.2 shows the geometry assumed for the Spherical Equations of radiative transfer. The angle \(\theta\) for which \(\mu=\cos\theta\) is defined with respect to the radius vector. Unlike the plane-parallel approximation the depth variable is the radius and increases outward.

Closing the Moment Equations and the Eddington Factor In Chapter 9 we observed [equation \ref{9.4.8}] that under conditions of near isotropy \(\mathrm{K}_\nu=\mathrm{J}_\nu/3\). This was the moment approximation needed to close the moment equations, and it is known as the diffusion approximation. However, such conditions do not prevail throughout the atmosphere, so it is common to assume that the two moments can be related by a scale factor, which has come to be known as the Eddington factor, defined as \[f_\nu(r)=\frac{K_\nu(r)}{J_\nu(r)}\label{10.4.9}\]

We can replace the radiation pressure by the Eddington factor and obtain \[\frac{\partial\left(f_\nu J_\nu\right)}{\partial r}+\frac{\left(3 f_\nu-1\right) J_\nu}{r}=\frac{-\rho\left(\kappa_\nu+\sigma_\nu\right) F_\nu}{4}\label{10.4.10}\]

for the second moment equation.

Equation \ref{10.4.10} combined with equation \ref{10.4.7} form a complete system for \(\mathrm{F}_\nu\) and \(\mathrm{J}_\nu\) subject to the appropriate boundary conditions. Of course, we have not fundamentally changed the problem since the Eddington factor is unknown and presumably a function of depth. It must be found so that any atmosphere produced is self-consistent under the constraint of radiative equilibrium. The Eddington factor basically measures the isotropy of the radiation field, since for isotropic radiation it is 1/3. Imagine a radiation field entirely directed along \(\mu=\pm1\). For such a field \(f_\nu=1\), while for a radiation field confined to a plane that is normal to this direction, \(f_\nu=0\). If we consider the normal radiation field emerging from a star, the temperature gradient normally produces limb-darkening, implying that the radiation field near the surface becomes more strongly directed along the normal to the atmosphere. Thus, we should expect the Eddington factor to increase as the surface approaches. This effect should be enhanced for stars with large spherical atmospheres. Thus, for normal stellar atmospheres \[\frac{1}{3} \leq f_\nu(r)<1\label{10.4.11}\]

b. An Approach to Solution of the Spherical Radiative Transfer Problem

Sphericality Factor This factor is introduced purely for mathematical convenience and as such has no major physical importance. However, it does tend to make the spherical moment equations resemble their plane-parallel counterparts. We define \[\ln \left(r^2 q_\nu\right) \equiv \int_{r_c}^r \frac{3 f_\nu-1}{x f_\nu} d x+\ln r_c^2\label{10.4.12}\]

so that \[\frac{-4}{\rho\left(\kappa_\nu+\sigma_\nu\right) q_\nu} \frac{\partial\left(r^2 q_\nu f_\nu J_\nu\right)}{\partial r}=r^2 F_\nu\label{10.4.13}\]

The parameter \(\mathrm{r_c}\) is the deepest radius for which the problem is to be solved. Given \(\mathrm{F}_\nu\), we can find the sphericality factor \(q_\nu\) by numerically integrating equation \ref{10.4.12}. Using this definition of \(q_\nu\), we may rewrite the second moment equation \ref{10.4.10}, as \[\frac{-4}{\rho\left(\kappa_\nu+\sigma_\nu\right) q_\nu} \frac{\partial\left(r^2 q_\nu f_\nu J_\nu\right)}{\partial r}=r^2 F_\nu\label{10.4.14}\]

This form is suitable for combining with the first moment equation \ref{10.4.7}, to eliminate \(\mathrm{F}_\nu\) and get \[\frac{\partial^2\left(r^2 q_\nu f_\nu J_\nu\right)}{\partial \tau_\nu^2}=\frac{r^2 \epsilon_\nu\left(J_\nu-B_\nu\right)}{q_\nu}\label{10.4.15}\]

where \[\partial \tau_\nu=-q_\nu \rho\left(\kappa_\nu+\sigma_\nu\right) \partial r\label{10.4.16}\]

and \(\varepsilon_\nu\) has the same meaning as before [see equation \ref{10.1.8}]. We have now generated a second order differential equation for \(\mathrm{J}_\nu\) that is similar to the one obtained for the Feautrier method, and we solve it in a similar manner.

Boundary Conditions The boundary conditions are determined in much the same manner as for the Feautrier method. For the lower boundary we make the same assumptions of isotropy as were made for equation \ref{10.3.9}. Indeed, we multiply equation \ref{10.3.9} by \(\mu\) and integrate over all \(\mu\), to get \[F_\nu \simeq 2 \int_{-1}^{+1} \mu^2 \frac{d B_\nu}{d \tau_\nu} d \mu+2 \int_{-1}^{+1} \mu B_\nu d \mu=\frac{4}{3} \frac{d B_\nu}{d \tau_\nu}\label{10.4.17}\]

This and equation \ref{10.4.14} allow us to specify the derivative of \(\mathrm{J}_\nu\) at the lower boundary as \[\left.\frac{\partial\left(r^2 q_\nu f_\nu J_\nu\right)}{\partial \tau_\nu}\right|_{r=r_c}=\left.\frac{r_c^2}{3} \frac{d B_\nu}{d \tau_\nu}\right|_{r=r_c}\label{10.4.18}\]

Again \(\mathrm{r_c}\) is the deepest point for which the solution is desired. Equation \ref{10.4.14} also sets the upper boundary condition at R as \[\left.4 \frac{\partial\left(r^2 q_\nu f_\nu J_\nu\right)}{\partial \tau_\nu}\right|_{r=R}=R^2 F_\nu=R^2 \frac{\int_0^1 \mu I(R, \mu) d \mu}{\int_0^1 I(R, \mu) d \mu} J_\nu(R)\label{10.4.19}\]

so that we again have a two-point boundary-value problem and a second-order differential equation for \(\mathrm{J}_\nu\) which we can solve by the same finite difference techniques that were used for the Feautrier method [see equations \ref{10.3.11} through \ref{10.3.16}].

The problem can now be solved, assuming we know the behavior of the Eddington factor with depth in the atmosphere. Unfortunately, to find this, we must know the angular distribution of the radiation field at all depths. Normally, we could appeal to the classical solution, for knowledge of \(\mathrm{J}_\nu\) would provide all the information needed to calculate the source function. But the classical solution was appropriate for only the plane-parallel approximation. To find the analog for spherical coordinates, we have to use the symmetry of a spherical atmosphere and perform still another coordinate transformation.

Impact Space and Formal Solution for the Spherical Equation of Radiative Transfer Consider a coordinate frame attached to the star so that the z axis points in the direction of the observer and passes through the center of the star (see Figure 10.3). Coordinates z and p designate all places within the star with p playing the role of an impact parameter for photons directed toward the observer parallel to z. The entire solution set \(\mathrm{I(\mu,r)}\) can be represented by the radiation streams \(\mathrm{I}_+(p,\mathrm{z})\) and \(\mathrm{I}-(p,\mathrm{z})\) by moving along surfaces of constant r.

Figure 10.3 describes 'impact space' for spherical transport. The z-axis points at the observer, while the p -coordinate is perpendicular to z and plays the role of an impact parameter for the photons directed toward the observer. The angle α denotes the angle between a line parallel to z, directed toward the observer, and a radius vector.

Thus any solution that gives us a complete representation of the specific intensity in the p-z plane will give a complete description of the radiation field. We can immediately write the equation of transfer for the special beams directed toward or away from the observer as \[\pm \frac{\partial I_{ \pm}(p, z)}{\partial z}=\rho\left(\kappa_\nu+\sigma_\nu\right)\left[S(p, z)-I_{ \pm}(p, z)\right]\label{10.4.20}\]

where the coordinate transformation from p-z coordinates to \(\mu-\mathrm{r}\) coordinates is \[r=\left(p^2+z^2\right)^{1 / 2} \quad|\mu|=\frac{z}{\left(p^2+z^2\right)^{1 / 2}}\label{10.4.21}\]

For simplicity we denote \[k_\nu=\rho\left(\kappa_\nu+\sigma_\nu\right)\label{10.4.22}\]

Equation \ref{10.4.20} is a linear first order equation that has a classical solution of \[\begin{aligned}
& I_{-}(p, z)=\int_z^{\left(R^2-p^2\right)^{1 / 2}} k_\nu(\xi) S_\nu(\xi) e^{-\tau(p, \zeta, z)} d \zeta \\
& I_{+}(p, z)=\int_0^z k_\nu(\xi) S_\nu(\xi) e^{-\tau(p, z, \zeta)} d \zeta+I_{-}(p, 0) e^{-\tau(p, z, 0)} \\
& \xi^2 \equiv p^2+\zeta^2 \quad \tau(p, a, b)=\int_a^b k(\xi) d \zeta
\end{aligned}\label{10.4.23}\]

While this is a complicated expression, it can be evaluated numerically as long as one has a representation of the source function. Thus, it is now possible to solve for the entire radiation field and recalculate the variable Eddington factor \(f_\nu\). Equation \ref{10.4.15} is then solved again for a new value of \(\mathrm{J}_\nu\) and hence \(S_\nu\). The entire procedure is repeated until a self-consistent solution is found. Rather than carry out the admittedly messy numerical integration, Mihalas³ describes a Feautrier-like method to calculate the intensities directly.

A method proposed by Schmid-Burgk⁴ assumes that the source function can be locally represented by a polynomial in the optical depth. This analytic function is then substituted into the formal solution in impact space so that the radiation field can be represented in terms of the undetermined coefficients of the source function's approximating polynomials. The moments of the radiation field can then be generated which depend only on these same coefficients. Thus, if one starts with an initial atmospheric structure and a guess for the source function, one can fit that source function to the local polynomial and thereby determine the approximating coefficients. These, in turn, can be used to generate the moments of the radiation field upon which an improved\nuersion of the source function rests. An excellent initial guess for the source function is \(S_\nu=\mathrm{B}_\nu\), and unless scattering completely dominates the opacity, the iteration process converges very rapidly.

It is clear that the spherical atmosphere poses significant difficulties over and above those found in the plane-parallel atmosphere. However, there are very few differences that are fundamental in nature. All present methods rely on the global symmetry of spherical stars, and it seems likely that those stars with atmospheres sufficiently extended to require the spherical treatment will also be subject to other forces, such as rotation, that further distort the atmospheres so that even this global symmetry is lost. However, such studies can offer insight into the severity of the effects that we can expect from the geometry.

We have only skimmed the surface of the methods and techniques devised to solve the equation of radiative transfer. The methods discussed merely comprise some of the more popular and successful methods currently in use. We have left to the studious reader the entire area of the "exact approximation" and the H-functions of Chandrasekhar¹ (pp. 105 to 126). No mention has been made of invariant embedding and the voluminous literature written for Linear two-point boundary-value problems. Many of these techniques have proved useful in solving specific radiative transfer problems, and those who would count themselves experts in this area should avail themselves of that literature. There is an entire field of study surrounding the transfer of radiation within spectral lines, some of which will be discussed later, but much of which will not be. This material is important for anyone interested in problems requiring line-transfer solutions. However, the methods presented here suffice for providing the solution to half of the task of constructing a normal stellar atmosphere, and next we turn to the solution of the other half of the problem.