10.3: Gray Atmosphere

For the better part of this century, theoretical astrophysicists have been concerned with the solution to an idealized radiative transfer problem known as the gray atmosphere. Although it is an idealized situation, it has some counterparts in nature. In addition, this problem possesses the virtue that a complete solution can be obtained for the radiation field without recourse to the physical details of the atmosphere. In this regard, the gray atmosphere model is rather like polytropic models for stellar interiors. As was the case for polytropes and stellar interiors, we may expect to gain significant insight into the properties of stellar atmospheres by understanding the solution to the gray atmosphere problem. The additional assumption required to turn our study of radiative transfer into that of a gray atmosphere is simple. Assume that the opacity, whether it is absorption or scattering, is independent of frequency. Thus, any frequency can be treated as any other frequency, as far as the radiative transfer is concerned. This independence of the radiative transfer from frequency has the interesting consequence that the mathematical solution to the equation of transfer for any frequency will be the solution for all frequencies, and thus must be the solution for the sum of all frequencies. Hence, the aspect of the solution that specifies the radiative flux also refers to the total flux, making the condition of radiative equilibrium relatively simple to apply. Since all aspects of the mathematical description are independent of frequency, we drop the subscript n for the balance of this discussion.

Knowing what we do about the physical processes of absorption, it is reasonable to ask if the gray atmosphere is anything more than an interesting mathematical exercise. Certainly bound-bound transitions are anything but gray. However, there are some bound-free transitions that exhibit only weak frequency dependence over substantial regions of the spectrum. If those regions of the spectrum correspond to that part of the spectrum containing most of the radiant flux, then the atmosphere will be very similar to a gray atmosphere. Absorption due to the H-minus ion is relatively frequency-independent throughout the visible part of the spectrum and in some stars is the dominant source of opacity. However, the premier example of a gray opacity source is electron scattering. Thomson scattering by free electrons is frequency-independent by definition, and for stars hotter than about 25,000 K, it is the dominant source of opacity throughout the range of wavelengths encompassing the maximum flow of energy. Thus, the early O and B stars have atmospheres that, to a very high degree, may be regarded as gray.

Since frequency dependence has been removed from the problem, we may write the equation of radiative transfer for a plane-parallel static atmosphere as \[\mu \frac{d I(\tau, \mu)}{d \tau}=I(\tau, \mu)-S(\tau)\label{10.2.1}\]

where, for isotropic coherent scattering, the source function is \[S(\tau)=\frac{\kappa B+\sigma J}{\kappa+\sigma}\label{10.2.2}\]

Now the independence of the opacity on frequency makes the condition of radiative equilibrium given by equation \ref{9.4.4} particularly simple. \[\int_0^{\infty} \kappa_\nu \rho\left(B_\nu-J_\nu\right) d v=0=B-J\label{10.2.3}\]

or simply \[B=J\label{10.2.4}\]

Substitution of this result into equation \ref{10.2.2} yields \[S(\tau)=B[T(\tau)]=J(\tau)\label{10.2.5}\]

The fact that the mean intensity is equal to the Planck function and that either can be taken to be the source function has the interesting result that the solution to the gray atmosphere is independent of the relative roles of scattering and absorption. Thus, the radiation field for a pure absorbing gray atmosphere, where the source function is clearly the Planck function, will be indistinguishable from the radiation field of a pure scattering gray atmosphere. In addition, since there is a general independence on frequency, the spectral energy distribution will be that resulting from a gray atmosphere where the source function is the Planck function.

The gray atmosphere implies that all the development of Chapters 9 and 10 will apply at each frequency. This is indeed the easiest way to obtain equations \ref{10.2.1} through \ref{10.2.4}. But there is much more. The integral equation for the source function [equation \ref{10.1.14}] and that for the moments of the radiation field [equations \ref{10.1.17}] become \[\begin{aligned}
& B(\tau)=\frac{1}{2} \int_0^{\infty} B(t) E_1|t-\tau| d t \quad J(\tau)=\frac{1}{2} \int_0^{\infty} J(t) E_1|t-\tau| d t \\
& H(\tau)=\frac{1}{2} \int_\tau^{\infty} B(t) E_2(t-\tau) d t-\frac{1}{2} \int_0^\tau B(t) E_2(\tau-t) d t \\
& K(\tau)=\frac{1}{2} \int_0^{\infty} B(t) E_3|t-\tau| d t
\end{aligned}\label{10.2.6}\]

