# 2n-Stream Conservative Scattering

W. A. van Wijngaarden  Department of Physics and Astronomy, York University,
Canada  W. Happer  Department of Physics, Princeton University, USA

###### Abstract

We show how to use matrix methods of quantum mechanics to efficiently and
accurately calculate axially symmetric radiation transfer in clouds, with
conservative scattering of arbitrary anisotropy. Analyses of conservative
scattering, where the single scattering albedo is $\tilde{\omega}=1$ and no
energy is exchanged between the radiation and scatterers, began with work by
Schwarzschild, Milne, Eddington and others[1] on radiative transfer in stars.
There the scattering is isotropic or nearly so. It has been difficult to
extend traditional methods to highly anisotropic scattering, like that of
sunlight in Earth’s clouds. The $2n$-stream method described here is a
practical way to handle highly anisotropic, conservative scattering. The basic
ideas of the $2n$-stream method are an extension of Wick’s seminal work on
transport of thermal neutrons by isotropic scattering[2] to scattering with
arbitrary anisotropy. How to do this for finite absorption and
$\tilde{\omega}<1$ was described in a previous paper[3]. But those methods
fail for conservative scattering, when $\tilde{\omega}=1$. Here we show that
minor modifications to the fundamental $2n$-scattering theory for
$\tilde{\omega}<1$ make it suitable for $\tilde{\omega}=1$.

Keywords:  conservative radiative transfer, multiple scattering, phase
functions, equation of transfer

##  1 Introduction

In a previous paper, 2n-Stream Radiative Transfer[3] we outlined how to use
matrix methods of quantum mechanics to accurately and efficiently analyze
radiative transfer with highly anisotropic scattering, like Earth’s cloud
particulates for sunlight. To facilitate subsequent discussions, we will refer
to this paper as WH and to equation ($x$) in it as (WH-$x$).

The theory outlined in WH works as long as there is some absorption and the
single-scattering albedo $\tilde{\omega}$ is less than 1. Some absorption and
emission is always present in nature, so the limitation of WH to
$\tilde{\omega}<1$, is consistent with reality. However, in many important
situations, the absorption can be small compared to elastic scattering. The
limit of no absorption at all and unit single-scattering albedo,
$\tilde{\omega}=1$, is called conservative scattering[4]. Examples of nearly
conservative scattering are Thomson scattering of radiation by free electrons
in stars or Mie scattering of sunlight in clouds. There is a large body of
theoretical literature focused on conservative scattering[1, 4]. This paper
shows how to extend the simple and computationally efficient methods of WH to
the limiting case of conservative scattering.

##  2 Vector notation

As in WH, we use Dirac notation and we assume axially symmetric radiation.
Then the time-independent, steady-state version of the equation of transfer
(WH-58) becomes

$\left(\hat{\mu}\frac{d}{d\tau}+\hat{\eta}\right)|I(\tau)\\}=B(\tau)\hat{\eta}|0)$
(1)

The abstract, $2n\times 1$ vector $|I(\tau)\\}$ approximates the monochromatic
intensity, $I=I(\mu,\tau)$, at the vertical optical depth $\tau$ above the
bottom of a cloud, and with direction cosine $\mu$ with respect to the zenith.
The $2n\times 2n$ direction-cosine matrix is denoted by the symbol $\hat{\mu}$
and the $2n\times 2n$ decay-efficiency matrix is denoted by $\hat{\eta}$. We
will discuss both $\hat{\mu}$ and $\hat{\eta}$ in more detail below. The right
side of (1) describes the isotropic emission of radiation by scattering
particles. The Planck intensity, $B(\tau)$, depends on the radiation frequency
and on the local temperature, $T(\tau)$, as shown in (WH-7). The symbol $|0)$
that multiplies $B(\tau)\hat{\eta}$ is the monopole basis vector $|l)$ of
(WH-35) for $l=0$.

