10.2: Classical Solution to the Equation of Radiative Transfer and Integral Equations for the Source Function

There are basically two schools of approach to the solution of the equation of transfer. One involves the solution of an integral equation for the source function, while the other deals directly with the differential equation of transfer. Both have their merits and drawbacks. Since both are widely used, we give examples of each. Both involve the classical solution, so that we begin the discussion with that solution.

a. Classical Solution of the Equation of Transfer for the Plane-Parallel Atmosphere

The equation of transfer is a linear differential equation, which implies that a formal solution exists for the radiation field in terms of the source function. This linear property is a marked difference from the situation in stellar interiors where the structure equations were all highly nonlinear. Although under some conditions the solution [i.e., \(\mathrm{I}_\nu\left(\tau_\nu, \mu\right)\)] itself is involved in the source function, this involvement is still linear. Let us consider a fairly general equation of radiative transfer for a plane-parallel atmosphere, but one where we may neglect time-dependent effects and the presence of the potential gradient on the radiation field. \[\mu \frac{d I_\nu\left(\tau_\nu, \mu\right)}{d \tau_\nu}=I_\nu\left(\tau_\nu, \mu\right)-S_\nu\left(\tau_\nu, \mu\right)\label{10.1.1}\]

Since this equation is linear in \(\mathrm{I}_\nu(\tau_\nu,\mu)\), we may write the complete solution as the sum of the solution to the homogeneous equation plus any particular solution. So let us choose as homogeneous and particular solutions \[I_h\left(\mu, \tau_\nu\right)=c_1 e^{-\tau_\nu / \mu}+c_2 e^{+\tau_\nu / \mu}, \quad I_p\left(\mu, \tau_\nu\right)=f\left(\tau_\nu\right) e^{\tau_\nu / \mu}\label{10.1.2}\]

Substitution into the equation of transfer places constraints on c₁ and \(f’(\tau_\nu)\), namely \[c_1=0 \quad f^{\prime}\left(\tau_\nu\right)=\frac{-S_\nu\left(\tau_\nu\right) e^{-\tau_\nu / \mu}}{\mu}\label{10.1.3}\]

While we have assumed that the geometry of the atmosphere is plane-parallel, we have not yet specified the extent of the atmosphere. For the moment, let us assume that the atmosphere consists of a finite slab of thickness. \(\tau_0\) (see Figure 10.1). The general classical solution for the plane-parallel slab is then \[I_\nu\left(\mu, \tau_\nu\right)=c_2 e^{\tau_\nu / \mu}-\int S(t) e^{-\left(t-\tau_\nu\right) / \mu} \frac{d t}{\mu}\label{10.1.4}\]

0) directed streams of radiation. The boundary conditions necessary for the solution are specified at \(\tau_v=0\), and \(\tau_v=\tau_0\)." style="width: 783px; height: 487px;" width="783px" height="487px" src="/@api/deki/files/66328/Screenshot_2026-01-28_at_12.28.32%25E2%2580%25AFAM.png">

Figure 10.1 shows the geometry for a plane-parallel slab. Note that there are inward (µ<0) and outward (µ>0) directed streams of radiation. The boundary conditions necessary for the solution are specified at \(\tau_\nu=0\), and \(\tau_\nu=\tau_0\).

Since the equation of transfer is a first order linear equation, only one constant must be specified by the boundary conditions. However, even though the depth variable \(\tau_\nu\) is the only independent variable that appears in a derivative, we must always remember that \(\mathrm{I}_\nu(\tau_\nu,\mu)\) is a function of the angular variable \(\mu\). Thus in general, the constant of integration \(c_2\) will depend on the direction taken by the radiation. For radiation flowing outward in the atmosphere (that is, \(\mu > 0\)), the constant \(c_2\) will be set equal to the radiation field at the base of the atmosphere [that is, \(\mathrm{I}_\nu(\mu,\tau_0)\)] and the integral will include the contribution from the source function from all depths ranging from \(\tau_0\) to the point of interest \(\tau_\nu\). If we were concerned about radiation flowing into the atmosphere (that is, \(\mu<0\)), then the integral in equation \ref{10.1.4} would cover the interval from 0 to \(\tau_\nu\) and \(c_2\) would be chosen equal to the incident radiation field \([\mathrm{I}_\nu(-\mu,0)]\).