Solution of these equations, combined with the classical solution to the equation of transfer, yields a complete description for the radiation field at all depths in the atmosphere. The method of solution for the gray atmosphere equation of transfer is also illustrative of the methods of solution for the more general nongray problem.

a. Solution of Schwarzschild-Milne Equations for the Gray Atmosphere

In general, an accurate solution of these equations must be accomplished numerically because the solution, even for the gray atmosphere, is not analytic everywhere. Particular care must be taken with these equations because the first exponential integral behaves badly as its argument approaches zero. Specifically \[\operatorname{Lim}_{x \rightarrow 0} E_1(x)=-\ln x \rightarrow \infty\label{10.2.7}\]

Thus, the kernel of first two of equations \ref{10.2.6} has a singularity when \(\mathrm{t}=\tau\). However, this singularity is integrated over, and the integral is finite and well behaved. For years this singularity was regarded as an insurmountable barrier, and interest in the solution of the integral equations of radiative transfer languished in favor of more direct methods applicable to the differential equation of transfer itself. However, the singularity of the kernel is not an essential one and may be easily removed. Simply adding and subtracting the solution \(\mathrm{B}(\tau)\) from the right-hand side of the first of equations \ref{10.2.6} yields \[B(\tau)=\frac{1}{2} \int_0^{\infty}[B(t)-B(\tau)] E_1|t-\tau| d t+\frac{1}{2} B(\tau) \int_0^{\infty} E_1|t-\tau| d t\label{10.2.8}\]

The integrand of the first of these integrals is now well-behaved for all values of (t) since \([\mathrm{B(t)-B(\tau)}]\) will go to zero faster than the exponential integral diverges as \(\mathrm{t}\rightarrow\tau\). The only condition placed on the solution is that \(\mathrm{B}(\tau)\) satisfy a Lipschitz condition which is a weaker condition than requiring the solution to be continuous. The second integral is analytic and can be evaluated by using the properties of exponential integrals given in equations \ref{10.1.12}. This yields a slightly different integral equation, but one that has a well behaved integrand: \[B(\tau)=\frac{\int_0^{\infty}[B(t)-B(\tau)] E_1|t-\tau| d t}{E_2(\tau)}\label{10.2.9}\]

A simple way to deal with this type of integral equation is to replace the integral with some standard numerical quadrature formula. While Simpson's rule enjoys a great popularity, a gaussian-type quadrature scheme offers much greater accuracy for the same number of points of evaluation of the integrand. When the integral is so replaced, we obtain \[B(\tau)=\frac{\sum_{i=1}^n\left[B\left(t_i\right)-B(\tau)\right] E_1\left|t_i-\tau\right| W_i}{E_2(\tau)}\label{10.2.10}\]

which is a functional equation for \(\mathrm{B}(\tau)\) in terms of the solution at a discrete set of points \(t_\mathrm{i}\). The quantities \(W_\mathrm{i}\) are just the weights of the quadrature scheme appropriate for the various points \(\mathrm{t_i}\). Evaluating the functional equation for \(\mathrm{B}(\tau)\) with \(\tau\) equal to each value of \(t_\mathrm{j}\), and rearranging terms, we can obtain a system of linear algebraic equations for the solution at the specific points \(t_\mathrm{i}\): \[\sum_{k=1}^n B\left(t_k\right)\left[\sum_{i=1}^n \frac{\left(\delta_{i k}-\delta_{j k}\right) E_1\left|t_i-t_j\right| W_i}{E_2\left(t_j\right)}-\delta_{k j}\right]=0 \quad j=1, \ldots, n\label{10.2.11}\]