The right multipole basis vector, $|l)$ and its conjugate left basis vector
$(\kern-2.59995pt(l|$, satisfy the orthogonality relation of (WH-32) and the
completeness relation of (WH-33)

$\displaystyle(\kern-2.59995pt(l|l^{\prime})=\delta_{ll^{\prime}},\quad\hbox{and}\quad\sum_{l=0}^{2n-1}|l)(\kern-2.59995pt(l|=\hat{1}.$
(2)

Here $\hat{1}$ denotes a $2n\times 2n$ unit matrix, $|l)$ can be thought of as
a $2n\times 1$ column vector and $(\kern-2.59995pt(l|$ as a $1\times 2n$ row
vector. Using (2), one can write $|I\\}$ in terms of the $2n$ multipole
amplitudes $(\kern-2.59995pt(l|I(\tau)\\}$,

$|I(\tau)\\}=\sum_{l=0}^{2n-1}|l)(\kern-2.59995pt(l|I(\tau)\\},$ (3)

where

$(\kern-2.59995pt(l|I(\tau)\\}=I_{l}(\tau)=\frac{1}{2}\int_{-1}^{1}d\mu\,P_{l}(\mu)I(\mu,\tau).$
(4)

As shown in (WH-31) and (WH-35), the projections of the left and right
multipole basis vectors $(\kern-2.59995pt(l|$ and $|l)$ onto the continuous,
right and left $\mu$-space basis vectors $|\mu\rangle$ and $\langle\mu|$ are

$(\kern-2.59995pt(l|\mu\rangle=\frac{1}{2}P_{l}(\mu),$ (5)

and

$\langle\mu|l)=P_{l}(\mu)(2l+1).$ (6)

Instead of expanding the intensity vector $|I(\tau)\\}$ on $2n$ multipole
basis vectors $|l)$ and $(\kern-2.59995pt(l|$, as in (3), we can expand
$|I(\tau)\\}$ on the $2n$ stream basis vectors, $|\mu_{i})$ and
$(\kern-2.59995pt(\mu_{i}|$, the right and left eigenvectors of the direction
cosine matrix $\hat{\mu}$, defined by (WH-105) as

$\hat{\mu}|\mu_{i})=\mu_{i}|\mu_{i})\quad\hbox{and}\quad(\kern-2.59995pt(\mu_{i}|\hat{\mu}=(\kern-2.59995pt(\mu_{i}|\mu_{i}.$
(7)

In analogy to (2) the right and left stream basis vectors $|\mu_{i})$ and
$(\kern-2.59995pt(\mu_{i}|$ satisfy the orthogonality relation of (WH-86) and
the completeness relation of (WH-87)

$\displaystyle(\kern-2.59995pt(\mu_{i}|\mu_{i^{\prime}})=\delta_{ii^{\prime}},\quad\hbox{and}\quad\sum_{i=1}^{2n}|\mu_{i})(\kern-2.59995pt(\mu_{i}|=\hat{1},$
(8)

In analogy to (3) we can use the completeness relation of (8) to expand the
intensity vector as

$|I(\tau)\\}=\sum_{i=1}^{2n}|\mu_{i})(\kern-2.59995pt(\mu_{i}|I(\tau)\\}.$ (9)

The coefficient $(\kern-2.59995pt(\mu_{i}|I(\tau)\\}$ of (7) is the product of
the Gauss-Legendre weight, $w_{i}$, and the value, $I(\mu_{i},\tau)$, of the
intensity at the Gauss-Legendre direction cosine, $\mu_{i}$.

$(\kern-2.59995pt(\mu_{i}|I(\tau)\\}=w_{i}I(\mu_{i},\tau)=\int_{-1}^{1}d\mu(\kern-2.59995pt(\mu_{i}|\mu\rangle
I(\mu,\tau).$ (10)

Experiments are normally designed to measure unweighted intensities,
$I(\mu_{i},\tau)=w_{i}^{-1}(\kern-2.59995pt(\mu_{i}|I(\tau)\\}$, for which a
common unit is W/(m2 cm-1 sr) or watts per square meter $\times$ wave number
$\times$ steradian. Weighted intensities,
$(\kern-2.59995pt(\mu_{i}|I(\tau)\\}$, are the natural quanties for
mathematical calculations. They are the elements of the $2n\times 1$ intensity
vector, $|I(\tau)\\}$ in $\mu_{i}$-space.

According to (WH-17) and (WH-36), the $2n\times 2n$ elements of the direction-
cosine matrix $\hat{\mu}$ in multipole space are

$\displaystyle(\kern-2.59995pt(l|\hat{\mu}|l^{\prime})$ $\displaystyle=$
$\displaystyle\frac{(l)\delta_{l^{\prime},l-1}+(l+1)\delta_{l^{\prime},l+1}}{2l+1}$
(18) $\displaystyle=$
$\displaystyle\left[\begin{array}[]{lllllll}0&1/1&0&0&0&0&\cdots\\\
1/3&0&2/3&0&0&0&\cdots\\\ 0&2/5&0&3/5&0&0&\cdots\\\ 0&0&3/7&0&4/7&0&\cdots\\\
0&0&0&4/9&0&5/9&\cdots\\\ 0&0&0&0&5/11&0&\cdots\\\
\vdots&\vdots&\vdots&\vdots&\vdots&\vdots&\ddots\\\ \end{array}\right]$

In the second line of (18) the row labels $l$ and column labels $l^{\prime}$
include the integers $0,1,2\ldots,$ $2n-1$. The elements of all but the last
row of the matrix $(\kern-2.59995pt(l|\hat{\mu}|l^{\prime})$ sum to 1. Using
(2) with (18) we note the identities

$(\kern-2.59995pt(0|\hat{\mu}=\sum_{l=0}^{2n-1}(\kern-2.59995pt(0|\hat{\mu}|l)(\kern-2.59995pt(l|=(\kern-2.59995pt(1|,$
(19)

and

$(\kern-2.59995pt(0|\hat{\mu}^{2}=(\kern-2.59995pt(1|\hat{\mu}=\sum_{l=0}^{2n-1}(\kern-2.59995pt(1|\hat{\mu}|l)(\kern-2.59995pt(l|=\frac{(\kern-2.59995pt(0|+2(\kern-2.59995pt(2|}{3}.$
(20)

As shown in connection with (WH-106), the $2n$ different eigenvalues $\mu_{i}$
of (7) are the roots of the Legendre polynomial $P_{2n}$, that is,

$P_{2n}(\mu_{i})=0.$ (21)

We will order the roots such that

$-1<\mu_{1}<\mu_{2}<\cdots<\mu_{2n}<1.$ (22)

We can use the eigenvalues $\mu_{i}$ and the eigenvectors $|\mu_{i})$ and
$(\kern-2.59995pt(\mu_{i}|$ of (7) to write the direction cosine matrix,
$\hat{\mu}$ as

$\hat{\mu}=\sum_{i=1}^{2n}\mu_{i}|\mu_{i})(\kern-2.59995pt(\mu_{i}|.$ (23)

The projections of the left and right multipole basis vectors
$(\kern-2.59995pt(l|$ and $|l)$ onto the right and left stream bases,
$|\mu_{i})$ and $(\kern-2.59995pt(\mu_{i}|$, onto the multipole bases were
give by (WH-84) and (WH-85) as

$(\kern-2.59995pt(l|\mu_{i^{\prime}})=\frac{1}{2}P_{l}(\mu_{i^{\prime}}).$
(24)

and

$(\kern-2.59995pt(\mu_{i}|l)=w_{i}(2l+1)P_{l}(\mu_{i}).$ (25)

We can multiply (24) on the right by $(\kern-2.59995pt(\mu_{i}|$, sum over $i$
and use (8) to find the useful identity

$\displaystyle(\kern-2.59995pt(l|$ $\displaystyle=$
$\displaystyle\sum_{i=1}^{2n}(\kern-2.59995pt(l|\mu_{i})(\kern-2.59995pt(\mu_{i}|$
(26) $\displaystyle=$
$\displaystyle\frac{1}{2}\sum_{i=1}^{2n}P_{l}(\mu_{i})(\kern-2.59995pt(\mu_{i}|$

The projections onto the continuous direction-cosine bases $|\mu\rangle$ and
$\langle\mu|$ of (WH-26) were given by

$\langle\mu|\mu_{i})=\sum_{l=0}^{2n-1}\langle\mu|l)(\kern-2.59995pt(l|\mu_{i})=\sum_{l=0}^{2n-1}\frac{2l+1}{2}P_{l}(\mu)P_{l}(\mu_{i}).$
(27)

and

$(\kern-2.59995pt(\mu_{i}|\mu\rangle=\sum_{l=0}^{2n-1}(\kern-2.59995pt(\mu_{i}|l)(\kern-2.59995pt(l|\mu\rangle=w_{i}\sum_{l=0}^{2n-1}\frac{2l+1}{2}P_{l}(\mu_{i})P_{l}(\mu).$
(28)

One of many equivalent formulas for the weights, $w_{i}$, was given by (WH-72)
as

$\frac{1}{w_{i}}=\sum_{l=0}^{2n-1}\frac{2l+1}{2}P_{l}^{2}(\mu_{i}).$ (29)

Numerical values of the Gauss-Legendre weights, $w_{i}$, and cosines, $\mu$,
are given in pages 916 - 919 of Abramowitz and Stegun[5].

Eq.(1) also contains the efficiency matrix $\hat{\eta}$, defined by (WH-54) in
terms of the scattering-phase matrix $\hat{p}$ and the single-scattering
albedo $\tilde{\omega}$ as

$\hat{\eta}=\hat{1}-\frac{\hat{\omega}}{2}\hat{p}.$ (30)

After a collision between a photon and a cloud particulate or atmospheric
molecule, the fraction of photons scattered, as opposed to being absorbed, is
the single-scattering albedo $\tilde{\omega}$. Since $\tilde{\omega}$ is a
probability, it is constrained to have the values

$0\leq\hat{\omega}\leq 1.$ (31)

The fraction of photons absorbed and converted to heat is $1-\tilde{\omega}$.

Following (WH-40), we write the scattering-phase matrix $\hat{p}$ of (30) as

$\hat{p}=2\sum_{l=0}^{\infty}p_{l}|l)(\kern-2.59995pt(l|.$ (32)

For randomly oriented scatterers, $\hat{p}$ must be diagonal in multipole
space, with matrix elements
$(\kern-2.59995pt(l|\hat{p}|l^{\prime})=2p_{l}\delta_{ll^{\prime}}$. The
coefficients $p_{l}$ are subject to the constraints of (WH-14) and (WH-15),

$p_{0}=1,\quad\hbox{and}\quad|p_{l}|<1,\quad\hbox{for}\quad l>0.$ (33)

We can use (5) and (6) to write the matrix elements for axially symmetric
scattering of radiation with direction cosine $\mu^{\prime}$ to radiation with
direction cosine $\mu$ as

$\displaystyle p(\mu,\mu^{\prime})$ $\displaystyle=$
$\displaystyle\langle\mu|\hat{p}|\mu^{\prime}\rangle$ (34) $\displaystyle=$
$\displaystyle
2\sum_{l=0}^{\infty}p_{l}\langle\mu|l)(\kern-2.59995pt(l|\mu^{\prime}\rangle$
$\displaystyle=$
$\displaystyle\sum_{l=0}^{\infty}p_{l}(2l+1)P_{l}(\mu)P_{l}(\mu^{\prime}),$

the same as (WH-11). We define the phase function for vertical incident
radiation, with $\mu^{\prime}=1$ as

$\displaystyle p(\mu)$ $\displaystyle=$ $\displaystyle p(\mu,1)$ (35)
$\displaystyle=$ $\displaystyle\sum_{l=0}^{\infty}p_{l}(2l+1)P_{l}(\mu)$

The context should prevent any confusion between the phase function
$p(\mu,\mu^{\prime})$ for scattering from $\mu^{\prime}$ to $\mu$ and the
phase function $p(\mu)=p(\mu,1)$ for scattering from $\mu^{\prime}=1$ to
$\mu$.

A simple example is the phase function for polarization-averaged Rayleigh
scattering

$p(\mu)=\frac{3}{4}(1+\mu^{2})=P_{0}(\mu)+\frac{1}{2}P_{2}(\mu).$ (36)

The non-zero coefficients of the multipole expansion (33) are

$p_{0}=1,\quad\hbox{and}\quad p_{2}=\frac{1}{10}.$ (37)

Other important examples are the maximum-forward-scattering phases,
$p(\mu)=\varpi^{\\{p\\}}(\mu)$ of (WH-134), which can be written as

$\displaystyle\varpi^{\\{p\\}}(\mu)$ $\displaystyle=$
$\displaystyle\frac{2(1+\mu)}{p(p+1)}\left[\frac{dP_{p}}{d\mu}(\mu)\right]^{2}$
(38) $\displaystyle=$
$\displaystyle\sum_{l=0}^{2p-1}P_{l}(\mu)(2l+1)\varpi^{\\{p\\}}_{l}.$

Representative multipole coefficients of (38),

$\displaystyle\varpi^{\\{p\\}}_{l}=\frac{1}{2}\int_{-1}^{1}d\mu\,P_{l}(\mu)\varpi^{\\{p\\}}(\mu),$
(39)

are listed in Table 1 of WH.

Substituting (32) into (30) we see that the eigenvectors of $\hat{\eta}$ are
the multipole bases $|l)$ and $(\kern-2.59995pt(l|$,

$\hat{\eta}=\sum_{l=0}^{\infty}\eta_{l}|l)(\kern-2.59995pt(l|,\quad\hbox{with}\quad\hat{\eta}|l)=\eta_{l}|l),\quad\hbox{and}\quad(\kern-2.59995pt(l|\hat{\eta}=(\kern-2.59995pt(\eta|\eta_{l}.$
(40)

In view of (33) and (31), we see that the values of the efficiency eigenvalues
$\eta_{l}$ must be nonnegative

$\eta_{l}=1-\tilde{\omega}p_{l}\geq 0$ (41)

For future reference, we note the special cases of (40) and (41)

$\hat{\eta}|0)=\eta_{0}|0)=(1-\tilde{\omega})|0)\quad\hbox{and}\quad(\kern-2.59995pt(0|\hat{\eta}=(\kern-2.59995pt(0|\eta_{0}=(\kern-2.59995pt(0|(1-\tilde{\omega}).$
(42)

##  3 Conservation of energy

According to (WH-19), the volume energy density of the radiation, $u(\tau)$,
at the optical depth $\tau$ is proportional to the monopole moment
$I_{0}(\tau)$ of the intensity,

$u(\tau)=\frac{2\pi}{c}\int_{-1}^{1}d\mu
I(\mu,\tau)=\frac{4\pi}{c}I_{0}(\tau)=(\kern-2.59995pt(0|I(\tau)\\}.$ (43)

Here the energy density vector $|u(\tau)\\}$ is related to the intensity
vector $|I(\tau)\\}$ by

$|u(\tau)\\}=\frac{4\pi}{c}|I(\tau)\\},$ (44)

where $c$ is the speed of light.

In (WH-20) the vertical flux $Z(\tau)$ of radiation is shown to be
proportional to the dipole moment $I_{1}(\tau)$ of the intensity

$\displaystyle Z(\tau)=2\pi\int_{-1}^{1}d\mu\,\mu I(\mu,\tau)=4\pi
I_{1}=(\kern-2.59995pt(0|Z(\tau)\\}.$ (45)

The flux vector $|Z(\tau)\\}$ is related to the intensity vector $|I(\tau)\\}$
by

$|Z(\tau)\\}=4\pi\hat{\mu}|I(\tau)\\},$ (46)

where $\hat{\mu}$ is the direction cosine matrix of (18). Unlike the energy
density $u(\tau)$ of (43) which is always non-negative, $u(\tau)\geq 0$, the
flux $Z(\tau)$ can be positive for net upward radiation or negative for net
downward radiation.

For conservative scattering, the second moment of the intensity is of special
significance. Following §11 (86) of Chandrasekhar[4], we write the second
moment in terms of the “$K$-integral,”

$\displaystyle
K(\tau)=\frac{1}{2}\int_{-1}^{1}d\mu\,\mu^{2}I(\mu,\tau)=\frac{I_{0}+2I_{2}}{3}=(\kern-2.59995pt(0|K(\tau)\\}.$
(47)

Here the $K$-integral vector $|K(\tau)\\}$ is related to the intensity vector
$|I(\tau)\\}$ by

$|K(\tau)\\}=\hat{\mu}^{2}|I(\tau)\\}.$ (48)

Multiplying (1) on the left by $(\kern-2.59995pt(0|$, using (4) and (42) with
(19) we find

$\frac{d}{d\tau}Z(\tau)+(1-\tilde{\omega})cu(\tau)=4\pi(1-\tilde{\omega})B(\tau).$
(49)

For time-independent conditions, thermal radiation may heat or cool the
atmosphere because of an excess of absorption over emission or vice versa.
Over most of Earth’s troposphere radiation cooling is approximately balanced
by convective heating. The heating rate $R(\tau)$ per unit of optical depth is
defined by

$R(\tau)\,d\tau=Z(\tau)-Z(\tau+d\tau)=-d\tau\frac{dZ(\tau)}{d\tau},\quad\hbox{or}\quad
R(\tau)=-\frac{dZ(\tau)}{d\tau}.$ (50)

The heating $R(\tau)d\tau$ of the infinitesimally thin layer between the
optical depths $\tau$ and $d\tau$ is the excess of the flux $Z(\tau)$ through
the bottom surface of the layer over the flux $Z(\tau+d\tau)$ through the top
surface. Using (49) with (50) we find that the heating rate can be written as

$R(\tau)=(1-\tilde{\omega})[cu(\tau)-4\pi B(\tau)].$ (51)

The absorption rate of energy is $(1-\tilde{\omega})cu(\tau)$ and the emission
rate is $(1-\tilde{\omega})4\pi B(\tau)$.

For conservative scattering, $\tilde{\omega}=1$ so (49) implies that
$dZ(\tau)/d\tau=0$. Then the flux $Z(\tau)$ of (45), and therefore the dipole
moment of the intensity, $I_{1}(\tau)$, are independent of $\tau$.

$Z(\tau)=Z(0),\quad\hbox{and}\quad I_{1}(\tau)=I_{1}(0).$ (52)

Examples of the $\tau$-independence of $I_{1}(\tau)$ in clouds with
conservative scattering are shown as the vertical green lines of Figs. 1, 2, 9
and 10.

Multiplying (1) on the left by $(\kern-2.59995pt(1|$ and using (4) and (40)
with (48) we find

$\displaystyle 0$ $\displaystyle=$
$\displaystyle(\kern-2.59995pt(1|\hat{\mu}\frac{d}{d\tau}|I(\tau)\\}+\eta_{1}I_{1}(\tau)$
(53) $\displaystyle=$
$\displaystyle(\kern-2.59995pt(0|\hat{\mu}^{2}\frac{d}{d\tau}|I(\tau)\\}+\eta_{1}I_{1}(\tau)$
$\displaystyle=$ $\displaystyle\frac{dK(\tau)}{d\tau}+\eta_{1}I_{1}(\tau).$

Figure 1:  Moments of the intensity vector $|I(\tau)\\}$ inside a non-
absorbing, Rayleigh-scattering cloud of total optical depth $\tau_{c}=20$,
calculated with (126), versus optical depth $\tau$ above the bottom of the
cloud. The multipole moments are $I_{l}(\tau)=(\kern-2.59995pt(l|I(\tau)\\}$,
for $l=0,1,2$ and the value of the $K$-integral is given by (47). The input
intensity vector $|I^{\rm in}\\}=|\mu_{10})2/\mu_{10}$ has a direction cosine
$\mu_{10}=0.9739$, the maximum possible for $2n=10$ streams. The dipole moment
of the input intensity is $I_{1}^{\rm in}=(\kern-2.59995pt(1|I^{\rm in}\\}=1$.
Because of diffuse reflection, the first moment $I_{1}$ of the intensity
inside the cloud is much less, $I_{1}=0.0773$. Only 7.73% of the input flux is
difffusely transmitted through the cloud, and the remaining 92.27% is
reflected from the bottom. The dashed red line is a measure of how much the
intensity violates the Eddington approximation, $I_{2}=0$. See the text for
more details.

Figure 2:  Like Fig. 1, but for a forward scattering phase function
$p(\mu)=\varpi^{\\{p\\}}(\mu)$ of (WH-134) with $p=5$. According to (WH-138),
the peak forward scattering is $p(1)=p(p+1)=30$, compared to $p(1)=1.5$ for
Rayleigh scattering, in accordance with (WH-132). Because the incident photons
penetrate more deeply into the cloud before their first scattering, compared
to the nearly isotropic Rayleigh scattering of Fig. 1, and since the photons
have a harder time reversing direction by multiple scattering, the diffuse
transmission is substantially greater than for Fig. 1, 32.99% compared to
7.73%. The diffuse reflection is correspondingly less, 67.01% compared to
92.27%.

Since $I_{1}(\tau)$ is independent of $\tau$, in accordance with (52), we can
solve (53) to find that the $K$-integral changes linearly with optical depth
$\tau$.

$K(\tau)=K(0)-\tau\eta_{1}I_{1}(0).$ (54)

Examples of the linear dependence of $K(\tau)$ on $\tau$ for clouds with
conservative scattering are shown as the dashed blue lines of Figs. 1, 2, 9
and 10. Chandresekhar[4] discusses the basic features of conservative
scattering outlined above in his sections §8-10.

###  3.1 Propagation

From now on we will consider radiation in clouds that are cold enough to
neglect thermal emission. Sunlight in Earth’s clouds is an example of such a
situation. Then we can set $B(\tau)=0$ on the right side of (1). In analogy to
(23) we write the direction-secant matrix

$\displaystyle\hat{\varsigma}$ $\displaystyle=$ $\displaystyle\hat{\mu}^{-1}$
(55) $\displaystyle=$
$\displaystyle\sum_{i=1}^{2n}\varsigma_{i}|\mu_{i})(\kern-2.59995pt(\mu_{i}|,\quad\hbox{where}\quad\varsigma_{i}=\frac{1}{\mu_{i}},$

Multiplying (1) on the left by $\hat{\varsigma}$ from (55) and recalling that
the exponentiation rate matrix is

$\hat{\kappa}=\hat{\varsigma}\hat{\eta}$ (56)

we find

$\left(\frac{d}{d\tau}+\hat{\kappa}\right)|I(\tau)\\}=\hat{0}.$ (57)

If we multiply (57) on the left by $(\kern-2.59995pt(\mu_{i}|$, and assume
isotropic scattering with $\hat{p}=2|0)(\kern-2.59995pt(0|$, we find the
fundamental Eq. (7) of Wick’s pioneering paper[2] on $2n$-stream scattering,
which he used to analyze thermal neutron diffusion with isotropic scattering
cross sections.

The multipole-space representations (WH-44) of $\hat{\varsigma}$, and (40) of
$\hat{\eta}$, allow us to write the exponentiation rate matrix of (57) as

$\displaystyle(\kern-2.59995pt(l|\hat{\kappa}|l^{\prime})$ $\displaystyle=$
$\displaystyle(\kern-2.59995pt(l|\hat{\varsigma}|l^{\prime})\eta_{l^{\prime}}$
(65) $\displaystyle=$
$\displaystyle\left[\begin{array}[]{rrrrrrr}\eta_{0}&3\eta_{1}&0&-14\eta_{3}/3&0&88\eta_{5}/15&\cdots\\\
0&0&0&0&0&0&\cdots\\\ 0&0&0&7\eta_{3}/3&0&-44\eta_{5}/15&\cdots\\\
0&0&5\eta_{2}/3&0&0&0&\cdots\\\ 0&0&0&0&0&11\eta_{5}/5&\cdots\\\
0&0&-4\eta_{2}/3&0&9\eta_{4}/5&0&\cdots\\\
\vdots&\vdots&\vdots&\vdots&\vdots&\vdots&\ddots\\\ \end{array}\right].$

From inspection of (65), and for future reference, we note the identity

$(\kern-2.59995pt(0|\hat{\kappa}|l)+2(\kern-2.59995pt(2|\hat{\kappa}|l)=\eta_{0}\delta_{0l}+3\eta_{1}\delta_{1l}.$
(66)

One can use Bonnet’s recursion formula of (WH-18) to write a formal proof of
(66).

In the simple analysis of WH, eigenvectors of the propagation-length matrix,
$\hat{\lambda}=\hat{\kappa}^{-1}$, were used to describe how the intensity
vector $|I(\tau)\\}$ depends on $\tau$. This procedure fails for conservative
scattering because $\eta_{0}=1-\tilde{\omega}\to 0$ as $\tilde{\omega}\to 1$
and therefore all elements of the first column of (65) vanish,
$(\kern-2.59995pt(l|\hat{\kappa}|0)\to 0$, as $\tilde{\omega}\to 1$. Then
$\det\hat{\kappa}\to 0$ and $\hat{\kappa}$ cannot be inverted. Here we show
how to get around this problem.

For further analysis we will assume that $\eta_{0}=0$ and we will find $2n$
independent solutions of (57), which we denote by
$|v_{1}),|v_{2}),|v_{3}),\ldots,|v_{2n})$. All of these $2n\times 1$
propagation basis vectors, $|v_{i})=|v_{i}(\tau)\big{)}$, depend on $\tau$ and
satisfy the equation

$\left(\frac{d}{d\tau}+\hat{\kappa}\right)|v_{i})=\hat{0},$ (67)

Then the intensity vector $|I(\tau)\\}$ can be written as

$|I(\tau)\\}=\sum_{i=1}^{2n}|v_{i})A_{i}$ (68)

The values of the $\tau$-independent amplitudes $A_{i}$ are determined by
boundary conditions. We can represent the bases $|v_{i})$ as $\tau$-dependent,
$2n\times 1$ column vectors with multipole amplitudes
$(\kern-2.59995pt(l|v_{i})$,

$|v_{i})=\left[\begin{array}[]{r}(\kern-2.59995pt(0|v_{i})\\\
(\kern-2.59995pt(1|v_{i})\\\ (\kern-2.59995pt(2|v_{i})\\\ \vdots\\\
(\kern-2.59995pt(2n-1|v_{i})\end{array}\right].$ (69)

One simple solution to (67) which can be found by inspection of (65) is

$|v_{1}\big{)}=3\eta_{1}\tau|0)-|1)=\left[\begin{array}[]{r}3\eta_{1}\tau\\\
-1\\\ 0\\\ \vdots\end{array}\right].$ (70)

For $|v_{1})$ the energy density $u$ of (43) is

$u=\frac{4\pi(\kern-2.59995pt(0|v_{1}\big{)}}{c}=\frac{12\pi\eta_{1}\tau}{c}.$
(71)

The flux, $Z(\tau)$ of (45), is

$Z(\tau)=4\pi(\kern-2.59995pt(1|v_{1}\big{)}=-4\pi.$ (72)

The $K$-integral $K(\tau)$ (47) of $|v_{2n})$ can be evaluated with (20) and
is

$K(\tau)=(\kern-2.59995pt(0|\hat{\mu}^{2}|v_{1})=\frac{(\kern-2.59995pt(0|v_{1})+2(\kern-2.59995pt(2|v_{1})}{3}=\eta_{1}\tau.$
(73)

Thus, $|v_{1})$ represents a $\tau$-independent downward flux,
$Z(\tau)=-4\pi$. The energy density $u(\tau)$ and the $K$-integral $K(\tau)$
have maximum values at the cloud top, where $\tau=\tau_{c}$, and they decrease
linearly to zero at the cloud bottom where $\tau=0$. The basis vector
$|v_{1})$ of (70) is mentioned as §10(80) by Chandrasekhar[4]. Our multipole
coefficient $p_{1}$ is his $\varpi_{1}/3$, since Chandrasekhar does not
include a statistical weight in his multipole expansion of the scattering
phase, §3(33).

An analog of (70), which also satisfies (67), is

$|v_{2n}\big{)}=3\eta_{1}(\tau_{c}-\tau)|0)+|1)=\left[\begin{array}[]{r}3\eta_{1}(\tau_{c}-\tau)\\\
1\\\ 0\\\ \vdots\end{array}\right].$ (74)

For $|v_{2n})$ the energy density $u(\tau)$ of (43) is

$u(\tau)=\frac{4\pi(\kern-2.59995pt(0|v_{2n}\big{)}}{c}=\frac{12\pi\eta_{1}(\tau_{c}-\tau)}{c}.$
(75)

The flux, $Z(\tau)$ of (45) is

$Z(\tau)=4\pi(\kern-2.59995pt(1|v_{2n}\big{)}=4\pi.$ (76)

The $K$-integral $K(\tau)$ (47) of $|v_{2n})$ can be evaluated with (20) and
is

$K(\tau)=(\kern-2.59995pt(0|\hat{\mu}^{2}|v_{2n})=\frac{(\kern-2.59995pt(0|v_{2n})+2(\kern-2.59995pt(2|v_{2n})}{3}=\eta_{1}(\tau_{c}-\tau).$
(77)

Thus, $|v_{2n})$ represents a $\tau$-independent upward flux, $Z(\tau)=4\pi$.
The energy density $u(\tau)$ and the $K$-integral $K(\tau)$ have maximum
values at the cloud bottom, where $\tau=0$, and they decrease linearly to zero
at the cloud top where $\tau=\tau_{c}$. Since the bases $|v_{1})$ and
$|v_{2n})$ contain only isotropic S waves with $l=0$ and nearly isotropic P
waves with $l=1$, we will call them the two quasi-isotropic bases, to
distinguish them from the remaining basis vectors
$|v_{2}),|v_{3}),\ldots,|v_{2n-1})$, which we will call directional bases.
Examples of quasi-isotropic and directional bases are shown in Fig. 3 and 5.

Numerical studies suggest that the directional bases $|v_{i})$, with
$i=2,3,\ldots,2n-1$, can be written as the sum of an isotropic part
$|v_{i}^{I})$ and a directional part $|v_{i}^{D})$,

$|v_{i})=|v_{i}^{I})+|v_{i}^{D}).$ (78)

The isotropic part is pure S-wave,

$|v_{i}^{I})=|0)(\kern-2.59995pt(0|v_{i}).$ (79)

The directional part is waves of multipolarity $l\geq 2$,

$|v_{i}^{D})=\sum_{l=2}^{2n-1}|l)(\kern-2.59995pt(l|v_{i})$ (80)

As we will prove below, the directional bases have no P-wave part,

$(\kern-2.59995pt(1|v_{i})=0.$ (81)

Therefore the directional bases $|v_{i})$, for $i=2,3,\ldots,2n-1$, carry no
flux, $Z=4\pi(\kern-2.59995pt(1|v_{i})=0$.

Then (67) becomes

$\left(\frac{d}{d\tau}+\left[\begin{array}[]{rr}\hat{\kappa}^{I}&\hat{\kappa}^{X}\\\
\hat{0}&\hat{\kappa}^{D}\end{array}\right]\right)\left[\begin{array}[]{r}|v_{i}^{I})\\\
|v_{i}^{D})\end{array}\right]=\hat{0}.$ (82)

From inspection of (65) we see that

$\hat{\kappa}^{I}=\left[\begin{array}[]{rr}0&3\eta_{1}\\\
0&0\end{array}\right]=3\eta_{1}|0)(\kern-2.59995pt(1|.$ (83)

The $2\times(2n-2)$ matrix that couples directional components $|v_{i}^{D})$
of $|v_{i})$ to quasi-isotropic components is

$\hat{\kappa}^{X}=\left[\begin{array}[]{ccccc}0&-14\eta_{3}/3&0&88\eta_{5}/15&\cdots\\\
0&0&0&0&\cdots\end{array}\right]=\sum_{l=2}^{2n-1}|0)(\kern-2.59995pt(0|\hat{\kappa}|l)(\kern-2.59995pt(l|.$
(84)

The $(2n-2)\times(2n-2)$ matrix that couples directional components of
$|v_{i})$ to each other is

$\hat{\kappa}^{D}=\left[\begin{array}[]{rrrrr}0&7\eta_{3}/3&0&-44\eta_{5}/15&\cdots\\\
5\eta_{2}/3&0&0&0&\cdots\\\ 0&0&0&11\eta_{5}/5&\cdots\\\
-4\eta_{2}/3&0&9\eta_{4}/5&0&\cdots\\\ \vdots&\vdots&\vdots&\vdots&\ddots\\\
\end{array}\right]=\sum_{l\geq 2,l^{\prime}\geq
2}|l)(\kern-2.59995pt(l|\kappa|l^{\prime})(\kern-2.59995pt(l^{\prime}|.$ (85)

The lowest row of the block-matrix equation (82) can be written as the
$(2n-2)\times(2n-2)$ matrix equation,

$\left(\frac{d}{d\tau}+\kappa^{D}\right)|v_{i}^{D})=\hat{0}.$ (86)

The matrix $\kappa^{D}$ will have $(2n-2)$ right eigenvectors $|\lambda_{i})$
and non-zero eigenvalues $\kappa_{i}$, defined, aside from normalization by

$\hat{\kappa}^{D}|\lambda_{i})=\kappa_{i}|\lambda_{i}).$ (87)

It will be convenient to label the eigenvectors $|\lambda_{i})$ with the
inverses $\lambda_{i}$ of the eigenvalues $\kappa_{i}$,

$\lambda_{i}=\frac{1}{\kappa_{i}}$ (88)

We will call the inverse eigenvalues penetration lengths, and we will assume
that they are ordered such that

$\lambda_{2}<\lambda_{3}<\cdots<\lambda_{2n-1}.$ (89)

Since $\kappa^{D}$ is odd under reflection,
$\hat{r}\kappa^{D}\hat{r}=-\kappa^{D}$, where the reflection operator
$\hat{r}$ of (WH-116) can be written as

$\hat{r}=\sum_{l=0}^{2n-1}(-1)^{l}|l)(\kern-2.59995pt(l|,$ (90)

the eigenvectors can be chosen to be reflection conjugates of each other

$\hat{r}|\lambda_{i})=|\lambda_{r(i)}),\quad\hbox{where}\quad r(i)=2n-i+1.$
(91)

The corresponding eigenvalues must be equal and opposite,

$\lambda_{i}=-\lambda_{r(i)}.$ (92)

We guess that the isotropic and directional parts of the bases $|v_{i})$ are

$\left[\begin{array}[]{r}|v_{i}^{I})\\\
|v_{i}^{D})\end{array}\right]=\left[\begin{array}[]{c}-2|0)(\kern-2.59995pt(2|\lambda_{i})\\\
|\lambda_{i})\end{array}\right]e^{-\kappa_{i}(\tau-\tau_{i})},$ (93)

or

$(\kern-2.59995pt(l|v_{i})=\left[\begin{array}[]{r}-2(\kern-2.59995pt(2|\lambda_{i})\\\
0\\\ (\kern-2.59995pt(2|\lambda_{i})\\\ (\kern-2.59995pt(3|\lambda_{i})\\\
\vdots\\\
(\kern-2.59995pt(2n-1|\lambda_{i})\end{array}\right]e^{-\kappa_{i}(\tau-\tau_{i})}.$
(94)

The amplitude of the isotropic, S-wave part of $|v_{i})$ is
$(\kern-2.59995pt(0|v_{i})=-2(\kern-2.59995pt(2|\lambda_{i})$.

In (93) and (94) the reference optical depths $\tau_{i}$ are

$\tau_{i}=\left\\{\begin{array}[]{rl}\tau_{c},&\mbox{if $i=2,3,4,\ldots,n$,
}\\\ 0,&\mbox{if $i=n+1,n+2,n+3,\ldots,2n-1$.}\end{array}\right.$ (95)

To prove that the bases $|v_{i})$ of (94) are solutions of (67) we first note
that differentiating (94) gives

$\frac{d}{d\tau}|v_{i})=-\kappa_{i}|v_{i}).$ (96)

Noting from (83) and (94) that $\hat{\kappa}^{I}|v_{i}^{I})=0$ we see that

$\hat{\kappa}|v_{i})=\left[\begin{array}[]{rr}\hat{\kappa}^{I}&\hat{\kappa}^{X}\\\
\hat{0}&\hat{\kappa}^{D}\end{array}\right]\left[\begin{array}[]{c}|v_{i}^{I})\\\
|v_{i}^{D})\end{array}\right]=\left[\begin{array}[]{c}\hat{\kappa}^{X}|v_{i}^{D})\\\
\hat{\kappa}^{D}|v_{i}^{D})\end{array}\right].$ (97)

The bottom element on the right of (97) is

$\hat{\kappa}^{D}|v_{i}^{D})=\hat{\kappa}^{D}|\lambda_{i})e^{-\kappa_{i}(\tau-\tau_{i})}=\kappa_{i}|\lambda_{i})e^{-\kappa_{i}(\tau-\tau_{i})}=\kappa_{i}|v_{i}^{D}).$
(98)

The top element on the right of (97) is

$\displaystyle\hat{\kappa}^{X}|v_{i}^{D})$ $\displaystyle=$
$\displaystyle\sum_{l=2}^{2n-1}|0)(\kern-2.59995pt(0|\hat{\kappa}|l)(\kern-2.59995pt(l|\lambda_{i})e^{-\kappa_{i}(\tau-\tau_{i})}$
(99) $\displaystyle=$
$\displaystyle-2\sum_{l=2}^{2n-1}|0)(\kern-2.59995pt(2|\hat{\kappa}^{D}|l)(\kern-2.59995pt(l|\lambda_{i})e^{-\kappa_{i}(\tau-\tau_{i})}$
$\displaystyle=$
$\displaystyle-2|0)(\kern-2.59995pt(2|\hat{\kappa}^{D}|\lambda_{i})e^{-\kappa_{i}(\tau-\tau_{i})}$
$\displaystyle=$
$\displaystyle-2\kappa_{i}|0)(\kern-2.59995pt(2|\lambda_{i})e^{-\kappa_{i}(\tau-\tau_{i})}$
$\displaystyle=$ $\displaystyle\kappa_{i}|v^{I}_{i}).$

The top line of (99) comes from (84). To get the second line we noted that for
$l\geq 2$ we can use (66) to write
$(\kern-2.59995pt(0|\hat{\kappa}|l)=-2(\kern-2.59995pt(2|\hat{\kappa}|l)=-2(\kern-2.59995pt(2|\hat{\kappa}^{D}|l)$.
The sum on $|l)(\kern-2.59995pt(l|$ then becomes a unit operator in
directional space, which leads to the third line. The fourth line follows from
(87). The final line follows from (93). Using (99) and (98) in (97) we find

$\hat{\kappa}|v_{i})=\kappa_{i}|v_{i}).$ (100)

Adding (100) to (96) we find (67).

The eigenvalue equations (87) do not determine how the eigenvectors
$|\lambda_{i})$ are normalized. It will be convenient to let the top element
have the value 1/2,

$-2(\kern-2.59995pt(2|\lambda_{i})=\frac{1}{2}.$ (101)

Then we can write the directional basis vectors $|v_{i})$ of (94), with
$i=2,3,\ldots,2n-1$, as

$|v_{i})=\left[\begin{array}[]{r}(\kern-2.59995pt(0|v_{i})\\\
(\kern-2.59995pt(1|v_{i})\\\ (\kern-2.59995pt(2|v_{i})\\\
(\kern-2.59995pt(3|v_{i})\\\ \vdots\\\
(\kern-2.59995pt(2n-1|v_{i})\end{array}\right]=\left[\begin{array}[]{r}0.5\\\
0\\\ -0.25\\\ (\kern-2.59995pt(3|\lambda_{i})\\\ \vdots\\\
(\kern-2.59995pt(2n-1|\lambda_{i})\end{array}\right]e^{-\kappa_{i}(\tau-\tau_{i})}.$
(102)

For the directional basis vectors $|v_{i})$ of (102), the energy density
$u(\tau)$, given by (43), decays exponentially with distance into the cloud
from the reference optical depth $\tau_{i}$,

$u=\frac{4\pi}{c}(\kern-2.59995pt(0|v_{i})=\frac{2\pi}{c}e^{-\kappa_{i}(\tau-\tau_{i})}.$
(103)

The flux, $Z(\tau)$, given by (45), is

$Z=4\pi(\kern-2.59995pt(1|v_{i})=0.$ (104)

The $K$-integral $K(\tau)$ given by (47) is also zero for the directional
bases,

$K(\tau)=\frac{(\kern-2.59995pt(0|v_{2n})+2(\kern-2.59995pt(0|v_{2n})}{3}=0.$
(105)

Calculations of radiation transfer in stars sometimes use the Eddington
approximation,

$(\kern-2.59995pt(l|I(\tau)\\}=I_{l}(\tau)=0,\quad\hbox{for}\quad l\geq 2.$
(106)

An example of how valid the Eddington criterion (106) is for a realistic cloud
can be seen in Fig. 1. For optical depths from about $\tau=2$ above the bottom
of the cloud to about 2 optical depths below the top of the cloud at
$\tau=20$, the Eddington approximation is very well satisfied, and
$2I_{2}(\tau)/I_{0}(\tau)\approx 0$. But the figure shows that the Eddington
approximation is not good near the bottom of the cloud, where the highly
directional input radiation is being isotropized by multiple scattering, or
near the top of the cloud where isotropic radiation from inside the cloud is
being transformed into purely upward radiation from the top surface. Near the
bottom and top of the cloud the intensity $|I(\tau)\\}$ contains large
fractions of directional bases (102), which do not satisfy the Eddington
criterion (106), since they are normalized such that
$(\kern-2.59995pt(2|v_{i})=-1/4$. But the directional bases decay
exponentially to zero with increasing distances from the top and bottom of the
cloud. Eddington’s criterion is very well satisfied near the center of an
optically thick cloud, where the radiation is mostly due to the quasi-
isotropic bases $|v_{2n})$ and $|v_{1})$, of (70) and (74), which identically
satisfy the Eddington criterion, $I_{l}(\tau)=0$ for $l\geq 2$.

Fig. 1 also shows the monopole moment, $I_{0}(\tau)$, or average intensity.
Isotropization of near-normal-incidence radiation causes the average intensity
to maximize just above the bottom of the cloud. Qualitatively similar
increases of nearly isotropic neutron thermal fluxes are observed because of
slowing down of fast neutrons in efficient neutron reflectors. [6]

The simplest model with both quasi-isotropic and directional basis vectors has
$n=2$ stream pairs. Then the basis vectors are

$(\kern-2.59995pt(l|v_{i})=\left[\begin{array}[]{r r r
r}3\eta_{1}\tau&0.5\,e^{-\kappa_{2}(\tau-\tau_{c})}&0.5\,e^{-\kappa_{3}\tau}&3\eta_{1}(\tau_{c}-\tau)\\\
-1&0&0&1\\\
0&-0.25\,e^{-\kappa_{2}(\tau-\tau_{c})}&-0.25\,e^{-\kappa_{3}\tau}&0\\\
0&(\kern-2.59995pt(3|\lambda_{2})\,e^{-\kappa_{2}(\tau-\tau_{c})}&(\kern-2.59995pt(3|\lambda_{3})\,e^{-\kappa_{3}\tau}&0\end{array}\right]$
(107)

In accordance with the symmetry (92),

$-\kappa_{2}=\kappa_{3}=\frac{\sqrt{35\eta_{2}\eta_{3}}}{3}.$ (108)

In accordance with the symmetry (91),

$(\kern-2.59995pt(3|\lambda_{2})=-(\kern-2.59995pt(3|\lambda_{3})=\frac{1}{4}\sqrt{\frac{5\eta_{2}}{7\eta_{3}}}.$
(109)

Figure (3) shows representative projections

$\displaystyle\langle\mu|v_{i}(\tau)\big{)}$ $\displaystyle=$
$\displaystyle\sum_{l=0}^{2n-1}\langle\mu|l)(\kern-2.59995pt(l|v_{i}(\tau)\big{)}$
(110) $\displaystyle=$
$\displaystyle\sum_{l=0}^{2n-1}P_{l}(\mu)(2l+1)(\kern-2.59995pt(l|v_{i}(\tau)\big{)}.$

of the basis vectors $|v_{i}(\tau)\big{)}$ of (69) onto the continuous
$\mu$-space basis $\langle\mu|$, with $\langle\mu|l)$ given by (6). As can be
seen in Fig. 3, for $|\lambda_{i}|\leq 1$, the directional bases
$\langle\mu|v_{i})$ have their largest values for $\mu\approx\lambda_{i}$. A
basis vector with large amplitudes for $\mu\approx 1$ represents nearly
vertical radiation which penetrates the cloud better than radiation from a
basis vector representing nearly horizontal radiation, for which the largest
values of $\langle\mu|v_{i})$ occur for $\mu\approx 0$. Horizontal radiation
must traverse a longer slant path through the cloud than vertical radiation,
and is attenuated more. So basis vectors $|v_{i})$ representing the most
horizontal radiation have the shortest penetration lengths $\lambda_{i}$.

Figure 3: Representative projections $\langle\mu|v_{i}(\tau))$ of (110) at the
cloud bottom where $\tau=0$. The cloud thickness is $\tau_{c}=20$. Rayleigh
scattering is modeled with $2n=32$ streams. The value of the quasi-isotropic
basis vector $\langle\mu|v_{32}(0)\big{)}$ of (74), shown in the left panel,
varies linearly with the direction cosines $\mu$. The middle panel shows a
nearly vertical directional basis vector $\langle\mu|v_{31}(0)\big{)}$ with a
penetration length $\lambda_{i}=0.981$. The right panel shows a nearly
horizontal directional basis vector $\langle\mu|v_{17}(0)\big{)}$ with a
penetration length $\lambda_{17}=0.051$, indicated by the horizontal green
line. See the text for more details.

Figure 4: The $2n-2=30$ penetration lengths $\lambda_{i}$ of (88) for a
conservative scattering model with $2n=32$ streams, and with the phase
functions $p(\mu)=\varpi^{\\{p\\}}(\mu)$ of (38). These phase functions give
the maximum forward scattering possible for functions constructed from the
first $2p$ Legendre polynomials. The possible values of $p$ are
$1,2,3,\ldots,16$. The values of $\lambda_{i}$ labeled with $p=0$ are for
purely isotropic scattering, with $p(\mu)=1$.

For phase functions $\varpi^{\\{p\\}}(\mu)$ of (38) for $p\gg 1$ and strong
forward scattering, the largest penetration lengths $|\lambda_{i}|$ of (38)
can exceed 1, as can be seen in Fig. 4. The directional basis vectors
$|v_{i})$, corresponding to $\lambda_{i}>1$, are somewhat delocalized, but are
largest near $\mu=1$. The basis vectors corresponding to $\lambda_{i}<-1$ are
largest near $\mu=-1$. Figure 5 shows representative projections for strongly
forward-peaked phase functions.

Figure 5:  Like Fig. 3 but for a strongly peaked forward-scattering phase
function $p(\mu)=\varpi^{\\{16\\}}(\mu)$, constructed from the first $2p=32$
Legendre polynomials and with a forward scattering phase of $p(1)=p(p+1)=272$.
The cloud thickness is $\tau_{c}=20$, and there are $2n=32$ streams. The value
of the quasi-isotropic basis vector $\langle\mu|v_{32})$ of (74), shown in the
left panel, varies linearly with $\mu$. The middle panel shows a nearly
vertical directional basis vector $\langle\mu|v_{31}(0)\big{)}$ with a
penetration length $\lambda_{i}=5.926$. The right panel shows a nearly
horizontal directional basis vector $\langle\mu|v_{17}(0)\big{)}$ with a
penetration length $\lambda_{17}=0.085$, indicated by the horizontal green
line. See the text for more details.

##  4 The scattering matrix

For a cloud with a bottom at $\tau=0$ and a top at $\tau=\tau_{c}$, the
incoming and outgoing intensity vectors are

$|I^{\\{\rm in\\}}\\}=\mathcal{M}_{\bf u}|I(0)\\}+\mathcal{M}_{\bf
d}|I(\tau_{c})\\},$ (111)

and

$|I^{\\{\rm out\\}}\\}=\mathcal{M}_{\bf d}|I(0)\\}+\mathcal{M}_{\bf
u}|I(\tau_{c})\\}.$ (112)

as discussed in (WH-174) and (WH-175). The projection matrices
$\mathcal{M}_{\bf d}$ and $\mathcal{M}_{\bf u}$ for downward streams with
indices $j\leq n$ and upward streams with indices $k>n$, respectively, are
given by (WH-88) as

$\mathcal{M}_{\bf
d}=\sum_{j=1}^{n}|\mu_{j})(\kern-2.59995pt(\mu_{j}|,\quad\hbox{and}\quad\mathcal{M}_{\bf
u}=\sum_{k=n+1}^{2n}|\mu_{k})(\kern-2.59995pt(\mu_{k}|.$ (113)

The projection matrices have the simple algebra

$\displaystyle\mathcal{M}_{\bf d}+\mathcal{M}_{\bf u}$ $\displaystyle=$
$\displaystyle\hat{1},$ (114) $\displaystyle\mathcal{M}^{2}_{\bf
d}=\mathcal{M}_{\bf d}\quad$ and $\displaystyle\quad\mathcal{M}^{2}_{\bf
u}=\mathcal{M}_{\bf u},$ (115) $\displaystyle\mathcal{M}_{\bf
d}\mathcal{M}_{\bf u}=\hat{0}\quad$ and $\displaystyle\quad\mathcal{M}_{\bf
u}\mathcal{M}_{\bf d}=\hat{0},$ (116)

For future reference, we use (114) – (116) to show that the inverses of (111)
and (112) are

$|I(0)\\}=\mathcal{M}_{\bf u}|I^{\\{\rm in\\}}\\}+\mathcal{M}_{\bf
d}|I^{\\{\rm out\\}}\\}.$ (117) $|I(\tau_{c})\\}=\mathcal{M}_{\bf d}|I^{\\{\rm
in\\}}\\}+\mathcal{M}_{\bf u}|I^{\\{\rm out\\}}\\}.$ (118)

Multiplying (112) on the left by $(\kern-2.59995pt(\mu_{i}|$, and using (68)
we see that the intensity of the $i$th outgoing stream is

$\displaystyle(\kern-2.59995pt(\mu_{i}|I^{\\{\rm out\\}}\\}$ $\displaystyle=$
$\displaystyle\sum_{i^{\prime}}\left[(\kern-2.59995pt(\mu_{i}|\mathcal{M}_{\bf
d}|v_{i^{\prime}}(0)\big{)}+(\kern-2.59995pt(\mu_{i}|\mathcal{M}_{\bf
u}|v_{i^{\prime}}(\tau_{c})\big{)}\right]A_{i^{\prime}}$ (119)
$\displaystyle=$
$\displaystyle\sum_{i^{\prime}}\mathcal{O}_{ii^{\prime}}A_{i^{\prime}}$

The elements of the outgoing matrix , with $j\leq n$ and $k>n$, are

$\displaystyle\mathcal{O}_{ji^{\prime}}$ $\displaystyle=$
$\displaystyle(\kern-2.59995pt(\mu_{j}|v_{i^{\prime}}(0)\big{)}=\sum_{l}(\kern-2.59995pt(\mu_{j}|l)(\kern-2.59995pt(l|v_{i^{\prime}}(0)\big{)}$
$\displaystyle\mathcal{O}_{ki^{\prime}}$ $\displaystyle=$
$\displaystyle(\kern-2.59995pt(\mu_{k}|v_{i^{\prime}}(\tau_{c})\big{)}=\sum_{l}(\kern-2.59995pt(\mu_{k}|l)(\kern-2.59995pt(l|v_{i^{\prime}}(\tau_{c})\big{)}\
$ (120)

Figure 6: Intensities inside the Rayleigh-scattering cloud of Fig. 1. The ray
lengths are proportional to the unweighted intensities
$I(\mu_{i},\tau)=w_{i}^{-1}(\kern-2.59995pt(\mu_{i}|I(\tau)\\}$, calculated
with (126), at the optical depth $\tau$. At the bottom of the cloud, with
$\tau=0$, there is an upward incident stream $|I^{\\{\rm
in\\}}=|\mu_{10})2/\mu_{10}$, indicated by dotted rays, and diffusely
reflected downward streams, indicated by the continuous rays. Within a few
optical depths above the bottom, the intensity isotropizes due to multiple
scattering, and there is little trace of the directional incident stream. The
red color of the outgoing rays at the top of the cloud, with
$\tau=\tau_{c}=20$, indicates that they have been magnified by a factor of 10
for clarity. Only upward streams exist from the top of the cloud since we
assumed only upward incident light. If the transmission at the top of the
cloud or the reflection at the bottom were Lambertian, the rays would have the
same lengths. But in both cases there is substantial “limb darkening.” The
nearly vertical rays are stronger than the nearly horizontal rays.

Figure 7: Unweighted intensities
$I(\mu_{i},\tau)=w_{i}^{-1}(\kern-2.59995pt(\mu_{i}|I(\tau)\\}$ like those of
Fig. 6 but for a forward scattering phase function
$p(\mu)=\varpi^{\\{p\\}}(\mu)$ of (WH-134) with $p=5$. According to (WH-138),
the peak forward scattering is $p(1)=p(p+1)=30$. Because the incident photons
penetrate more deeply into the cloud before their first scattering, compared
to the nearly isotropic Rayleigh scattering of Fig. 6, and because the photons
have a harder time reversing direction by multiple scattering, the diffuse
reflection is substantially less, and the diffuse transmission is
substantially greater than for Fig. 6, as discussed in the caption of Fig. 2.
Even at an optical depth $\tau=5$ above the bottom of the cloud, the intensity
has not fully isotropized. Both the diffuse transmission and reflection are
strongly limb darkened.

Figure 8: Unweighted intensities
$I(\mu_{i},\tau)=w_{i}^{-1}(\kern-2.59995pt(\mu_{i}|I(\tau)\\}$ like those of
Fig. 7 but with the most horizontal incident radiation possible for $2n=10$
streams. At the bottom of the cloud, with $\tau=0$, there is a nearly-
horizontal incident stream $|I^{\\{\rm in\\}}\\}=|\mu_{6})2/\mu_{6}$,
indicated with the dotted rays, and diffusely reflected downward streams. The
direction cosine of the input stream is $\mu_{6}=0.1489$. The diffusely
transmitted light at the top of the cloud is limb darkened, with more intense
vertical than horizontal radiation. But at the bottom, the diffuse reflection
is strongly limb brightened since the incident light is so nearly horizontal,
and there is strong forward scattering. The transmission of the cloud, 13.41%,
is much less than for a nearly vertical incident beam like that of Fig. 7 or
Fig. 2, where the transmission is 32.99%.

Figure 9:  Like Fig. 2 but with the most horizontal incident radiation
possible for $2n=10$ streams, as in Fig. 8. The input intensity vector is
$|I^{\\{\rm in\\}}\\}=|\mu_{6})2/\mu_{6}$, the direction cosine of the input
stream is $\mu_{6}=0.1489$, and the input dipole moment of the intensity
vector is $I^{\\{\rm in\\}}_{1}=(\kern-2.59995pt(1|I^{\\{\rm in\\}}\\}=1$. The
energy density, $u=4\pi I_{0}/c$, peaks at the bottom surface of the cloud,
rather than several optical depths inside, as for vertical incident radiation.
The diffuse transmission has been reduced from the 32.99% of Fig. 2 to 13.41%.

Figure 10: Intensity moments inside a cloud with a forward scattering phase
function $p(\mu)=\varpi^{\\{p\\}}(\mu)$ of (38), constructed from the first
$2p=10$ Legendre polynomials. There is half-isotropic input radiation at the
bottom, $|I^{\rm in}\\}=\mathcal{M}_{\bf
u}|0)/(\kern-2.59995pt(1|\mathcal{M}_{\bf u}|0)$, normalized such that
$I_{1}^{\\{\rm in\\}}=(\kern-2.59995pt(1|I^{\rm in}\\}=1$. The transmission,
$I_{1}(\tau)$, is 26.14%, intermediate between the transmission, 32.99% of the
most vertical possible stream of Fig. 2, and the transmission 13.41% of the
most horizontal possible beam of Fig. 9. The radiation stream directions
inside the cloud are shown by the red rays of Fig. 13.

Figure 11: The diffuse reflectances $R_{\bf u}$ of (147) for upward incident
streams, $|I^{\\{\rm in\\}}\\}=|\mu_{i})$, with $i=n+1,n+2,\ldots,2n$,
incident on the bottom of the cloud. Results are shown for three cloud
thicknesses, $\tau_{c}=20,5$ and $1$. The Rayleigh-scattering phase function
$p(\mu)$ of (36) is assumed. Results for $2n=10$ streams (large circles) and
$2n=32$ streams (small black squares) hardly differ. The reflectance increases
as the incident beams become more horizontal and as the cloud thickness
increases. For conservative scattering the transmittance is $T_{\bf
u}=1-R_{\bf u}$, in accordance with (152).

Figure 12: Like Fig. 11 but for a forward-scattering phase function
$p(\mu)=\varpi^{\\{p\\}}(\mu)$ of (38), with $p=5$. The forward-scattering
amplitude, $p(1)=p(p+1)=30$, is twenty times larger than that for Rayleigh
scattering, $p(1)=1.5$. Because of the forward scattering, the stream
reflectances $R_{\bf u}$ are less and transmissivities, $T_{\bf u}=1-R_{\bf
u}$, are greater than for Fig. 11.

Figure 13: A cloud with the maximum forward scattering phase function
$p(\mu)=\varpi^{\\{p\\}}(\mu)$ of (38) constructed from $p=5$ pairs of
Legendre polynomials. The incoming intensity vector is isotropic, $|I^{\\{\rm
in\\}}\\}=|0)$. The optical thickness is $\tau_{c}=20$ and the radiation is
modeled with $2n=10$ streams. Unweighted intensities $I^{\\{\bf
u\\}}(\mu_{i},\tau)=w_{i}^{-1}(\kern-2.59995pt(\mu_{i}|I^{\\{\bf
u\\}}(\tau)\\}$, originating as upward streams onto the bottom of the cloud
are shown as red rays and are given by (LABEL:iir14). Unweighted intensities
$I^{\\{\bf d\\}}(\mu_{i},\tau)=w_{i}^{-1}(\kern-2.59995pt(\mu_{i}|I^{\\{\bf
d\\}}(\tau)\\}$ originating as downward streams onto the top of the cloud are
shown as blue rays and are given by (LABEL:iir16). The top and bottom input
rays are dotted. 73.86% of the radiation incident on the bottom is diffusely
reflected with 26.24% reaching the top, and vice versa for radiation incident
on the top. The moments of the red rays are shown in Fig. 10. The total
radiation remains isotropic and of equal intensity everywhere inside and
outside the cloud, as discussed in the text.

Figure 14: Like Fig. 13 but for a thinner cloud, with optical thickness
$\tau_{c}=1$. A smaller fraction, 15.73%, of the “red” radiation incident onto
the bottom of the cloud is reflected, with a larger fraction, 84.27%,
transmitted at the top, and vice versa for the “blue” radiation incident onto
the top. The more vertical streams are less attenuated than the more
horizontal streams, which leads to substantial limb darkening of the red
emergent radiation at the top, and the blue emergent radiation at the bottom.

Multiplying (111) on the left by $(\kern-2.59995pt(\mu_{i}|$, we see that the
intensity of the $i$th incoming stream is

$\displaystyle(\kern-2.59995pt(\mu_{i}|I^{\\{\rm in\\}}\\}$ $\displaystyle=$
$\displaystyle\sum_{i^{\prime}}\left[(\kern-2.59995pt(\mu_{i}|\mathcal{M}_{\bf
u}|v_{i^{\prime}}(0)\big{)}+(\kern-2.59995pt(\mu_{i}|\mathcal{M}_{\bf
d}|v_{i^{\prime}}(\tau_{c})\big{)}\right]A_{i^{\prime}}$ (121)
$\displaystyle=$
$\displaystyle\sum_{i^{\prime}}\mathcal{I}_{ii^{\prime}}A_{i^{\prime}}$

The elements of the incoming matrix, with $j\leq n$ and $k>n$, are

$\displaystyle\mathcal{I}_{ji^{\prime}}$ $\displaystyle=$
$\displaystyle(\kern-2.59995pt(\mu_{j}|v_{i^{\prime}}(\tau_{c})\big{)}=\sum_{l}(\kern-2.59995pt(\mu_{j}|l)(\kern-2.59995pt(l|v_{i^{\prime}}(\tau_{c})\big{)}$
$\displaystyle\mathcal{I}_{ki^{\prime}}$ $\displaystyle=$
$\displaystyle(\kern-2.59995pt(\mu_{k}|v_{i^{\prime}}(0)\big{)}=\sum_{l}(\kern-2.59995pt(\mu_{k}|l)(\kern-2.59995pt(l|v_{i^{\prime}}(0)\big{)}\
$ (122)

Inverting (121) we find

$A_{i^{\prime}}=\sum_{i^{\prime\prime}}\mathcal{I}^{-1}_{i^{\prime}i^{\prime\prime}}(\kern-2.59995pt(\mu_{i^{\prime\prime}}|I^{\\{\rm
in\\}}\\}$ (123)

Substituting (123) into (119) we find

$(\kern-2.59995pt(\mu_{i}|I^{\\{\rm
out\\}}\\}=\sum_{i^{\prime\prime}}\mathcal{S}_{ii^{\prime\prime}}(\kern-2.59995pt(\mu_{i^{\prime\prime}}|I^{\\{\rm
in\\}}\\},\quad\hbox{or}\quad|I^{\\{\rm out\\}}\\}=\mathcal{S}|I^{\\{\rm
in\\}}\\}.$ (124)

The elements of the scattering matrix are

$\mathcal{S}_{ii^{\prime\prime}}=\sum_{i^{\prime}}\mathcal{O}_{ii^{\prime}}\mathcal{I}^{-1}_{i^{\prime}i^{\prime\prime}},\quad\hbox{or}\quad\mathcal{S}=\mathcal{O}\mathcal{I}^{-1}.$
(125)

For radiation transfer with finite absorption, $\tilde{\omega}<1$, which was
considered in WH, there is no need to introduce fundamentally different
formulas for the quasi-isotropic bases $|v_{1})$ and $|v_{2n})$, which both
have a linear dependence on optical depth, and the directional bases,
$|v_{2})$, $|v_{3}),\ldots|v_{2n-1})$, which evolve exponentially with $\tau$.
With finite absorption, even the quasi-isotropic modes $|\lambda_{1})$ and
$|\lambda_{2n})$ evolve exponentially with $\tau$, albeit more and more slowly
as $\tilde{\omega}\to 1$. As shown by (WH-173), the left bases,
$(\kern-2.59995pt(\lambda_{1}|$ and $(\kern-2.59995pt(\lambda_{2n}|$ diverge
as $(1-\tilde{\omega})^{-1/2}$ for the conservative scattering limit,
$\tilde{\omega}\to 1$.

###  4.1 Intracloud intensities

Using (123) with (68), we write the intensity vector for optical depths $\tau$
inside a cloud as

$\displaystyle|I(\tau)\\}$ $\displaystyle=$
$\displaystyle\sum_{i}^{2n}|v_{i^{\prime}}(\tau)\big{)}A_{i^{\prime}}$ (126)
$\displaystyle=$
$\displaystyle\sum_{i^{\prime}i^{\prime\prime}}|v_{i^{\prime}}(\tau)\big{)}\mathcal{I}^{-1}_{i^{\prime}i^{\prime\prime}}(\kern-2.59995pt(\mu_{i^{\prime\prime}}|I^{\\{\rm
in\\}}\\}.$

In (126) and (127) the summation indices $i^{\prime}$ and $i^{\prime\prime}$
take on the values $1,2,3,\ldots 2n$. The projections of (126) onto the stream
bases are

$\displaystyle(\kern-2.59995pt(\mu_{i}|I(\tau)\\}=\sum_{l=0}^{2n-1}(\kern-2.59995pt(\mu_{i}|l)(\kern-2.59995pt(l|I(\tau)\\}$
(127)

Examples of (127) are shown in Figs. 6, 7 and 8 where the lengths of the rays
are equal to the unweighted values
$w_{i}^{-1}(\kern-2.59995pt(\mu_{i}|I(\tau)\\}$ of the stream intensities.
These correspond to the intensity moments of Figs. 1, 2 and 9.

###  4.2 Identities for $\mathcal{S}$ and $\Omega$

The cloud albedo matrix $\Omega$, defined in (WH-216) is related to the
scattering matrix by the similarity transformation

$\Omega=(\hat{\mu}_{\bf u}-\hat{\mu}_{\bf d})\mathcal{S}(\hat{\varsigma}_{\bf
u}-\hat{\varsigma}_{\bf d}).$ (128)

The upward and downward parts of the direction-cosine matrix (23) are

$\hat{\mu}_{\bf u}=\hat{\mu}\mathcal{M}_{\bf u}=\mathcal{M}_{\bf
u}\hat{\mu},\quad\hbox{and}\quad\hat{\mu}_{\bf d}=\hat{\mu}\mathcal{M}_{\bf
d}=\mathcal{M}_{\bf d}\hat{\mu},$ (129)

and the upward and downward parts of the direction-secant matrix (55) are

$\hat{\varsigma}_{\bf u}=\hat{\varsigma}\mathcal{M}_{\bf u}=\mathcal{M}_{\bf
u}\hat{\varsigma},\quad\hbox{and}\quad\hat{\varsigma}_{\bf
d}=\hat{\varsigma}\mathcal{M}_{\bf d}=\mathcal{M}_{\bf d}\hat{\varsigma},$
(130)

According to (52), the dipole moment of the intensity, $I_{1}(\tau)$ is
independent of the optical depth $\tau$ above the bottom of a cloud, so we can
write

$\displaystyle 0$ $\displaystyle=$ $\displaystyle I_{1}(\tau_{c})-I_{1}(0)$
(131) $\displaystyle=$ $\displaystyle(\kern-2.59995pt(1|\left(\mathcal{M}_{\bf
d}|I^{\\{\rm in\\}}\\}+\mathcal{M}_{\bf u}|I^{\\{\rm
out\\}}\\}-\mathcal{M}_{\bf u}|I^{\\{\rm in\\}}\\}-\mathcal{M}_{\bf
d}|I^{\\{\rm out\\}}\\}\right)$ $\displaystyle=$
$\displaystyle(\kern-2.59995pt(0|\hat{\mu}\left(\mathcal{M}_{\bf
d}-\mathcal{M}_{\bf u}+[\mathcal{M}_{\bf u}-\mathcal{M}_{\bf
d}]\mathcal{S}\right)|I^{\\{\rm in\\}}\\}$ $\displaystyle=$
$\displaystyle(\kern-2.59995pt(0|\left(\hat{\mu}_{\bf d}-\hat{\mu}_{\bf
u}+[\hat{\mu}_{\bf u}-\hat{\mu}_{\bf d}]\mathcal{S}\right)|I^{\\{\rm
in\\}}\\}.$

In the second line of (131) we used (4), (117) and (118). The third line
follows from (19) and (124). Eq. (131) is true for any $|I^{\\{\rm in\\}}\\}$
so it implies that

$(\kern-2.59995pt(0|(\hat{\mu}_{\bf u}-\hat{\mu}_{\bf
d})\mathcal{S}=(\kern-2.59995pt(0|(\hat{\mu}_{\bf u}-\hat{\mu}_{\bf d}).$
(132)

Multiplying both sides of (132) on the right by $(\hat{\varsigma}_{\bf
u}-\hat{\varsigma}_{\bf d})$ and using (128) we find,

$(\kern-2.59995pt(0|\Omega=(\kern-2.59995pt(0|$ (133)

The identity (133) can serve as a useful consistency check for numerical
calculations.

In analogy to (133), for conservative scattering the $S$-matrix satisfies the
identity

$S|0)=|0),$ (134)

which can also serve as a useful check of numerical calculations. A simple
proof of (134) follows from Kirchhoff’s law of (WH-279)

$\mathcal{E}=\hat{1}-\mathcal{S}$ (135)

Here $\mathcal{E}$ is the isothermal emissivity matrix of the cloud.
Kirchoff’s law remains valid as we approach the conservative scattering limit
$\tilde{\omega}\to 1$, when the thermal emission of the cloud must approach
zero,

$\mathcal{E}|B\\}=\mathcal{E}|0)B=(\hat{1}-\mathcal{S})|0)B\to
0,\quad\hbox{as}\quad\tilde{\omega}\to 1.$ (136)

Here $B$ is the Planck intensity of an isothermal cloud, given by (WH-7). The
limit (136) proves (134).

##  5 Transmission, reflection and absorption

For greater clarity, we will consider transmission and reflection by a cloud
with a single scattering albedo $\tilde{\omega}<1$ and therefore some
absorption. In the conservative scattering limit of $\tilde{\omega}\to 1$ the
absorption vanishes.

As shown by (WH-212) and (WH-214) the incoming and outgoing flux vectors
$|Z^{\\{\rm in\\}}\\}$ and $|Z^{\\{\rm out\\}}\\}$ are related to the incoming
and outgoing intensities of (111) and (112) by

$\displaystyle|Z^{\\{\rm in\\}}\\}$ $\displaystyle=$ $\displaystyle
4\pi(\hat{\mu}_{\bf u}-\hat{\mu}_{\bf d})|I^{\\{\rm in\\}}\\},$ (137)
$\displaystyle|Z^{\\{\rm out\\}}\\}$ $\displaystyle=$ $\displaystyle
4\pi(\hat{\mu}_{\bf u}-\hat{\mu}_{\bf d})|I^{\\{\rm out\\}}\\}.$ (138)

The albedo matrix $\Omega$ gives the proportionality of $|Z^{\\{\rm
out\\}}\\}$ to $|Z^{\\{\rm in\\}}\\}$,

$|Z^{\\{\rm out\\}}\\}=\Omega|Z^{\\{\rm in\\}}\\}.$ (139)

We will define the absorption probability of the cloud as

$A=1-\langle\Omega\rangle$ (140)

As shown in (WH-217) the expectation value $\langle\Omega\rangle$ of the cloud
albedo depends on the input flux vector $|Z^{\\{\rm in\\}}\\}$, and is given
by

$\displaystyle\langle\Omega\rangle$ $\displaystyle=$
$\displaystyle\frac{Z^{\\{\rm out\\}}}{Z^{\\{\rm
in\\}}}=\frac{(\kern-2.59995pt(0|\Omega|Z^{\\{\rm
in\\}}\\}}{(\kern-2.59995pt(0|Z^{\\{\rm in\\}}\\}}.$ (141)

Substituting (141) into (140) we find that the absorption probability can be
written as

$A=A_{\bf d}f_{\bf d}+A_{\bf u}f_{\bf u}.$ (142)

The fractions of downward and upward input flux are

$f_{\bf d}=\frac{(\kern-2.59995pt(0|\mathcal{M}_{\bf d}|Z^{\\{\rm
in\\}}\\}}{(\kern-2.59995pt(0|Z^{\\{\rm in\\}}\\}},\quad\hbox{and}\quad f_{\bf
u}=\frac{(\kern-2.59995pt(0|\mathcal{M}_{\bf u}|Z^{\\{\rm
in\\}}\\}}{(\kern-2.59995pt(0|Z^{\\{\rm in\\}}\\}}.$ (143)

The absorption probabilities for the downward and upward fractions of the
input flux are

$A_{\bf d}=\frac{(\kern-2.59995pt(0|(1-\Omega)\mathcal{M}_{\bf d}|Z^{\\{\rm
in\\}}\\}}{(\kern-2.59995pt(0|\mathcal{M}_{\bf d}|Z^{\\{\rm
in\\}}\\}}=1-\frac{(\kern-2.59995pt(0|\Omega\mathcal{M}_{\bf d}|Z^{\\{\rm
in\\}}\\}}{(\kern-2.59995pt(0|\mathcal{M}_{\bf d}|Z^{\\{\rm in\\}}\\}}$ (144)

and

$A_{\bf u}=\frac{(\kern-2.59995pt(0|(1-\Omega)\mathcal{M}_{\bf u}|Z^{\\{\rm
in\\}}\\}}{(\kern-2.59995pt(0|\mathcal{M}_{\bf u}|Z^{\\{\rm
in\\}}\\}}=1-\frac{(\kern-2.59995pt(0|\Omega\mathcal{M}_{\bf u}|Z^{\\{\rm
in\\}}\\}}{(\kern-2.59995pt(0|\mathcal{M}_{\bf u}|Z^{\\{\rm in\\}}\\}}$ (145)

From inspection of (133) we see that the absorption probabilities of (144) and
(145) vanish in the conservative scattering limit

$A_{\bf d}\to 0,\quad\hbox{and}\quad A_{\bf u}\to
0,\quad\hbox{as}\quad\tilde{\omega}\to 1.$ (146)

We define the transmissivity $T_{\bf u}$ and reflectivity $R_{\bf u}$ of the
upward input flux as

$T_{\bf u}=\frac{(\kern-2.59995pt(0|\mathcal{M}_{\bf u}\Omega\mathcal{M}_{\bf
u}|Z^{\\{\rm in\\}}\\}}{(\kern-2.59995pt(0|\mathcal{M}_{\bf u}|Z^{\\{\rm
in\\}}\\}},\quad\hbox{and}\quad R_{\bf
u}=\frac{(\kern-2.59995pt(0|\mathcal{M}_{\bf d}\Omega\mathcal{M}_{\bf
u}|Z^{\\{\rm in\\}}\\}}{(\kern-2.59995pt(0|\mathcal{M}_{\bf u}|Z^{\\{\rm
in\\}}\\}}.$ (147)

In like manner, the transmissivity $T_{\bf d}$ and reflectivity $R_{\bf d}$ of
the downward input flux are

$T_{\bf d}=\frac{(\kern-2.59995pt(0|\mathcal{M}_{\bf d}\Omega\mathcal{M}_{\bf
d}|Z^{\\{\rm in\\}}\\}}{(\kern-2.59995pt(0|\mathcal{M}_{\bf d}|Z^{\\{\rm
in\\}}\\}},\quad\hbox{and}\quad R_{\bf
d}=\frac{(\kern-2.59995pt(0|\mathcal{M}_{\bf u}\Omega\mathcal{M}_{\bf
d}|Z^{\\{\rm in\\}}\\}}{(\kern-2.59995pt(0|\mathcal{M}_{\bf d}|Z^{\\{\rm
in\\}}\\}}.$ (148)

Summing the transmission and reflection probabilities of (147) and using (114)
with (144) we find

$\displaystyle T_{\bf u}+R_{\bf u}$ $\displaystyle=$
$\displaystyle\frac{(\kern-2.59995pt(0|(\mathcal{M}_{\bf u}+\mathcal{M}_{\bf
d})\Omega\mathcal{M}_{\bf u}|Z^{\\{\rm
in\\}}\\}}{(\kern-2.59995pt(0|\mathcal{M}_{\bf u}|Z^{\\{\rm in\\}}\\}}$ (149)
$\displaystyle=$ $\displaystyle\frac{(\kern-2.59995pt(0|\Omega\mathcal{M}_{\bf
u}|Z^{\\{\rm in\\}}\\}}{(\kern-2.59995pt(0|\mathcal{M}_{\bf u}|Z^{\\{\rm
in\\}}\\}}$ $\displaystyle=$ $\displaystyle 1-A_{\bf u},$

or

$T_{\bf u}+R_{\bf u}+A_{\bf u}=1.$ (150)

The probabilities $T_{\bf u}$ for transmission, $R_{\bf u}$ for reflection,
and $A_{\bf u}$ for absorption of upward incoming flux sum to 1. In like
manner, one can show that for downward incoming flux

$T_{\bf d}+R_{\bf d}+A_{\bf d}=1.$ (151)

For conservative scattering, when absorption probabilities $A_{\bf d}$ and
$A_{\bf u}$ vanish, in accordance with (146), we can write (150) and (151) as

$T_{\bf d}+R_{\bf d}=1\quad\hbox{and}\quad T_{\bf u}+R_{\bf u}=1.$ (152)

###  5.1 Examples

The transmission and reflection of clouds depends on the cloud thickness
$\tau_{c}$, on the scattering phase function $p(\mu)$ of (35), and on the
directions of the incoming radiation. The directions are described by the
input intensity vector $|I^{\\{\rm in\\}}\\}$ of (111) or by its equivalent,
the input flux vector $|Z^{\\{\rm in\\}}\\}=4\pi(\hat{\mu}_{\bf
u}-\hat{\mu}_{\bf d})|I^{\\{\rm in\\}}\\}$ of (137).

For a model with $2n=10$ streams, Fig. 1 and Fig. 6 show details of how the
radiation penetrates a cloud of optical depth $\tau_{c}=20$ for the most
vertical possible input stream, with $|I^{\\{\rm in\\}}\\}=|\mu_{2n})$ and for
Rayleigh scattering.

In like manner, Fig. 2 and Fig. 7 show details of how the radiation penetrates
a cloud of optical depth $\tau_{c}=20$ for the most vertical possible input
stream, with $|I^{\\{\rm in\\}}\\}=|\mu_{2n})$, and for the maximum possible
forward scattering, with a phase function, $p(\mu)=\varpi^{\\{n\\}}(\mu)$ for
radiation described with $2n=10$ streams.

Fig. 8 and Fig. 9 show what happens if the input stream is the most horizontal
upward stream possible, with $|I^{\\{\rm in\\}}\\}=|\mu_{n+1})$. Input streams
that are nearly normal to the surfaces of the cloud, have significantly larger
diffuse transmission and smaller diffuse reflection than input streams that
are more nearly horizontal, and therefore must traverse longer slant distances
to get from one side of the cloud to the other.

Fig. 10 shows the moments of a half isotropic input stream
$|I^{\\{in\\}}\\}=\mathcal{M}_{\bf u}|0)$ incident on the bottom of a cloud
with the same optical depth, $\tau_{c}=20$ as those of Fig. 7 for a near
vertical stream and Fig. 8 for a near horizontal stream. The reflectivities
and transmissivities are intermediate between those of near vertical and near
horizontal streams. The intensity directions for Fig. 10 are shown as the red
rays of Fig. 13. The Eddington criterion (106) is much better satisfied near
the bottom of the cloud, because of the near isotropy of the intensity there.

Single-stream reflectivities calculated with (147) for upward incoming streams
$|I^{\\{\rm in\\}}\\}=|\mu_{i})$, with $i=n+1,n+2,\ldots,2n$, for various
cloud thicknesses are shown in Fig. 11 for Rayleigh scattering and in Fig. 12
for forward scattering.

##  6 Isotropic incident radiation

An interesting special case of conservative scattering is that of isotropic
input radiation, when upward radiation is incident with equal intensities from
all directions onto the bottom of the cloud, and downward radiation is
incident with equal intensities onto the top of the cloud. Then we can write
the incoming intensity vector (111) as

$|I^{\\{\rm in\\}}\\}=|0).$ (153)

According to (124) and (134), the outgoing intensity vector is

$\displaystyle|I^{\\{\rm out\\}}\\}$ $\displaystyle=$
$\displaystyle\mathcal{S}|I^{\\{\rm in\\}}\\}$ (154) $\displaystyle=$
$\displaystyle\mathcal{S}|0)$ $\displaystyle=$ $\displaystyle|0).$

For conservative scattering of an isotropic input intensity vector, the output
intensity vector is also isotropic. From (117) and (114) we see that the
intensity vector at the bottom of the cloud is

$\displaystyle|I(0)\\}$ $\displaystyle=$ $\displaystyle\mathcal{M}_{\bf
u}|I^{\\{\rm in\\}}\\}+\mathcal{M}_{\bf d}|I^{\\{\rm out\\}}\\}\noindent$
(155) $\displaystyle=$ $\displaystyle(\mathcal{M}_{\bf u}+\mathcal{M}_{\bf
d})|0)\noindent$ (156) $\displaystyle=$ $\displaystyle|0).$ (157)

In like manner, we find from (118) and (114) that the intensity vector at the
top of the cloud is

$|I(\tau_{c})\\}=|0).$ (158)

Inside the cloud, with $0\leq\tau\leq\tau_{c}$ we assume that the general
solution to the intensity vector, the expansion (68), includes only the two
quasi-isotropic bases $|v_{1})$ and $|v_{2n})$. We choose the amplitudes to be

$A_{i}=\frac{1}{3\eta_{1}\tau_{c}}(\delta_{i,1}+\delta_{i,2n}),$ (159)

so the intensity vector becomes

$\displaystyle|I(\tau)\\}$ $\displaystyle=$
$\displaystyle\sum_{i=1}^{2n}|v_{i}(\tau)\big{)}A_{i}$ (160) $\displaystyle=$
$\displaystyle\frac{|v_{1}(\tau)\big{)}+|v_{2n}(\tau)\big{)}}{3\eta_{1}\tau_{c}}$
$\displaystyle=$ $\displaystyle|0).$

One can derive the last line of (160) from the previous line by inspection of
(70) and (74). For $\tau=0$ and $\tau=\tau_{c}$, the intracloud intensity
vector (160) is consistent with (157) and (158). For isotropic incident
radiation, the intensity inside the cloud is also isotropic and has the same
magnitude as the incident radiation.

We can think of the intracloud intensity vector (160) as a part $|I^{\\{\bf
u\\}}(\tau)\\}$ generated by the upward part of the input intensity vector of
(153), $\mathcal{M}_{\bf u}|I^{\\{\rm in\\}}\\}=\mathcal{M}_{\bf u}|0)$,
incident on the bottom of the cloud and a part $|I^{\\{\bf d\\}}(\tau)\\}$
generated by the downward input intensity vector $\mathcal{M}_{\bf
d}|I^{\\{\rm in\\}}\\}=\mathcal{M}_{\bf d}|0)$. incident on the top of the
cloud. We use (126) to write the intensity vector from the upward incoming
intensity vector as

$\displaystyle|I^{\\{\bf u\\}}(\tau)\\}$ $\displaystyle=$
$\displaystyle\sum_{i^{\prime}i^{\prime\prime}}|v_{i^{\prime}}(\tau)\big{)}\mathcal{I}^{-1}_{i^{\prime}i^{\prime\prime}}(\kern-2.59995pt(\mu_{i^{\prime\prime}}|\mathcal{M}_{\rm
u}|0)$ $\displaystyle=$
$\displaystyle\sum_{i^{\prime}k}|v_{i^{\prime}}(\tau)\big{)}\mathcal{I}^{-1}_{i^{\prime}k}(\kern-2.59995pt(\mu_{k}|0)$
$\displaystyle=$
$\displaystyle\sum_{i^{\prime}k}|v_{i^{\prime}}(\tau)\big{)}\mathcal{I}^{-1}_{i^{\prime}k}w_{k}.$

Here the first summation index $i^{\prime}$ takes on the values
$i^{\prime}=1,2,3,\ldots,2n$. From the definition (113) of $\mathcal{M}_{\bf
u}$, we see that the second summation index takes on the values
$i^{\prime\prime}=k=n+1,n+2,n+3,\ldots,2n$. From (25) we see that
$(\kern-2.59995pt(\mu_{k}|0)=w_{k}$. The projections
$(\kern-2.59995pt(\mu_{i}|I^{\\{\bf u\\}}(\tau)\\}$ are shown as the red rays
in Figs. 13 and 14.

In like manner we use (126) to write the intensity vector from the downward
incoming intensity as

$\displaystyle|I^{\\{\bf d\\}}(\tau)\\}$ $\displaystyle=$
$\displaystyle\sum_{i^{\prime}i^{\prime\prime}}|v_{i^{\prime}}(\tau)\big{)}\mathcal{I}^{-1}_{i^{\prime}i^{\prime\prime}}(\kern-2.59995pt(\mu_{i^{\prime\prime}}|\mathcal{M}_{\rm
d}|0)$ $\displaystyle=$
$\displaystyle\sum_{i^{\prime}j}|v_{i^{\prime}}(\tau)\big{)}\mathcal{I}^{-1}_{i^{\prime}j}w_{j}.$

Here the first summation index $i^{\prime}$ takes on the values
$i^{\prime}=1,2,3,\ldots,2n$. From the definition (113) of $\mathcal{M}_{\bf
d}$, we see that the second summation index takes on the values
$i^{\prime\prime}=j=1,2,3,\ldots,n$. The projections
$(\kern-2.59995pt(\mu_{i}|I^{\\{\bf d\\}}(\tau)\\}$ are shown as the blue rays
in Figs. 13 and 14.

In summary, for the isotropic incoming intensity vector (153), the intensities
inside a cloud remain isotropic and equal to the input intensity at any
optical depth $\tau$ above the bottom, for any scattering phase $p(\mu)$, and
for any cloud thickness $\tau_{c}$. For real absorbing clouds of sufficiently
large optical thickness $\tau_{c}\gg 1$, the intensity attenuates to
negligibly small values near the middle.

##  7 Thin Clouds

For infinitesimally thin clouds of optical depth $\tau_{c}=d\tau$ the
scattering matrix simplifies to

$\mathcal{S}=\hat{1}-d\tau(\hat{\varsigma}_{\bf u}-\hat{\varsigma}_{\bf
d})\hat{\eta}.$ (163)

in accordance with (WH-184). The efficiency matrix $\hat{\eta}$ was given by
(30). From (128) we find that the cloud albedo matrix becomes

$\displaystyle\Omega$ $\displaystyle=$ $\displaystyle(\hat{\mu}_{\bf
u}-\hat{\mu}_{\bf d})\mathcal{S}(\hat{\varsigma}_{\bf u}-\hat{\varsigma}_{\bf
d})$ (164) $\displaystyle=$
$\displaystyle\hat{1}-d\tau\hat{\eta}(\hat{\varsigma}_{\bf
u}-\hat{\varsigma}_{\bf d}).$

From (145) we can write the absorption probability for the upward part of the
incoming radiation as

$\displaystyle A_{\bf u}$ $\displaystyle=$ $\displaystyle
1-\frac{(\kern-2.59995pt(0|\\{\hat{1}-d\tau\hat{\eta}(\hat{\varsigma}_{\bf
u}-\hat{\varsigma}_{\bf d})\\}\mathcal{M}_{\bf u}|Z^{\\{\rm
in\\}}\\}}{(\kern-2.59995pt(0|\mathcal{M}_{\bf u}|Z^{\\{\rm in\\}}\\}}$ (165)
$\displaystyle=$ $\displaystyle
d\tau\frac{(\kern-2.59995pt(0|\hat{\eta}\hat{\varsigma}_{\bf u}|Z^{\\{\rm
in\\}}\\}}{(\kern-2.59995pt(0|\mathcal{M}_{\bf u}|Z^{\\{\rm in\\}}\\}}$
$\displaystyle=$ $\displaystyle
d\tau(1-\tilde{\omega})\frac{(\kern-2.59995pt(0|\hat{\varsigma}_{\bf
u}|Z^{\\{\rm in\\}}\\}}{(\kern-2.59995pt(0|\mathcal{M}_{\bf u}|Z^{\\{\rm
in\\}}\\}}$

Here we noted from (42) that
$(\kern-2.59995pt(0|\hat{\eta}=(1-\tilde{\omega})(\kern-2.59995pt(0|$. From
(165) we see that the absorption $A_{\bf u}$ of a thin cloud depends on the
direction of the input flux, as described by $|Z^{\\{\rm in\\}}\\}$, but is
independent of the scattering phase function.

As a simple example, suppose that the incoming flux vector consists of a
single upward stream of index $k>n$,

$|Z^{\\{\rm in\\}}\\}=4\pi\mu_{k}|\mu_{k})$ (166)

Substituting (166) into (165) we find that the absorption is

$\displaystyle A_{\bf u}$ $\displaystyle=$ $\displaystyle
d\tau(1-\tilde{\omega})\frac{(\kern-2.59995pt(0|\hat{\varsigma}_{\bf
u}|\mu_{k})}{(\kern-2.59995pt(0|\mathcal{M}_{\bf u}|\mu_{k})}$ (167)
$\displaystyle=$ $\displaystyle
d\tau\varsigma_{k}(1-\tilde{\omega})\frac{(\kern-2.59995pt(0|\mu_{k})}{(\kern-2.59995pt(0|\mu_{k})}$
$\displaystyle=$ $\displaystyle d\tau\varsigma_{k}(1-\tilde{\omega}).$

According to (167), for an input intensity vector $|I^{\\{\rm
in\\}}\\}=|\mu_{k})$, corresponding to the input flux vector (166), the
absorption $A_{\bf u}$ is proportional to $d\tau\varsigma_{k}$, the slant
distance of the input stream through the cloud, and to the probability
$(1-\tilde{\omega})$ that a photon collision with a cloud particulate leads to
absorption rather than scattering.

As another example, suppose that the incoming intensity vector at the bottom
of the cloud is half isotropic, $|I^{\\{\rm in\\}}\\}=\mathcal{M}_{\bf u}|0)$,
so that the incoming flux vector (137) becomes

$|Z^{\\{\rm in\\}}\\}=4\pi\hat{\mu}_{\bf u}|0).$ (168)

Substituting (168) into (165) and noting that
$(\kern-2.59995pt(0|\mathcal{M}_{\bf u}|0)=1/2$, as follows from (179) below,
we find

$\displaystyle A_{\bf u}$ $\displaystyle=$ $\displaystyle
d\tau(1-\tilde{\omega})\frac{(\kern-2.59995pt(0|\mathcal{M}_{\bf
u}|0)}{(\kern-2.59995pt(0|\hat{\mu}_{\bf u}|\mu_{k})}$ (169) $\displaystyle=$
$\displaystyle\frac{d\tau(1-\tilde{\omega})}{2(\kern-2.59995pt(0|\hat{\mu}_{\bf
u}|0)}$

The matrix element in the denominator of (169) can be written as

$\displaystyle(\kern-2.59995pt(0|\hat{\mu}_{\bf u}|0)$ $\displaystyle=$
$\displaystyle\frac{1}{2}\sum_{k=n+1}^{2n}w_{k}\mu_{k}$ (170)
$\displaystyle\to$ $\displaystyle\frac{1}{2}\int_{0}^{1}d\mu\,\mu$
$\displaystyle=$ $\displaystyle\frac{1}{4},\quad\hbox{as}\quad n\to\infty.$

Here we made use of the limit of the Gauss-Legendre quadrature,

$\lim_{n\to\infty}\sum_{i=1}^{2n}w_{i}f(\mu_{i})=\int_{-1}^{1}d\mu f(\mu).$
(171)

Substituting (170) into (169) we find

$\lim_{n\to\infty}A_{\bf u}=2\,d\tau(1-\tilde{\omega}).$ (172)

For quasi-isotropic illumination from one side, the effective thickness of a
thin cloud is $2d\tau$.

Figure 15: Reflection efficiencies for $2n=32$ streams of half-isotropic
radiation (168), incident onto the bottom of a thin cloud. The small blue
circles, $a^{\\{p\\}}_{f}$, given by (183) are for the maximum forward
scattering phase function $p(\mu)=\varpi^{\\{p\\}}(\mu)$ of (38). The small
red circles, given by (184), are for the the maximum backward scattering phase
function $p(\mu)=\varpi^{\\{p\\}}(-\mu)$. The values of
$a^{\\{p\\}}_{f}=a^{\\{p\\}}_{b}=1$ for $p=0$ are for isotropic scattering. As
one would expect, strongly peaked forward-scattering phase functions,
constructed from many pairs $p$ of Legendre polynomials, suppress the thin
cloud reflection efficiency $a^{\\{p\\}}_{f}$, compared to isotropic
scattering. Strongly backward-peaked phase functions can nearly double the
reflection efficiency $a^{\\{p\\}}_{b}$. For $2n=10$ streams the coefficients
$a^{\\{p\\}}_{f}$ and $a^{\\{p\\}}_{b}$, shown as large circles, are almost
the same as for $2n=32$ streams. The largest difference is for $p=5$, where
for $2n=10$ streams we find $a_{f}^{\\{5\\}}=0.2759$ and for $2n=32$ streams
we find $a_{f}^{\\{5\\}}=0.2968$.

From (147) we see that the reflected fraction of the half-isotropic input flux
vector (168) is

$\displaystyle R_{\bf u}$ $\displaystyle=$
$\displaystyle\frac{(\kern-2.59995pt(0|\mathcal{M}_{\bf d}\Omega\hat{\mu}_{\bf
u}|0)}{(\kern-2.59995pt(0|\hat{\mu}_{\bf u}|0)}$ (173) $\displaystyle=$
$\displaystyle\frac{(\kern-2.59995pt(0|\mathcal{M}_{\bf d}(\hat{\mu}_{\bf
u}-d\tau\hat{\eta}\mathcal{M}_{\bf u})|0)}{(\kern-2.59995pt(0|\hat{\mu}_{\bf
u}|0)}$ $\displaystyle=$
$\displaystyle-d\tau\frac{(\kern-2.59995pt(0|\mathcal{M}_{\bf
d}\hat{\eta}\mathcal{M}_{\bf u}|0)}{(\kern-2.59995pt(0|\hat{\mu}_{\bf u}|0)}$
$\displaystyle=$ $\displaystyle-d\tau\frac{(\kern-2.59995pt(0|\mathcal{M}_{\bf
d}(\hat{1}-\tilde{\omega}\hat{p}/2)\mathcal{M}_{\bf
u}|0)}{(\kern-2.59995pt(0|\hat{\mu}_{\bf u}|0)}$ $\displaystyle=$
$\displaystyle d\tau\frac{\tilde{\omega}(\kern-2.59995pt(0|\mathcal{M}_{\bf
d}\hat{p}\mathcal{M}_{\bf u}|0)}{2(\kern-2.59995pt(0|\hat{\mu}_{\bf u}|0)}$
$\displaystyle=$ $\displaystyle
d\tau\tilde{\omega}\sum_{l=0}^{2n-1}p_{l}\frac{(\kern-2.59995pt(0|\mathcal{M}_{\bf
d}|l)(\kern-2.59995pt(l|\mathcal{M}_{\bf
u}|0)}{(\kern-2.59995pt(0|\hat{\mu}_{\bf u}|0)}.$

To get the second line of (173) from the first we used (164) and the fact that
$(\hat{\varsigma}_{\bf u}-\hat{\varsigma}_{\bf d})\hat{\mu}_{\bf
u}=\mathcal{M}_{\bf u}$. The third line follows from the second since
$\mathcal{M}_{\bf d}\hat{\mu}_{\bf u}=\hat{0}$. To get the fourth line from
the third we used the expression (30) for the efficiency matrix $\hat{\eta}$.
To get the fifth line from the fourth, we noted from (116) that
$\mathcal{M}_{\bf d}\mathcal{M}_{\bf u}=\hat{0}$. To get the last line from
the fifth, we used the multipole expansion (32) of the phase matrix $\hat{p}$.

Using (24) and (25) with (113) we write

$\displaystyle(\kern-2.59995pt(l|\mathcal{M}_{\bf u}|0)$ $\displaystyle=$
$\displaystyle\sum_{k=n+1}^{2n}(\kern-2.59995pt(l|\mu_{k})(\kern-2.59995pt(\mu_{k}|0)$
(174) $\displaystyle=$
$\displaystyle\frac{1}{2}\sum_{k=n+1}^{2n}w_{k}P_{l}(\mu_{k}).$

In like manner we write

$\displaystyle(\kern-2.59995pt(0|\mathcal{M}_{\bf d}|l)$ $\displaystyle=$
$\displaystyle\sum_{j=1}^{n}(\kern-2.59995pt(0|\mu_{j})(\kern-2.59995pt(\mu_{j}|l)$
(175) $\displaystyle=$
$\displaystyle\frac{2l+1}{2}\sum_{j=1}^{n}w_{j}P_{l}(\mu_{j})$
$\displaystyle=$
$\displaystyle(-1)^{l}\frac{2l+1}{2}\sum_{k=n+1}^{2n}w_{k}P_{l}(\mu_{k})$
$\displaystyle=$
$\displaystyle(-1)^{l}(2l+1)(\kern-2.59995pt(l|\mathcal{M}_{\bf u}|0).$

For reflection-conjugate pairs of indices $j$ and $k$, with $k=2n+1-j$, we
noted the symmetries $\mu_{k}=-\mu_{j}$, $w_{j}=w_{k}$ and
$P_{l}(\mu_{j})=(-1)^{l}P_{l}(\mu_{k})$, to get the third line from the
second. We used (174) to write the last line of (175).

Using (175), we write (173) as

$\displaystyle R_{\bf u}$ $\displaystyle=$
$\displaystyle\frac{d\tau\tilde{\omega}}{(\kern-2.59995pt(0|\hat{\mu}_{\bf
u}|0)}\sum_{l=0}^{2n-1}(-1)^{l}(2l+1)p_{l}(\kern-2.59995pt(l|\mathcal{M}_{\bf
u}|0)^{2}.$ (176)

Using the symmetries mentioned above, we see from (174) that

$(\kern-2.59995pt(l|\mathcal{M}_{\bf
u}|0)=(-1)^{l}(\kern-2.59995pt(l|\mathcal{M}_{\bf d}|0).$ (177)

Using (2), (114) and (177) we find

$\displaystyle\delta_{l0}$ $\displaystyle=$
$\displaystyle(\kern-2.59995pt(l|0)=(\kern-2.59995pt(l|\hat{1}|0)$ (178)
$\displaystyle=$ $\displaystyle(\kern-2.59995pt(l|\mathcal{M}_{\bf
u}+\mathcal{M}_{\bf d}|0)$ $\displaystyle=$
$\displaystyle(\kern-2.59995pt(l|\mathcal{M}_{\bf u}|0)[1+(-1)^{l}].$

For even values of the multipole index $l$ we can write (178)

$(\kern-2.59995pt(l|\mathcal{M}_{\bf
u}|0)=\frac{\delta_{l0}}{2},\quad\hbox{for}\quad l=0,2,4,\ldots.$ (179)

We can therefore write the reflectivity (176) as the product of a reflectivity
$R_{\bf u}^{\\{\rm iso\\}}$ of a thin cloud with isotropic scattering and an
anisotropic scattering efficiency $a$,

$R_{\bf u}=R_{\bf u}^{\\{\rm iso\\}}a.$ (180)

Using (169), we find that the reflectivity for isotropic scattering, with
$p_{l}=\delta_{l0}$, is

$R_{\bf u}^{\\{\rm
iso\\}}=\frac{d\tau\tilde{\omega}}{4(\kern-2.59995pt(0|\hat{\mu}_{\bf
u}|0)}=\frac{\tilde{\omega}}{2(1-\tilde{\omega})}A_{\bf
u}\approx\tilde{\omega}d\tau$ (181)

The anisotropic scattering efficiency is

$\displaystyle
a=1-4\sum_{l^{\prime}}(2l^{\prime}+1)p_{l^{\prime}}(\kern-2.59995pt(l^{\prime}|\mathcal{M}_{\bf
u}|0)^{2},\quad\hbox{where}\quad l^{\prime}=1,3,5,\ldots,2n-1.$ (182)

Formulas for downward radiation can be obtained by replacing the index $\\{\bf
u\\}$ by $\\{\bf d\\}$ in (180)–(182).

Examples of (180) are shown in Fig. 15 for the maximum forward-scattering
phase functions, $p(\mu)=\varpi^{\\{p\\}}(\mu)$ of (38). We denote the
anisotropic efficiencies for maximum forward scattering by

$\displaystyle
a^{\\{p\\}}_{f}=1-4\sum_{l^{\prime}}(2l^{\prime}+1)\varpi^{\\{p\\}}_{l^{\prime}}(\kern-2.59995pt(l^{\prime}|\mathcal{M}_{\bf
u}|0)^{2},\quad\hbox{where}\quad l^{\prime}=1,3,5,\ldots,2n-1.$ (183)

The anisotropic efficiencies for maximum backward scattering, where the
scattering phase is $p(\mu)=\varpi^{\\{p\\}}(-\mu)$, are

$\displaystyle
a^{\\{p\\}}_{b}=1+4\sum_{l^{\prime}}(2l^{\prime}+1)\varpi^{\\{p\\}}_{l^{\prime}}(\kern-2.59995pt(l^{\prime}|\mathcal{M}_{\bf
u}|0)^{2},\quad\hbox{where}\quad l^{\prime}=1,3,5,\ldots,2n-1.$ (184)

We see that

$\displaystyle a^{\\{p\\}}_{f}+a^{\\{p\\}}_{b}=2.$ (185)

As one can see from Fig. 15, the numerical coefficients $a^{\\{p\\}}_{f}$ and
$a^{\\{p\\}}_{b}$ are nearly the same for radiation transfer represented by
$n=5$ stream pairs, shown as large circles, and $n=16$ stream pairs, shown as
small circles. The number of Legendre-polynomial pairs $p$ in the phase
functions $\varpi^{\\{p\\}}(\mu)$ cannot exceed $n$, that is $p\leq n$.

##  8 Summary

We have shown how to use matrix methods of quantum mechanics to efficiently
and accurately calculate axially symmetric radiation transfer in clouds, with
conservative multiple scattering, and with arbitrary single-scattering
anisotropies or phase functions.

At a given optical thickness, $\tau$, we model the radiation with the $2n$
unweighted intensities $I(\mu_{i},\tau)$ of (10), which are values of the
intensity $I(\mu,\tau)$ at the Gauss-Legendre direction cosines $\mu=\mu_{i}$.
The $\mu_{i}$ are the roots of the Legendre polynomial $P_{2n}(\mu)$.
According to (7), the $\mu_{i}$ are also the eigenvalues of the $2n\times 2n$
direction-cosine matrix $\hat{\mu}$ of (18). As shown in (10), it is
convenient to write the unweighted intensities as
$I(\mu_{i},\tau)=w_{i}^{-1}(\kern-2.59995pt(\mu_{i}|I(\tau)\\}$ where the
Gauss-Legendre weights can be defined by (29). Here the weighted intensity
$(\kern-2.59995pt(\mu_{i}|I(\tau)\\}$ denotes the projection of an abstract
intensity vector $|I(\tau)\\}$ onto the left stream vector
$(\kern-2.59995pt(\mu_{i}|$ of (7).

For clouds with non-zero absorption and $\tilde{\omega}<1$, WH showed how to
use matrix methods to model $2n$-stream radiative transfer with arbitrary
scattering anisotropy. In WH the intensity vector $|I(\tau)\\}$ of (1) was
expanded onto $2n$ penetration-modes $|\lambda_{i})$ which depend on optical
depth $\tau$ as $e^{-\tau/\lambda_{i}}$. Here $\lambda_{i}$ is the penetration
length and $|\lambda_{i})$ is the right eigenvector of the penetration-length
matrix $\hat{\lambda}=\hat{\kappa}^{-1}$ the inverse of the exponentiation-
rate matrix $\hat{\kappa}$ of (56). When $\tilde{\omega}\to 1$, and absorption
vanishes entirely, two of the penetration lengths diverge to infinity,
$\lambda_{1}\to-\infty$ and $\lambda_{2n}\to\infty$. Then it is necessary to
replace the two infinitely-penetrating modes $|\lambda_{1})$ and
$|\lambda_{2n})$ with the quasi-isotropic bases of $|v_{1}\big{)}$ of (70) and
$|v_{2n})$ of (74) that depend linearly, not exponentially on $\tau$. For
$i=2,3,4,\ldots 2n-1$ the directional bases $|v_{i})$ are the limits of the
penetration modes of WH for $\tilde{\omega}\to 1$, that is
$|\lambda_{i})\to|v_{i}(\tau_{i})\big{)}$, where the reference optical depths
$\tau_{i}$ were given by (95). Quasi-isotropic and directional bases are
illustrated in Fig. 3 and 5. The dependence of the penetration lengths
$\lambda_{i}$ on the sharpness of the forward scattering phase function
$\varpi^{\\{p\\}}(\mu)$ of (38) is shown in Fig. 4.

For conservative scattering, the flux $Z$ and the proportional dipole moment
$I_{1}$ of the intensity are independent of optical depth $\tau$ in the cloud,
as discussed in connection with (52). Also, the $K$-integral, or the second
moment $K$ of the intensity, given by (47), varies linearly with $\tau$, as
shown by (54). The invariance of $Z$ and the linear variation of $K$ with
$\tau$ are illustrated with the moment diagrams of Figs. 1, 2, 9 and 10.

As directional radiation penetrates a cloud, multiple conservative scattering
causes the radiation to isotropize within a few optical depths of the input
surface. The isotropization is illustrated with the ray diagrams of Figs. 6,
7, and 8. The Eddington criterion (106), that the higher multipole moments,
$I_{2}(\tau),I_{3}(\tau),I_{4}(\tau),\ldots$, are negligible compared to the
first two moments, $I_{0}(\tau)$ and $I_{1}(\tau)$, is very well satisfied
near the centers of optically thick clouds, as one can see from the dashed red
curves of Figs. 1, 2, 9 and 10.

The transmission and reflection of clouds depends strongly on the angle of
incidence of incoming radiation. For the same thickness, $\tau_{c}=20$, and
phase function, $p(\mu)=\varpi^{\\{5\\}}(\mu)$, a cloud transmits about three
times more near-normal-incidence radiation, shown in Fig. 2 or Fig. 7, than
near-horizontal-incidence radiation of Fig. 8 or Fig. 9. Fig. 11 shows how the
reflectance of a cloud with the Rayleigh-scattering phase function of (36)
depends on the direction cosine $\mu_{i}$ of the input radiation and on the
cloud thickness $\tau_{c}$. Fig. 12 shows that a cloud with the more strongly
forward scattering phase function $p(\mu)=\varpi^{\\{5\\}}(\mu)$ of (38)
reflects less (and transmits more).

For a cloud illuminated from the bottom with equal intensity from all upward
directions and illuminated from the top with the same equal intensity from all
downward directions, the intensity that results from multiple scattering is
isotropic and of equal intensity inside and outside the cloud. This
“conservation of isotropy” is independent of the scattering phase function as
is illustrated in Fig. 13 for forward scattering in a cloud with a relatively
large optical depth, $\tau_{c}=20$ and in Fig. 14 for a thinner cloud with
$\tau_{c}=1$. Even for $\tau_{c}\to\infty$, the cloud will remain filled with
isotropic radiation with equal intensities at all depths. Real clouds of
sufficiently large optical depth $\tau_{c}$, and with finite absorption,
$1-\tilde{\omega}>0$ have negligible intensity at their centers,
$\tau=\tau_{c}/2$, under the same conditions of half-isotropic illumination
from the top and bottom. The input radiation from both the top and bottom is
absorbed before reaching the center.

Fig. 15 shows how the reflectance of a thin cloud with optical thickess
$d\tau\ll 1$ decreases with increasing numbers $p$ of Legendre-polynomial
pairs in the maximum forward-scattering phase functions
$p(\mu)=\varpi^{\\{p\\}}(\mu)$ of (38). The forward-scattering phases are
$\varpi^{\\{p\\}}(1)=p(p+1)$, so the integer $p=1,2,3,\ldots$ measures the
forward peaking of the phase functions. Fig. 15 also shows that increasingly
sharply peaked backscattering phase functions $p(\mu)=\varpi^{\\{p\\}}(-\mu)$
increase the reflectance toward a limiting value, $\mathcal{R}_{\bf u}\to
2\mathcal{R}^{\\{\rm iso\\}}_{\bf u}$ that is twice that of an isotropically
scattering phase function, $\mathcal{R}^{\\{\rm iso\\}}_{\bf u}$, given by
(181).

The number $2n$ of streams needed to make accurate calculations is often
surprsisingly small. An example is shown in Figs. 11 and 12, which show the
reflectance of clouds versus the angle of incidence of input radiation.
Calculations made with $2n=32$ streams are indicated with small circles, and
those made with $2n=10$ streams are shown as large circles. The results are
virtually the same. Another example is provided by Fig. 15. For $p=5$, the
maximum value possible for $n=5$, the results for $n=5$ and $n=16$ stream
pairs differ by about 7%.

Using the efficient matrix methods outlined above with the capabilities of
modern mathematical software makes calculations of radiation transfer fast and
simple. The figures of this paper were generated with a few dozen lines of
Matlab code.

## Acknowledgements

The Canadian Natural Science and Engineering Research Council provided
financial support of one of us. We are grateful to Professor C. A. de Lange
for valuable discussions of this work.

## References

  * [1] H. Frisch, Radiative Transfer, Springer Nature Switzerland AG (2022).
  * [2] G. C. Wick, Über ebene Diffusionsprobleme, Z. Physik, 121, 702 (1943).
  * [3] W. A. van Wijngaarden and W. Happer, 2n-Stream Radiative Transfer , https://arxiv.org/pdf/2205.09713.pdf
  * [4] S. Chandrasekhar, Radiative Transfer, Dover, New York (1960).
  * [5] M. Abramowitz and I. Stegun, Handbook of Mathematical Functions, Dover Publications, New York (1965).
  * [6] Neutron reflectors, https://ansn.iaea.org/Common/documents/Training/TRIGA%20Reactors%20%28Safety%20and%20Technology%29/chapter2/physics121.htm