At this point we encounter one of the notational problems that often leads to confusion in understanding the literature in radiative transfer. For most problems in stellar atmospheres, there is a significant difference between the radiation field represented by the inward-directed streams of radiation and that represented by those flowing outward. In modeling the normal stellar atmosphere, there is no incident radiation present so that the incident intensity \(\mathrm{I}_\nu(-\mu,0)=0\). However, the outward-directed streams always result from a lower boundary condition which is nonzero. Thus it is useful to distinguish between the inward- and outward-directed streams in some notational way. We have already used a standard method of indicating this difference; namely, we explicitly labeled the inward-directed streams by \(-\mu\). Thus, we usually regard the angular variable m as an intrinsically positive quantity that is bounded by \(0<\mu<1\). The sign of m must then be explicitly indicated, and we do this when we use this convention. Thus, to gain a physical understanding of the meaning of any solution for the radiation field, one must always keep in mind which streams of radiation are being considered.

The general classical solution for the two streams can then be written as \[\begin{array}{ll}
I_\nu\left(+\mu, \tau_\nu\right)=-\int_{\tau_0}^{\tau_\nu} S(t) e^{\left(\tau_\nu-t\right) / \mu} \frac{d t}{\mu}+I_\nu\left(+\mu, \tau_0\right) e^{\left(\tau_\nu-\tau_0\right) / \mu} & \mu \geq 0 \\
I_\nu\left(-\mu, \tau_\nu\right)=\int_0^{\tau_\nu} S(t) e^{-\left(\tau_\nu-t\right) / \mu} \frac{d t}{\mu}+I_\nu(-\mu, 0) e^{-\tau_\nu / \mu} &
\end{array}\label{10.1.5}\]

White \(\tau_\nu\) represents the vertical depth in the atmosphere increasing inward, \(\tau_\nu/\mu\) is the actual path along the direction taken by the radiation. In general, extinction by scattering or absorption will exponentially diminish the strength of the intensity by \(\mathrm{e}^{-\tau/\mu}\). Since the source function represents the local source of photons from all processes, and since it is attenuated by the optical distance along the path of the radiation, the integrand of the integral represents the local contribution of the source function to the value of the intensity at \(\tau_\nu\). The remaining term simply represents the local contribution to the specific intensity of the attenuated incident radiation.

One further complication must be dealt with before we can use this description of a stellar atmosphere. In general, stellar atmospheres can be regarded as being infinitely thick. Since the influence of the lower boundary diminishes as \(\mathrm{e}^{\tau-\tau_0}\), and since this optical depth will exceed several hundred within a few thousand kilometers of the surface for main sequence stars, we can take it to be infinity. In addition, we should require the radiative flux to be finite everywhere. This will force the constant \(c_2\) in equation \ref{10.1.4} to vanish. Furthermore, the surface is generally unilluminated. So we can write the classical solution for the semi-infinite plane-parallel atmosphere as \[\begin{aligned}
& I_\nu\left(+\mu, \tau_\nu\right)=-\int_{\infty}^{\tau_\nu} S(t) e^{+\left(\tau_\nu-t\right) / \mu} \frac{d t}{\mu} \\
& I_\nu\left(-\mu, \tau_\nu\right)=\int_0^{\tau_\nu} S(t) e^{-\left(\tau_\nu-t\right) / \mu} \frac{d t}{\mu}
\end{aligned}\label{10.1.6}\]

b. Schwarzschild-Milne Integral Equations

One reason that the equation of transfer admits such a simple solution compared to the equations of stellar structure is that we have confined most of the difficult physics to the source function. What is left is largely geometry and hence affords a simple solution. However, the classical solution does allow for the generation of the entire radiation field should it be possible to specify the source function. It also allows us to remove the explicit structure of the radiation field and to generate an expression for the source function itself. The result is an integral equation, that is, an equation where the unknown appears under the integral sign as well as outside it.

While much attention has been paid to the solution of differential equations, far less has been given to integral equations. However, it is very often numerically more efficient and accurate to solve an integral equation as opposed to the corresponding differential equation. Therefore, we spend some time and effort with these integral equations, for they provide a very productive path toward the solution of problems in radiative transfer.

Integral Equation for the Source Function In Chapter 9 we showed that, for coherent isotropic scattering, we could write a quite general expression for the source function [equation \ref{9.2.33}]. If we re-express that result in terms of the mean intensity, we get \[S_\nu\left(\tau_\nu\right)=\epsilon_\nu B_\nu\left[T\left(\tau_\nu\right)\right]+\left(1-\epsilon_\nu\right) J_\nu\left(\tau_\nu\right)\label{10.1.7}\]

where \[\epsilon_\nu=\frac{\kappa_\nu}{\kappa_\nu+\sigma_\nu}\label{10.1.8}\]