The term governed by the summation over i depends only on the type of quadrature scheme chosen, and so the equation \ref{10.2.11} represents n linear homogeneous algebraic equations that have the standard form \[\sum_{k=1}^n B\left(t_k\right) A_{k j}=0 \quad j=1, \ldots, n\label{10.2.12}\]

The fact that these equations are homogeneous points out an observation made earlier. For the gray atmosphere, the radiation field is decoupled from the values of the physical state variables. Thus, the homogeneous equations constitute an eigenvalue problem, and, as we see later, the eigenvalue is the value of the total radiative flux or alternately the effective temperature. One approach to the solution of equations \ref{10.2.12} would be to define a new set of variables \(\mathrm{B}(t_\mathrm{i})/\mathrm{B}(t_1)\) say, and to generate a system of inhomogeneous equations that can then be solved for the ratio of the source function to its value at one of the given points. Once the source function (or its ratio) has been found at the discrete points \(t_{\mathrm{i}}\), the solution can be obtained everywhere by substitution into equation 10.2.10. Since this is a functional equation, the results will have the same level of accuracy as that obtained for the values of \(\mathrm{B}(t_\mathrm{i})\). To achieve a level of accuracy significantly greater than that offered by the Eddington approximation, we will have to use a particularly accurate quadrature formula. Also the exponential nature of the exponential integral implies that the quadrature scheme should be chosen with great care.

b. Solutions for the Gray Atmosphere Utilizing the Eddington Approximation

We have already seen that the diffusion approximation yields moment equations from the equation of transfer given by equation \ref{9.4.11}. For the gray atmosphere, these take the particularly simple form \[\frac{d F}{d \tau}=0 \quad \frac{d J}{d \tau}=\frac{3}{4} F\label{10.2.13}\]

The first is a statement of radiative equilibrium which says that for a gray atmosphere \(\mathrm{F}_\nu\) is constant, and its integrated value can be related to the effective temperature. The second equation is immediately integrable, yielding a constant of integration. Thus, \[F(\tau)=\text { const }=\frac{\sigma T_e^4}{\pi} \quad J(\tau)=\frac{3}{4} F \tau+\text { const }\label{10.2.14}\]

Using the Eddington approximation as given by equation \ref{9.4.13}, we can evaluate the constant and arrive at the dependence of the mean intensity with depth in the atmosphere. \[J(0)=\frac{1}{2} F=\text { const } \quad J(\tau)=\frac{3}{4} F\left(\tau+\frac{2}{3}\right)\label{10.2.15}\]

Remembering that J = S = B for a gray atmosphere in radiative equilibrium, we find that the temperature of the atmosphere should vary as \[\left[\frac{T(\tau)}{T_e}\right]^4=\frac{3}{4}\left(\tau+\frac{2}{3}\right)\label{10.2.16}\]

Thus, we see that at large depths, where we should expect the diffusion approximation to yield accurate results, the source function becomes linear with depth. Also, when \(\tau=2/3\), the local temperature equals the effective temperature. So, in some real sense, we can consider the optical "surface" to be located at (\tau=2/3\). This is the depth from which the typical photon emerges from the atmosphere into the surrounding space. Only at depths less than 2/3 does the source function begin to depart significantly from linearity with depth. Unfortunately, this is the region in which most of the spectral lines that we see in stellar spectra are formed. Thus, we will have to pay special attention to that part of the atmosphere lying above optical depth 2/3.

We may check on the accuracy of the Eddington approximation by seeing how well it reproduces the surface boundary condition that it assumes. Using the definition for the mean intensity, the classical solution for the equation of transfer [equation \ref{10.1.5}], and the fact that the source function is J itself, we obtain \[J(0)=\frac{1}{2} \int_0^1 I(\mu, 0) d \mu=\frac{1}{2} \int_0^1 \int_0^{\infty} J(t) e^{-t / \mu} \frac{d t}{\mu} d \mu=\frac{7 F}{16}\label{10.2.17}\]

So the Eddington approximation fails to be self-consistent by about 1 part in 8 or 12.5 percent in reproducing the surface value for the flux. To improve on this result, we will have to take a rather more complicated approach to the radiative problem.

c. Solution by Discrete Ordinates: Wick-Chandrasekhar Method

The following method for the solution of radiative transfer problems has been extensively developed by Chandrasekhar¹ and we only briefly sketch it and its implications here. The method begins by noting that if one takes the source function to be the mean intensity J, then the equation of transfer can be written in terms of the specific intensity alone. However, the resulting equation is an integrodifferential equation. That is, the intensity, which is a function of the two variables \(\mu\) and \(\tau\), appears differentiated with respect to one of them and is integrated over the other. Thus, \[\mu \frac{d I(\mu, \tau)}{d \tau}=I(\mu, \tau)-\frac{1}{2} \int_{-1}^{+1} I\left(\mu^{\prime}, \tau\right) d \mu^{\prime}\label{10.2.18}\]

Now, as we did in the integral equation for the source function, we can replace the integral by a quadrature summation so that \[\mu \frac{d I(\mu, \tau)}{d \tau}=I(\mu, \tau)-\frac{1}{2} \sum_{j=1}^n I\left(\tau, \mu_j\right) a_j\label{10.2.19}\]

Here the \(a_\mathrm{j}\) values are the weights of the quadrature scheme. This is a functional differential equation for \(\mathrm{I}(\tau,\mu)\) in terms of the solution at certain discrete values of \(\mu_\mathrm{i}\). Chandrasekhar¹ is very explicit about using a gaussian quadrature scheme; a scheme that yields exact answers for polynomials of degree \(2n-1\) or less utilizes the zeros of the Legendre polynomials of degree n as defined in the interval -1 to +1. A more accurate procedure is to divide the integral in equation \ref{10.2.18} into two integrals, one from -1 to 0 and the other from 0 to +1, and to approximate these integrals separately. The reason for this is that, since there is no incident radiation, the intensity develops a discontinuity in \(\mu\) at \(\tau=0\). Numerical quadrature schemes rely on the function to be integrated, in this case \(\mathrm{I}(\mu,\tau)\), being well approximated by a polynomial throughout the range of the integral. Splitting the integral at the discontinuity allows the resulting integrals to be well approximated where the single integral cannot be. This procedure is sometimes called the double-gauss quadrature scheme. However, this “engineering detail” in no way affects the validity of the basic approach.

As we did with equation \ref{10.2.10}, we evaluate the functional equation of transfer [equation \ref{10.2.19}] at the same values of \(\mu\) as are used in the summation so that \[\mu_i \frac{d I\left(\tau, \mu_i\right)}{d \tau}=I\left(\tau, \mu_i\right)-\frac{1}{2} \sum_{j=1}^n I\left(\tau, \mu_j\right) a_j\label{10.2.20}\]

We now have a system of \(n\) homogeneous linear differential equations for the functions \(\mathrm{I}(\tau,\mu_\mathrm{i})\). Each of these functions represents the specific intensity along a particular direction specified by the value of \(\mu_\mathrm{i}\). Since the value \(\mu_\mathrm{i}=0\) represents the point of discontinuity in \(\mathrm{I}(\mu,\tau)\) at the surface, this value should be avoided. Thus, there will normally be as many negative values of \(\mu_\mathrm{i}\) as positive ones. To solve the problem, we must find n constants of integration for the \(n\) first-order differential equations.

Inspired by the general exponential attenuation of a beam of photons passing through a medium, let us assume a solution of the form \[I\left(\tau, \mu_i\right)=g_i e^{-k \tau}\label{10.2.21}\]