Now the role of the classical solution becomes evident. The source function contains the mean intensity \(\mathrm{J}_\nu(\tau_\nu)\), which can be generated from the classical solution that contains the source function itself. Thus, if we substitute the classical solution [equation \ref{10.1.6}] into the definition for \(\mathrm{J}_\nu(\tau_\nu)\) [equation \ref{9.3.2}], we get \[\begin{aligned}
J_\nu\left(\tau_\nu\right)= & \frac{1}{2} \int_0^1 I\left(+\mu^{\prime}, \tau_\nu\right) d \mu^{\prime}+\frac{1}{2} \int_0^1 I\left(-\mu^{\prime}, \tau_\nu\right) d\left(-\mu^{\prime}\right) \\
J_\nu\left(\tau_\nu\right)= & -\frac{1}{2} \int_0^1\left[\int_{\infty}^{\tau_\nu} S(t) e^{+\left(\tau_\nu-t\right) / \mu^{\prime}} \frac{d t}{\mu^{\prime}}\right] d \mu^{\prime} \\
& +\frac{1}{2} \int_0^1\left[\int_0^{\tau_\nu} S(t) e^{-\left(\tau_\nu-t\right) / \mu^{\prime}} \frac{d t}{\mu^{\prime}}\right] d \mu^{\prime}
\end{aligned}\label{10.1.9}\]

Now notice that the argument of the exponential is always negative and that the two integrals over t are contiguous. Thus, we can combine these integrals into a single integral that ranges from 0 to 4. In addition, \(\mathrm{t}\) and \(\mu\) are independent variables so that we may interchange the order of integration and get \[J_\nu\left(\tau_\nu\right)=\frac{1}{2} \int_0^{\infty} S(t)\left[\int_0^1 e^{-\left|\tau_\nu-t\right| / \mu^{\prime}} \frac{d \mu^{\prime}}{\mu^{\prime}}\right] d t\label{10.1.10}\]

The quantity in brackets is a well-known function in mathematical physics known as the exponential integral. It depends only on the independent variables of the problem and therefore can be regarded as a largely geometric function. Its formal definition is \[E_n(z)=\int_1^{\infty} \frac{e^{-z t} d t}{t^n}=\int_0^1 e^{-z / y} y^{n-2} d y=\int_0^1 \frac{e^{-z / y} y^{n-1} d y}{y}\label{10.1.11}\]

and when expressed by the final integral, it has the same form as the integral in brackets in equation \ref{10.1.10}. While the exponential integral may not be terribly familiar, it should be regarded with no more fear and trepidation than sines and cosines. There is an entire set of these functions where each member is denoted by n, and they have a single argument, which for our purposes will be confined to the real line. These functions (except for the first exponential integral at the origin) are well behaved and resemble \(e^{\mathrm{-x/(nx)}}\) for large x. Some useful properties of exponential integrals are \[\begin{aligned}
& n E_{n+1}(x)=e^{-x}-x E_n(x) \quad n>1 \\
& \frac{d E_{n+1}(x)}{d x}=-E_n(x) \\
& E_n(0)=\frac{1}{n-1} \quad E_n(\infty)=0
\end{aligned}\label{10.1.12}\]

Making use of the first exponential integral, we can rewrite our expression for the mean intensity [equation \ref{10.1.10}] as \[J_\nu\left(\tau_\nu\right)=\frac{1}{2} \int_0^{\infty} S(t) E_1\left|\tau_\nu-t\right| d t\label{10.1.13}\]

Combining this with equation \ref{10.1.7} for the source function, we arrive at the desired integral equation for the source function: \[S_\nu\left(\tau_\nu\right)=\epsilon_\nu B_\nu\left[T\left(\tau_\nu\right)\right]+\left(1-\epsilon_\nu\right) \frac{1}{2} \int_0^{\infty} S_\nu(t) E_1\left|\tau_\nu-t\right| d t\label{10.1.14}\]

Any function that multiplies the unknown in the integrand of an integral equation is called the kernel of the integral equation. Thus, the first exponential integral is the kernel of the integral equation for the source function. The connection between the physical state of the gas and the source function is contained in the term that makes the equation inhomogeneous, namely, the one involving the Planck function \(\mathrm{B}_\nu[T(\tau_\nu)]\). A solution of this equation, when combined with the classical solution, will yield the full solution to the radiative transfer problem since \(\mathrm{I}_\nu(\mu,\tau_\nu)\) will be specified for all values of \(\mu\) and \(\tau_\nu\).