Substitution of this form into this set of linear differential equations \ref{10.2.20}, will satisfy the equations if \[\frac{S_\nu\left(\tau_\nu\right)}{I_\nu(0,1)} \approx \frac{B_\nu\left[T\left(\tau_\nu\right)\right]}{I_\nu(0,1)}=\sum_{i=0}^n \beta_i \tau_\nu^i=\sum_{i=0}^n \beta_i\left(\frac{\tau_\nu}{\mu}\right)^i \mu^i\label{10.2.22}\]

and \(k\) satisfies the eigenvalue equation \[1=\sum_{j=1}^{n / 2} \frac{a_j}{1-\mu_j^2 k^2} \quad \mu_j>0\label{10.2.23}\]

Thus equation \ref{10.2.22} provides a constant of integration for every distinct value of \(k\). Since in all quadrature schemes the sum of the weights must equal the interval, \(k^2=0\) will satisfy equation \ref{10.2.23}. Thus, since equation \ref{10.2.23} is essentially polynomic in form there will be \(n/2 - 1\) distinct nonzero values of \(k^2\) and thus \(n - 2\) distinct nonzero values of \(k\) which we denote as \(\pm k_\alpha\). When these are combined with the value \(k = 0\), we are still missing one constant of integration. Wick, inspired by the Eddington approximation, suggested a solution of the form \[I\left(\tau, \mu_i\right)=b\left(\tau+q_i\right)\label{10.2.24}\]

Substitution of this form into equations \ref{10.2.20} also satisfies the equation of transfer provided that \[q_i=\mu_i+Q\label{10.2.25}\]

The product constant \(b\mathrm{Q}\) can be identified with the constant obtained from \(k^2=0\) so it cannot be regarded as a new constant of integration; but the term \(b\tau\) can be regarded as such and therefore completes the solution, so that \[I\left(\tau, \mu_i\right)=b\left[\sum_{\alpha=1}^m\left(\frac{L_{+\alpha} e^{-k_\alpha \tau}}{1+\mu_i k_\alpha}+\frac{L_{-\alpha} e^{+k_\alpha \tau}}{1-\mu_i k_\alpha}\right)+\mu_i+\tau+Q\right]\label{10.2.26}\]

where \[m=\frac{n}{2}-1\label{10.2.27}\]

and the values of \(\mu_\mathrm{i}\) range from -1 to +1. The constants \(\mathrm{L}_{\pm\alpha}\) are the constants that result from equation \ref{10.2.22} and the distinct values of \(k_\alpha\).

Moments of the Radiation Field from Discrete Ordinates We can generate the moments of the radiation field at a level of approximation which is consistent with the solution given by equation \ref{10.2.26} by using the same quadrature scheme for the evaluation of the integrals over m that was used to replace the integral in the integrodifferential equation of radiative transfer. Thus, \[J(\tau)=\frac{1}{2} \int_{-1}^{+1} I(\tau, \mu) d \mu=\frac{1}{2} \sum_{i=1}^n I\left(\tau, \mu_i\right) a_i\label{10.2.28}\]

We already have the values \(\mathrm{I}(\tau,\mu_\mathrm{i})\) required to evaluate the resulting sums. For the gaussian quadrature schemes suggested, the \({a_\mathrm{i}}’\mathrm{s}\) are symmetrically distributed in the interval -1 to +1, while the \({\mu_\mathrm{i}}’\mathrm{s}\) are antisymmetrically distributed. Making use of these facts, substituting the solution [equation \ref{10.2.26}] into equation \ref{10.2.28}, and manipulating, we get \[J(\tau)=b\left[\sum_{\alpha=1}^m\left(L_{+\alpha} e^{-k_\alpha \tau}+L_{-\alpha} e^{+k_\alpha \tau}\right)+\tau+Q\right]\label{10.2.29}\]

Following the same procedure for the flux, we get \[F(\tau)=\frac{4 b}{3}=\mathrm{const}\label{10.2.30}\]

so that the constant \(b\) of the Wick solution is related to the constant flux. All that remains to complete the solution is to determine the constants \(\mathrm{L}_{\forall \alpha}\) from the boundary conditions.

Application of Boundary Values to the Discrete Solution At no point in the derivation have we used of the fact that the atmosphere is assumed semi-infinite. So, in principle, the solution given by equation \ref{10.2.26} is correct for finite slabs. Some applications of the approach have been used in the study of planetary atmospheres, and so for generality let us consider the application to an atmosphere which has a finite thickness \(\tau_0\). For such an atmosphere, we must know the distribution of the intensity entering the atmosphere at the base \(\tau_0\) as well as that which is incident on the surface. Given that, it is a simple matter to equate the solution [equation \ref{10.2.26}] to the boundary values, and we get \[\begin{aligned}
I\left(-\mu_i, 0\right) & =\frac{3}{4} F\left[\sum_{\alpha=1}^m\left(\frac{L_{+\alpha}}{1-\mu_i k_\alpha}+\frac{L_{-\alpha}}{1+\mu_i k_\alpha}\right)-\mu_i+Q\right] \\
I\left(+\mu_i, \tau_0\right) & =\frac{3}{4} F\left[\sum_{\alpha=1}^m\left(\frac{L_{+\alpha} e^{-k_\alpha \tau_0}}{1+\mu_i k_\alpha}+\frac{L_{-\alpha} e^{+k_\alpha \tau_0}}{1-\mu_i k_\alpha}\right)+\mu_i+\tau_0+Q\right] \\
i & =1, \ldots, \frac{n}{2}
\end{aligned}\label{10.2.31}\]

These equations represent n equations in n unknowns. There are \(2n-2\) values of \({\mathrm{L}_{\pm\alpha}}’\mathrm{s}\), \(\mathrm{F}\), and \(\mathrm{Q}\) all specified by the n values of the boundary intensity. Here we explicitly incorporated the sign of \(\mu_\mathrm{i}\) into the equation so that all values of \(\mu_\mathrm{i}\) should be taken to be positive. Although the equations are effectively linear in the unknowns, note that the coefficients of those equations grow exponentially with optical depth. Indeed, since the nonzero values of \(k_\alpha\) are all greater than unity, that growth is quite rapid. In practice, it is virtually impossible to solve these equations for any value of \(\tau_0>100\). Indeed, if the order of approximation is large, the practical upper limit is nearer 10. This instability is inherent in all discrete ordinate methods used for finite atmospheres.

The reason is fairly straightforward. Each of the \({k_\alpha}’\mathrm{s}\) corresponds to a stream of radiation with a particular value of \(\mu_\mathrm{i}\). The total optical path for this radiation stream is \(\tau_0/\mu_\mathrm{i}\). Since the solution of equation \ref{10.2.26} is essentially a linear two-point boundary-value problem, the solution at one boundary is determined by the solution at the other boundary. If part of the solution at one boundary is optically remote from the other boundary, it will decouple from the solution, causing the solution to become singular or poorly determined. Physically, the photons from the remote boundary have been so randomized by scatterings or absorptions that all information pertaining to their direction of entrance into the atmosphere has been lost. In the case of the semi-infinite atmosphere, this has explicitly been taken into account, and the information from the lower boundary is contained in the finite and constant radiative flux.

We can see the effect of this constraint on the discrete solution by examining the behavior of the solution [equation \ref{10.2.31}] as \(\tau_0\rightarrow4\). Since we require the radiation field to remain finite as \(\tau_0\rightarrow4\), the \({\mathrm{L}_{-\alpha}}’\mathrm{s}\) must go to zero. Thus, the influence of the deep radiation field explicitly disappears from the solution, and the radiative flux becomes the eigenvalue of the problem. So the complete solution for the semi-infinite gray atmosphere for the method of discrete ordinates is \[\begin{aligned}
& B(\tau)=J(\tau)=\frac{3}{4} F\left(\sum_{\alpha=1}^m L_{+\alpha} e^{-k_\alpha \tau}+\tau+Q\right) \\
& I\left(-\mu_i, 0\right)=0=\frac{3}{4} F\left(\sum_{\alpha=1}^m \frac{L_{+\alpha}}{1-\mu_i k_\alpha}-\mu_i+Q\right) \\
& F=\text { const }=\frac{\sigma T_e^4}{\pi} \\
& I=\sum_{j=1}^{m+1} \frac{a_j}{1-\mu_j^2 k_\alpha^2}
\end{aligned}\label{10.2.32}\]