It is possible to understand equation \ref{10.1.14} from a physical standpoint. Now \(\varepsilon(\tau_\nu)\) is the fraction of locally generated photons that arise from thermal processes, so that the first term is simply the local contribution to the source function from thermal properties of the gas. The second term represents the contribution from scattering. We have already said that a fundamental aspect of stellar atmospheres is the dependence of the local radiation field on the global solution for the radiation field. Nowhere is this more clearly demonstrated than in this term. The scattering contribution to the source function is made up of contributions from the source function throughout the atmosphere. However, these contributions decline with increasing distance from the point of interest, and they decline roughly exponentially.

One may object that this integral equation is a very specialized equation since it relies on the source function's being expressible in terms of the mean intensity and therefore is valid only for isotropic scattering. However, consider the very general expression for the source function given by equation \ref{9.2.27}. As long as the angular dependence of the redistribution function is known, it will be possible to carry out the integrals over solid angle and express the source function as a combination of the moments of the radiation field. As long as this can be done, the appropriate moments can be generated from the classical solution for the equation of transfer which will, in turn, involve only the source function. Thus, the moments can be eliminated from the moment expression for the source function, yielding an integral equation. To be sure, this will be a more complicated integral equation, but it will still be solvable by the same techniques that we apply to equation \ref{10.1.14}. Thus, the existence of an integral equation for the source function is a quite general result and represents the separation of the depth dependence of the radiation field from the angular dependence, which can be obtained from the classical solution.

Integral Equations for Moments of the Radiation Field Useful as the integral equation for the source function is, it is often convenient to have similar expressions for the moments of the radiation field. We should not be surprised that such expressions exist since the angular moments are free, by definition, of the angular dependence characteristic of the classical solution. Indeed, we have already supplied the required expressions to obtain an integral equation for the mean intensity. We simply use equation \ref{10.1.7} to eliminate \(S_\nu(\mathrm{t})\) from equation \ref{10.1.13}, and we have \[\begin{aligned}
J_\nu\left(\tau_\nu\right)= & \frac{1}{2} \int_0^{\infty} \epsilon_\nu(t) B_\nu[T(t)] E_1\left|\tau_\nu-t\right| d t \\
& +\frac{1}{2} \int_0^{\infty}\left[1-\epsilon_\nu(t)\right] J_\nu(t) E_1\left|\tau_\nu-t\right| d t
\end{aligned}\label{10.1.15}\]

It is now clear how to develop similar expressions for the remaining moments, since equation \ref{10.1.13} was obtained by taking moments of the classical solution to the equation of transfer. Let us define an operator which is commonly used to represent this process. \[\Lambda_n\left(\tau_\nu\right)|G(t)| \equiv \int_0^{\infty}\left(\frac{\left(t-\tau_\nu\right)}{\left|t-\tau_\nu\right|}\right)^{n+1} E_n\left|\tau_\nu-t\right| G(t) d t\label{10.1.16}\]

The \(\Lambda_n\) operator is an integral operator which operates on a function by employing an exponential integral kernel. The term in large parentheses simply denotes the sign of the kernel throughout the region. With this integral operator, we can express the first three moments of the radiation field in terms of the source function as follows: \[\begin{aligned}
J_\nu\left(\tau_\nu\right) & =\frac{1}{2} \Lambda_1\left(\tau_\nu\right)\left|S_\nu(t)\right| \\
H_\nu\left(\tau_\nu\right) & =\frac{1}{2} \Lambda_2\left(\tau_\nu\right)\left|S_\nu(t)\right| \\
K_\nu\left(\tau_\nu\right) & =\frac{1}{2} \Lambda_3\left(\tau_\nu\right)\left|S_\nu(t)\right|
\end{aligned}\label{10.1.17}\]

Such equations are known as Schwarzschild-Milne type of equations and are extremely useful for the construction of model stellar atmospheres. For example, consider the condition of radiative equilibrium where it is necessary to know the radiative flux throughout the atmosphere, but not the complete radiation field. This information can be obtained directly with the aid of the flux equation of equations \ref{10.1.17} and the source function. Thus, determination of the source function provides a complete solution of the radiative transfer problem.

c. Limb-darkening in a Stellar Atmosphere

There is one property of the classical solution of the equation of transfer that we should address before moving on. If we consider the classical solution for the emergent intensity, we see that it basically represents the Laplace transform of the source function, namely \[I(\mu, 0)=\int_0^{\infty} S(t) e^{-t / \mu} \frac{d t}{\mu}=\frac{\mathscr{L}[S(t)]}{\mu}\label{10.1.18}\]

where \(\mathscr{L}[S(t)]\) is the Laplace transform of the source function. Thus determination of the angular distribution of the emergent intensity is equivalent to determining the behavior of the source function with depth. Since the source function is determined by the temperature, determination of the depth dependence of the source function is equivalent to determining the depth dependence of the temperature. This is of considerable significance for stars where this dependence can be measured directly for it provides a direct observational check on the models of those stellar atmospheres.