Table 10.1 contains some values of \(\mathrm{L}_{+\alpha}\), \(k_\alpha\), and \(\mathrm{Q}\) for various orders of approximation for the semi-infinite gray atmosphere for the single-gauss quadrature scheme. By analogy to the Eddington approximation, the source function is sometimes written as \[J(\tau)=\frac{3}{4} F[\tau+q(\tau)]\label{10.2.33}\]

where \[q(\tau)=Q+\sum_{\alpha=1}^m L_{+\alpha} e^{-k_\alpha \tau}\label{10.2.34}\]

is known as the Hopf function. It is clear that for the Eddington approximation the appropriate Hopf function would be \(q(\tau)=2/3\). The Eddington approximation also avoids the problem of the solution's becoming unstable with increasing depth, by the use of the diffusion approximation, which basically assumes that the radiation field has been directionally randomized.

Nonconservative Gray Atmospheres The notion of a nonconservative gray atmosphere may sound like a contradiction in terms, and if it were meant to apply to all frequencies, it would be. However, consider the case where the opacity is essentially gray over the part of the spectrum containing most of the emergent radiation, but radiative equilibrium does not apply because some energy is lost from the radiation field to perhaps convection. Or consider an atmosphere where the dominant opacity source is the scattering of light from a hot external source, but the atmosphere itself is so cold that the thermal emission can be neglected. Planetary atmospheres often fit into this category.

Table 10.1

Under these conditions, the equation of transfer becomes \[\mu \frac{d I(\tau, \mu)}{d \tau}=I(\tau, \mu)-\frac{p}{2} \int_{-1}^{+1} I\left(\tau, \mu^{\prime}\right) d \mu^{\prime}\label{10.2.35}\]

which, by the same methods used to generate equation \ref{10.2.23}, yields the eigenvalue equation \[1=p \sum_{j=1}^{n / 2} \frac{a_j}{1-\mu_j^2 k_\alpha^2} \quad \mu_j>0\label{10.2.36}\]

Here \(p\) is the scattering albedo, or the fraction of interacting photons that are scattered. Since \(p<1\) for a nonconservative atmosphere, there will now be \(n\) distinct \({k_\alpha}’\mathrm{s}\) and \(n\) distinct \(\mathrm{L_{\pm\alpha}}’\mathrm{s}\), so that the \(n\) values of the boundary radiation field completely specify the solution. The \(\mathrm{L}_{\forall \alpha} \text { 's }\) are specified by the boundary equations \[\begin{aligned}
& I\left(-\mu_i, 0\right)=\sum_{\alpha=1}^{n / 2}\left(\frac{L_{+\alpha}}{1-\mu_i k_\alpha}+\frac{L_{-\alpha}}{1+\mu_i k_\alpha}\right) \\
& I\left(+\mu_i, 0\right)=\sum_{\alpha=1}^{n / 2}\left(\frac{L_{+\alpha}}{1+\mu_i k_\alpha}+\frac{L_{-\alpha}}{1-\mu_i k_\alpha}\right)
\end{aligned}\label{10.2.37}\]

and the source function for the atmosphere is given by \[S(\tau)=J(\tau)=\sum_{\alpha=1}^{n / 2} L_{+\alpha} e^{-k_\alpha \tau}+L_{-\alpha} e^{+k_\alpha \tau}\label{10.2.38}\]

We need not consider the unilluminated semi-infinite atmosphere since all radiation moving up through a nonconservative semi-infinite atmosphere will eventually be lost before it emerges. Thus, only the finite slab or an illuminated semi-infinite nonconservative atmosphere will yield anything other than the trivial solution.