If we anticipate some later results and assume that the source function can be approximated by \[S(t)=a t+b\label{10.1.19}\]

then \[I(\mu, 0)=a \mu+b\label{10.1.20}\]

Thus, the coefficient \(a\) that multiplies the angular parameter \(\mu\) in the emergent intensity is a direct measure of the source function gradient, while the constant term \(b\) denotes the value of the source function at the boundary. The decrease in brightness as one approaches the limb of the apparent stellar disk implied by equation \ref{10.1.20} is called limb-darkening. Since for spherical stars the variation across the apparent disk is the same as the local angular dependence of the emergent intensity, measurement of the limb-darkening coefficient a yields a measurement of the source function gradient. This is of particular interest for the sun where such measurements are possible. Unfortunately, the poorest theoretical representation of the model atmosphere occurs near the surface, and this corresponds to just that region of the stellar disk (i.e., near the limb where \(\mu\rightarrow0\)) where confirmatory measurements are most difficult to make. Although we have made an approximation to the depth dependence of the source function in equation \ref{10.1.19}, the approximation is unnecessary and more rigorous studies of this depth dependence would deal directly with the Laplace transform itself as given by equation \ref{10.1.18}. We have now compiled methods by which we can theoretically relate the emergent intensity to the source function and provided a potential observational method to verify our result. However, before discussing methods for the solution for the integral equation for the source function [equation \ref{10.1.14}] we consider the solutions to a somewhat simpler problem, in order to gain an appreciation for the behavior of these solutions.

Empirical Determination of \(\boldsymbol{T(\tau_\nu)}\) for the Sun In the sun and some eclipsing binary stars, it is possible to determine the variation of the specific intensity across the apparent disk. If we approximate that variation by \[\frac{I_\nu(0, \mu)}{I_\nu(0,1)} \approx \sum_{i=0}^n a_i \mu^i\label{10.1.21}\]

we can use equation \ref{10.1.18} to obtain a power series representation of the source function with optical depth. Let us further assume that the source function can be represented by the Planck function, which in turn can be expanded in a power series in the optical depth so that \[\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.1.22}\]

Then the substitution of this power series representation into equation \ref{10.1.18} yields \[\sum_{i=0}^n \beta_i i!\mu^i=\sum_{i=0}^n a_i \mu^i\label{10.1.23}\]

Since for the sun, the \({a_\mathrm{i}}’\mathrm{s}\) and \(\mathrm{I}_\nu(0,1)\) may be determined from observation, the \({b_\mathrm{i}}’\mathrm{s}\) may be regarded as known. Thus, the temperature variation with monochromatic optical depth may be recovered from \[\frac{2 h \nu^3 / c^2}{e^{h \nu /\left[k T\left(\tau_\nu\right)\right]}-1}=I_\nu(0,1) \sum_{i=0}^n \beta_i \tau_\nu^i\label{10.1.24}\]

In the sun, the assumption that \(S_\nu\left(\tau_\nu\right)=B_\nu\left(\tau_\nu\right)\) is a particularly good one, so that for the sun the optical depth variation of the temperature can be determined with the same sort of accuracy that attends the determination of the limb-darkening.

Empirical Determination of \(\boldsymbol{\kappa(\tau_1)/\kappa(\tau_2)}\) for the Sun This type of analysis can be continued under the above assumptions to obtain the variation with optical depth of the ratio of two monochromatic absorption coefficients. Since by definition \[d \tau_\nu=-\kappa_\nu \rho d x\label{10.1.25}\]

the ratio of two monochromatic optical depths is \[\frac{d \tau_1}{d \tau_2}=\frac{\kappa_1\left(\tau_\nu\right)}{\kappa_2\left(\tau_\nu\right)}\label{10.1.26}\]

Differentiating equation \ref{10.1.22} with respect to temperature and substituting the result into equation \ref{10.1.26}, we get \[\frac{\kappa_1\left(\tau_\nu\right)}{\kappa_2\left(\tau_\nu\right)}=\frac{I_2(0,1)}{I_1(0,1)} \frac{d B\left(\nu_1\right) / d T}{d B\left(\nu_2\right) / d T} \frac{\sum_i i \beta_i\left(\tau_2\right) \tau_2^{(i-1)}}{\sum_i i \beta_i\left(\tau_1\right) \tau_1^{(i-1)}}\label{10.1.27}\]

Thus it is possible to determine the approximate wavelength dependence of the opacity for stars like the sun from the observed limb-darkening. Such observations provide a valuable check on the theory of stellar atmospheres.