# Radiation Transfer in Cloud Layers

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

###### Abstract

For $2n$-stream radiation transfer theory, a stack of $m=2,3,4,\ldots$ clouds
can be represented as an equivalent cloud. Individual clouds, indexed by
$c=1,2,3,\ldots,m$, are characterized by $2n\times 2n$ scattering matrices
$\mathcal{S}^{\\{c\\}}$, that describe how the cloud interacts with $2n$
streams of axially symmetric incoming radiation, propagating in upward and
downward Gauss-Legendre sample directions. Some of the radiation is
transmitted, some is absorbed and converted to heat, and some is scattered
into $2n$ outgoing streams along the directions of the incoming streams. The
clouds are also characterized by $2n\times 1$ thermal source vectors
$|\dot{J}^{\\{c\\}}\\}$ that describe the thermal emission of radiation along
the stream directions by cloud particulates and gas molecules. The $2n\times
2n$ scattering matrix for the equivalent cloud, $\mathcal{S}^{\\{\rm
ev\\}}=S^{\\{m\\}}\raisebox{-1.0pt}{*}S^{\\{m-1\\}}\raisebox{-1.0pt}{*}\cdots
S^{\\{2\\}}\raisebox{-1.0pt}{*}S^{\\{1\\}}$, is the Redheffer star product of
the scattering matrices $\mathcal{S}^{\\{c\\}}$ of the individual clouds. The
$2n\times 1$ thermal source vector for the equivalent cloud, $|\dot{J}^{\\{\rm
ev\\}}\\}=G^{[m,m\\}}|\dot{J}^{\\{m\\}}\\}+G^{[m,m-1\\}}|\dot{J}^{\\{m-1\\}}\\}+\cdots+G^{[m,2\\}}|\dot{J}^{\\{2\\}}\\}+G^{[m,1\\}}|\dot{J}^{\\{1\\}}\\}$,
is a linear combination of the thermal source vectors $|\dot{J}^{\\{c\\}}\\}$
of the individual clouds. The $2n\times 2n$ discrete Green’s matrices
$G^{[m,c\\}}$ can be constructed from the scattering matrices
$\mathcal{S}^{\\{c\\}}$ of individual clouds. The equivalent scattering matrix
$\mathcal{S}^{\\{\rm ev\\}}$ and the equivalent thermal source vector
$|\dot{J}^{\\{\rm ev\\}}\\}$ are analogous to the equivalent resistance and
the equivalent electromotive force (emf) of Thévenin’s theorem for a network
of electrical circuits. Illustrative numerical examples are given for single
clouds, 3-cloud stacks and 10-cloud stacks. These methods are useful for
modeling radiation transfer in Earth’s atmosphere, which can be represented by
layers of invisible clouds, consisting of clear air and greenhouse gases, or
visible clouds which also include condensed water, smog, etc.

Keywords:  radiative transfer, multiple scattering, reflection, absorption,
emission, scatttering phase functions, equation of transfer, Gauss-Legendre
quadrature, two-port networks

###### Contents

  1. 1 Introduction
  2. 2 Radiation Intensity and Flux
    1. 2.1 The $2n$-stream basis
    2. 2.2 The multipole basis
    3. 2.3 Equation of radiative transfer
    4. 2.4 Vertical flux
    5. 2.5 Outgoing and incoming radiation
    6. 2.6 Incident and thermal radiation
    7. 2.7 Identities for scattering matrices $\mathcal{S}$
    8. 2.8 Scattering matrix $\mathcal{S}$ for a homogeneous cloud
    9. 2.9 Continuous Green’s matrix $G(\tau^{\prime})$
    10. 2.10 Black clouds
    11. 2.11 Planck emissivities
    12. 2.12 Planck albedos
    13. 2.13 Half isotropic incoming intensity
    14. 2.14 Radiative heating and cooling of an isolated cloud
    15. 2.15 Radiative equilibrium
    16. 2.16 1-cloud examples
  3. 3 Cloud Stacks
    1. 3.1 Stack analysis
    2. 3.2 Thermal equilibrium
    3. 3.3 Formal solutions
    4. 3.4 Gap intensities and fluxes
    5. 3.5 Model cloud stacks
    6. 3.6 Radiative heating and cooling of cloud stacks
    7. 3.7 3-cloud examples
    8. 3.8 10-cloud examples
  4. 4 A Double Cloud
    1. 4.1 The matrix $(\hat{1}-X)^{-1}$ for $m=2$
    2. 4.2 Discrete Green’s matrices
    3. 4.3 Scattering coefficients
    4. 4.4 Star products
  5. 5 Summary

##  1 Introduction

In a previous paper, 2n-Stream Radiative Transfer [1], we outlined how to use
operator and matrix methods to accurately and efficiently analyze radiative
transfer with angular distributions of scattering, including highly forward
scattering of sunlight by cloud particulates and more nearly isotropic
Rayleigh scattering by gas molecules in Earth’s atmosphere. In a subsequent
paper, $2n$-Stream Conservative Scattering [2], we showed how to extend these
methods to the limit of conservative scattering, where photons can only be
transmitted or scattered, but not absorbed or emitted. A third paper of this
series, 2n-Stream Thermal Emission from Clouds  [3], was devoted to
Kirchhoff’s laws of transmission, absorption, scattering and emission of
thermal radiation by clouds. The $2n$-stream method used in these papers is a
generalization of the 2-stream method of Schuster [4] in his classic paper of
1905, Radiation Through a Foggy Atmosphere.

Our three previous papers on radiation transfer dealt mostly with homogeneous
clouds, where the single-scattering albedo $\tilde{\omega}$ and the scattering
phase function $p(\mu,\mu^{\prime})$ were spatially uniform. In the present
paper we will extend this discussion to stacks of individual clouds, $c$, each
with its own scattering matrix $\mathcal{S}^{\\{c\\}}$, and thermal source
vector $|\dot{J}^{\\{c\\}}\\}$. When scattering incoming radiation, a stack
acts like a single cloud with an equivalent scattering matrix
$\mathcal{S}^{\\{\rm
ev\\}}=S^{\\{m\\}}\raisebox{-1.0pt}{*}S^{\\{m-1\\}}\raisebox{-1.0pt}{*}\cdots
S^{\\{2\\}}\raisebox{-1.0pt}{*}S^{\\{1\\}}$, the Redheffer star product of the
scattering matrices $\mathcal{S}^{\\{c\\}}$ of the individual clouds. When
emitting thermal radiation, a stack acts like a single cloud with an
equivalent thermal source vector $|\dot{J}^{\\{\rm
ev\\}}\\}=G^{[m,m\\}}|\dot{J}^{\\{m\\}}\\}+G^{[m,m-1\\}}|\dot{J}^{\\{m-1\\}}\\}+\cdots+G^{[m,2\\}}|\dot{J}^{\\{2\\}}\\}+G^{[m,1\\}}|\dot{J}^{\\{1\\}}\\}$,
a linear combination of the thermal source vectors $|\dot{J}^{\\{c\\}}\\}$ of
the individual clouds. The $2n\times 2n$ discrete Green’s matrices
$G^{[m,c\\}}$ can be constructed from the individual scattering matrices
$\mathcal{S}^{\\{c\\}}$. Stacks of clouds have close analogies to networks of
two-port electrical circuits connected in series [5]. The equivalent
scattering matrix $\mathcal{S}^{\\{\rm ev\\}}$ and the equivalent thermal
source vector $|\dot{J}^{\\{\rm ev\\}}\\}$ are analogous to the equivalent
resistance and the equivalent electromotive force (emf) of a network of
electrical circuits that follow from Thévenin’s theorem [6].

An instructive history of radiative transfer theory has been given by
Mobley[7]. To our knowledge, the earliest paper where the $2n$-stream method
was used to analyze multiple scattering was published in 1943 by G. C. Wick in
connection with his analysis of neutron diffusion, Über ebene
Diffusionsprobleme [8]. Wick’s work has led to an extensive literature on
variants of the $2n$-stream computational methods, often described as the
discrete ordinate method (DOM). A useful review of DOM work has been given by
Ganapol[11] in 2015 as The Response Matrix Discrete Ordinates Solution to the
1D Radiative Transfer Equation. An early application of the $2n$ stream method
to describe radiation transport in clouds was published by Flannery, Roberge
and Rybicki [10] in 1980 as The Penetration of Diffuse Ultraviolet Radiation
into Stellar Clouds. We will frequently refer to the authoritative review of
radiative transfer in Chandrasekhar’s classic book, Radiative Transfer[12].
The book by Thomas and Stamnes[13], Radiation Transfer in the Atmosphere and
Ocean has extensive discussions of discrete ordinate methods.

##  2 Radiation Intensity and Flux

For axial symmetry about the zenith direction, we will characterize time-
independent, steady state radiation of spatial frequency $\nu$ at an altitude
$z$ above Earth’s surface with the monochromatic intensity $I(\nu,\mu,z)$,
also called the radiance. One can think of the intensity as streams of photons
making various angles, $\theta=\cos^{-1}\mu$, with the vertical.
$I(\nu,\mu,z)\,d\mu\,d\nu$ is the radiative flux carried by photons with
direction cosines between $\mu$ and $\mu+d\mu$ and with spatial frequencies
between $\nu$ and $\nu+d\nu$. A representative unit of $I(\nu,\mu,z)$ is W m-2
cm sr-1, where W= watts, is the unit of radiation power, m2 = square meters,
is the unit of irradiated area, cm-1 = waves per cm, is the unit of spatial
frequency of the radiation, and sr = steradian, is the unit of solid angle. In
2.1(3), Chandrasekhar[12] uses the symbol $I_{\nu}$ to denote our intensity
$I=I(\nu,\mu,z)$. For most of the remainder of this paper we will discuss only
monochromatic radiation and we will usually omit the frequency variable $\nu$
and write $I(\nu,\mu,z)=I(\mu,z)$.

The intensity $I(\mu,\tau)$ obeys the steady-state equation of transfer,

$\left(\mu\frac{\partial}{\partial\tau}+1\right)I(\mu,\tau)=(1-\tilde{\omega})B(\tau)+\frac{\tilde{\omega}}{2}\int_{-1}^{1}d\mu^{\prime}\,p(\mu,\mu^{\prime})I(\mu^{\prime},\tau).$
(1)

We neglect any variation of the intensity in horizontal spatial directions.
Chandrasekhar[12] writes (1) as §6(47). He uses the symbol ${\Im}_{\nu}$ to
denote both source terms on the right of (1). He writes the emissive part of
the source, which is proportional to the Planck intensity $B(\tau)$, as
§5(42). He writes the scattering part of the source, that is proportional to
the scattering phase $p(\mu,\mu^{\prime})$, as §5(41). The optical depth
$\tau=\tau(z)$ at an altitude $z$ above the bottom of the cloud is

$\tau=\int_{0}^{z}dz^{\prime}\,\alpha(z^{\prime}),$ (2)

The net attenuation rate, due to absorption and scattering, at the altitude
$z$ is $\alpha=\alpha(z)$. We see that $d\tau=\alpha\,dz$.

The source terms on the right of (1) are characterized by the single-
scattering albedo $\tilde{\omega}$, by the Planck intensity $B(\tau)$, and by
the scattering phase function $p(\mu,\mu^{\prime})$. Both $\tilde{\omega}$ and
$p(\mu,\mu^{\prime})$ may depend on altitude $z$, or equivalently, on the
vertical optical depth $\tau$. The single-scattering albedo $\tilde{\omega}$
is the probability that a photon, after a collision with a molecule or cloud
particulate, is elastically scattered into other directions. A fraction
$1-\tilde{\omega}$ of the photons is absorbed and converted to atmospheric
heat. We will neglect the very small effects of Raman scattering in Earth’s
atmosphere, where scattered photons emerge at substantially different
frequencies due to internal energy changes of the scattering molecule.
Chandrasekhar[12] uses the symbol $\varpi_{0}$, the first term in his
multipole expansion §3(33) of the scattering phase function, to denote the
single-scattering albedo $\tilde{\omega}$ of (1).

The phase function $p(\mu,\mu^{\prime})$ of (1) is a symmetric and nonnegative
function of the direction cosine, $\mu=\cos\theta$, of the scattered radiation
and the direction cosine, $\mu^{\prime}=\cos\theta^{\prime}$, of incident
radiation

$p(\mu,\mu^{\prime})=p(\mu^{\prime},\mu)\geq 0.$ (3)

The probability $dP$ that a collision with a cloud particulate scatters a
photon, propagating with direction cosine $\mu^{\prime}$, into a photon with a
direction cosine between $\mu$ and $\mu+d\mu$ is
$dP=\tilde{\omega}\,p(\mu,\mu^{\prime})d\mu/2$. In keeping with its
significance as a probability density, the phase function satisfies the
identity

$\frac{1}{2}\int_{-1}^{1}d\mu\,p(\mu,\mu^{\prime})=1.$ (4)

The Planck intensity $B=B(\tau)$ of (1) depends on the local temperature, $T$
of the scattering medium, and on the spatial frequency $\nu$ of the radiation,
as described by

$B=\frac{2h_{\rm P}c^{2}\nu^{3}}{e^{\nu c\,h_{\rm P}/(k_{\rm B}T)}-1}.$ (5)

We are using cgs units, where $h_{\rm P}$ is Planck’s constant and $k_{\rm B}$
is Boltzmann’s constant. We assume that the absolute air temperature, $T=T(z)$
may depend on altitude $z$. Then the Planck intensity, $B=B(\tau)$, may depend
on the altitude $z$ or equivalently, on the optical depth $\tau$. The
radiation wavelength is $\lambda=1/\nu$.

Figure 1: Sample directions of $2n=4$ streams of axially symmetric radiation.
The streams are centered on conical surfaces with opening angles $\theta_{i}$
to the zenith. In accordance with (6), the direction cosines of the streams,
$\mu_{i}=\cos\theta_{i}$, are the zeros of the Legendre polynomial $P_{2n}$,
in this example, $P_{2n}(\mu_{i})=P_{4}(\mu_{i})=0$. Figure 2: Representative
parameters for Gauss-Legendre quadratures and $2n=4$ streams. In accordance
with (6) the sample direction cosines, $\mu_{i}=\cos\theta_{i}$, of the
streams of Fig. 1 are the zeros of the Legendre polynomial,
$P_{2n}(\mu)=P_{4}(\mu)$, which is shown as the continuous red curve. The
zeros $\mu_{i}$, where $P_{2n}(\mu_{i})=0$, are marked with small red circles.
For this example the values are
$\mu_{1},\mu_{2},\mu_{3},\mu_{4}=-0.8611,-0.3400,\,0.3400,\,0.8611$. The
stream weights, calculated with (13), are
$w_{1},w_{2},w_{3},w_{4}=0.3479,\,0.6521,\,0.6521,\,0.3479$. See the text for
more discussion.

###  2.1 The $2n$-stream basis

The $2n$-stream method of reference[1] allows one to solve the integro-
differential equation of transfer (1) accurately and efficiently with modern
computer software. The angular dependence of the intensity is characterized
with $2n$ sample values, $I(\mu_{i},\tau)\geq 0$, along the directions of the
streams $i=1,2,3,\ldots,2n$. As sketched in Fig. 1, the $i$th stream makes an
angle $\theta_{i}=\cos^{-1}\mu_{i}$ to the zenith. The Gauss-Legendre [14]
direction cosines, $\mu_{i}$, are the zeros of the Legendre polynomial
$P_{2n}$ of degree $2n$,

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

We will choose the indices $i=1,2,3,\cdots,2n$ such that

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

Because the Legendre polynomial $P_{2n}$ is even, with
$P_{2n}(\mu)=P_{2n}(-\mu)$, the values of $\mu_{i}$ occur as equal and
opposite pairs,

$\mu_{i}=-\mu_{r(i)}.$ (8)

The index reflection function is

$r(i)=2n+1-i.$ (9)

For the ordering convention (7)

$\displaystyle\mu_{j}$ $\displaystyle<0$ $\displaystyle\quad\hbox{for}\quad
j=1,2,3,\ldots,n,$ (10) $\displaystyle\mu_{k}$ $\displaystyle>0$
$\displaystyle\quad\hbox{for}\quad k=n+1,n+2,n+3,\ldots,2n.$ (11)

The weighted sample values of the intensity, $w_{i}I(\mu_{i},\tau)$, will be
denoted with the symbol $(\kern-2.59995pt(\mu_{i}|I(\tau)\\}$,

$(\kern-2.59995pt(\mu_{i}|I(\tau)\\}=w_{i}I(\mu_{i},\tau)\geq 0.$ (12)

A formula for the Gauss-Legendre weights, $w_{i}>0$, was given by Eq. (72) of
reference [1] as

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

The weights $w_{i}$ sum to 2,

$\displaystyle\sum_{i=1}^{2n}w_{i}=2.$ (14)

The weights $w_{i}$ of (13) and the sample direction cosines $\mu_{i}$ defined
by (6) are illustrated in Fig. 2. The weights are almost, but not precisely,
the distances between local maxima and minima of $P_{2n}(\mu)$ for values of
$\mu$ in the interval $[-1,\,1]$. For example,
$dP_{4}(\mu)/d\mu=(5\mu/3)(7\mu^{2}-3)$. The solutions to $dP_{4}(\mu)/d\mu=0$
are $\mu=0,\pm\sqrt{3/7}$. But $\sqrt{3/7}=0.6547\neq w_{3}=0.6521$, the
number calculated with (13).

To simplify equations, we denote the intensity as an abstract vector
$|I\\}=|I(\tau)\\}$. We can represent the abstract vector $|I\\}$ with a
$2n\times 1$ column vector

$|I\\}=\sum_{i=1}^{2n}|\mu_{i})(\kern-2.59995pt(\mu_{i}|I\\}=\left[\begin{array}[]{c}(\kern-2.59995pt(\mu_{1}|I\\}\\\
(\kern-2.59995pt(\mu_{2}|I\\}\\\ \vdots\\\
(\kern-2.59995pt(\mu_{2n}|I\\}\end{array}\right].$ (15)

We will call the column vector on the right of (15) the $\mu$-space
representation of the abstract vector $|I\\}$. We will use other $2n\times 1$
arrays of numbers to represent the same abstract vector $|I\\}$ in other
bases, for example, the multipole basis $|l)$, discussed in Section 2.2.

The stream basis vectors $|\mu_{i})$ of (15) can be represented with the unit
column vectors

$|\mu_{1})=\left[\begin{array}[]{c}1\\\ 0\\\ \vdots\\\
0\end{array}\right],\quad|\mu_{2})=\left[\begin{array}[]{c}0\\\ 1\\\ \vdots\\\
0\end{array}\right],\quad\cdots,\quad|\mu_{2n})=\left[\begin{array}[]{c}0\\\
0\\\ \vdots\\\ 1\end{array}\right].$ (16)

Corresponding left basis vectors can be represented as the unit row vectors

$\displaystyle(\kern-2.59995pt(\mu_{1}|$ $\displaystyle=$
$\displaystyle[1\quad 0\quad\cdots\quad 0],$
$\displaystyle(\kern-2.59995pt(\mu_{2}|$ $\displaystyle=$
$\displaystyle[0\quad 1\quad\cdots\quad 0],$ $\displaystyle\vdots$
$\displaystyle(\kern-2.59995pt(\mu_{2n}|$ $\displaystyle=$
$\displaystyle[0\quad 0\quad\cdots\quad 1].$

We use a double left parenthesis, $(\kern-2.59995pt(\mu_{i}|$, as a reminder
that the row (or left) basis vectors need not be Hermitian conjugates of the
column (or right) basis vectors. The row basis vectors
$(\kern-2.59995pt(\mu_{i}|$ are like reciprocal lattice vectors of a crystal
[15]. The column basis vectors, $|\mu_{i})$ are like direct lattice vectors.
Just as low-symmetry crystals can have oblique, non-orthogonal lattice
vectors, the right basis vectors $|\mu_{i})$ need not be orthogonal to each
other, although they are orthonormal to the left eigenvectors
$(\kern-2.59995pt(\mu_{i}|$.

As discussed in connection with Eq. (105) of reference [1], the stream basis
vectors are right and left eigenvectors of the direction cosine matrix
$\hat{\mu}$

$\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}.$
(18)

The eigenvectors are chosen to have the orthonormality property

$(\kern-2.59995pt(\mu_{i}|\mu_{i^{\prime}})=\delta_{ii^{\prime}}.$ (19)

They have the completeness property

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

In (20) and elsewhere, we will use the symbol $\hat{1}$ to denote a square
identity matrix with ones along the main diagonal and zeros elsewhere. It has
the same dimensions as any square matrices that are added to, subtracted from,
or equated to it.

Multiplying (20) on the left or right by $\hat{\mu}$ and using (18) we find an
expression for the direction-cosine operator,

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

The direction secant matrix $\hat{\varsigma}$ is the inverse of the direction
cosine matrix $\hat{\mu}$. We can use (21) to write

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

The eigenvalues of the direction secant matrix are the inverses of the
eigenvalues $\mu_{i}$ of the direction cosine matrix

$\varsigma_{i}=\frac{1}{\mu_{i}}.$ (23)

The direction secant matrix $\hat{\varsigma}$ has the same left and right
eigenvectors as the direction cosine matrix $\hat{\mu}$,

$\hat{\varsigma}|\mu_{i})=\varsigma_{i}|\mu_{i}),\quad\hbox{and}\quad(\kern-2.59995pt(\mu_{i}|\hat{\varsigma}=(\kern-2.59995pt(\mu_{i}|\varsigma_{i}.$
(24)

In accordance with Eq. (88) of reference [1] it is convenient to use the
stream basis vectors to define a projection matrix $\mathcal{M}_{\bf d}$ for
downward streams with indices $j\leq n$ and $\mu_{j}<0$, and a projection
matrix $\mathcal{M}_{\bf u}$ for upward streams with indices $k>n$ and
$\mu_{k}>0$.

$\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}|.$ (25)

The projection matrices of (25) have the simple algebra

$\displaystyle\mathcal{M}_{\bf d}+\mathcal{M}_{\bf u}$ $\displaystyle=$
$\displaystyle\hat{1}$ (26) $\displaystyle\mathcal{M}^{2}_{\bf
d}=\mathcal{M}_{\bf d}\quad$ and $\displaystyle\quad\mathcal{M}^{2}_{\bf
u}=\mathcal{M}_{\bf u}$ (27) $\displaystyle\mathcal{M}_{\bf d}\mathcal{M}_{\bf
u}=\breve{0}\quad$ and $\displaystyle\quad\mathcal{M}_{\bf u}\mathcal{M}_{\bf
d}=\breve{0}.$ (28)

Here and elsewhere, the symbol $\breve{0}$ denotes a matrix, not necessarily
square, for which all the elements are zero. The dimensions of $\breve{0}$ are
the same as the dimensions of other matrices to which it is added, subtracted
or equated.

In accordance with Eqs. (109) and (110) of reference [1], we write the
direction-cosine matrix $\hat{\mu}$ as the sum of a downward part
$\hat{\mu}_{\bf d}$ and an upward part $\hat{\mu}_{\bf u}$

$\hat{\mu}=\hat{\mu}_{\bf d}+\hat{\mu}_{\bf u}.$ (29)

Expressions for the downward and upward parts are

$\displaystyle\hat{\mu}_{\bf d}$ $\displaystyle=$
$\displaystyle\mathcal{M}_{\bf d}\hat{\mu}=\hat{\mu}\mathcal{M}_{\bf
d}=\sum_{j=1}^{n}\mu_{j}|\mu_{j})(\kern-2.59995pt(\mu_{j}|,$ (30)
$\displaystyle\hat{\mu}_{\bf u}$ $\displaystyle=$
$\displaystyle\mathcal{M}_{\bf u}\hat{\mu}=\hat{\mu}\mathcal{M}_{\bf
u}=\sum_{k=n+1}^{2n}\mu_{k}|\mu_{k})(\kern-2.59995pt(\mu_{k}|.$ (31)

In like manner, we write the direction-secant matrix $\hat{\varsigma}$ as the
sum of a downward part $\hat{\varsigma}_{\bf d}$ and an upward part
$\hat{\varsigma}_{\bf u}$

$\hat{\varsigma}=\hat{\varsigma}_{\bf d}+\hat{\varsigma}_{\bf u}.$ (32)

Expressions for the downward and upward parts are

$\displaystyle\hat{\varsigma}_{\bf d}$ $\displaystyle=$
$\displaystyle\mathcal{M}_{\bf
d}\hat{\varsigma}=\hat{\varsigma}\mathcal{M}_{\bf
d}=\sum_{j=1}^{n}\varsigma_{j}|\mu_{j})(\kern-2.59995pt(\mu_{j}|,$ (33)
$\displaystyle\hat{\varsigma}_{\bf u}$ $\displaystyle=$
$\displaystyle\mathcal{M}_{\bf
u}\hat{\varsigma}=\hat{\varsigma}\mathcal{M}_{\bf
u}=\sum_{k=n+1}^{2n}\varsigma_{k}|\mu_{k})(\kern-2.59995pt(\mu_{k}|.$ (34)

It will often be convenient to write intensity vectors, flux vectors, source
vectors, etc. as sums of downward and upward parts, for example,

$|I\\}=|I_{\bf d}\\}+|I_{\bf u}\\}=\left[\begin{array}[]{c}|I_{\bf d}\\}\\\
|I_{\bf u}\\}\end{array}\right],$ (35)

where

$|I_{\bf d}\\}=\mathcal{M}_{\bf d}|I\\}\quad\hbox{and}\quad|I_{\bf
u}\\}=\mathcal{M}_{\bf u}|I\\}.$ (36)

###  2.2 The multipole basis

Describing the angular distribution of the axially symmetric intensity
$I(\mu,\tau)$ with the $2n$ sample values, $I(\mu_{i},\tau)$, at the Gauss-
Legendre direction cosines $\mu_{i}$ of (6), is equivalent to approximating
the intensity as a superposition of the first $2n$ Legendre polynomials,

$I(\mu,\tau)=\sum_{l=0}^{2n-1}(2l+1)P_{l}(\mu)I_{l}(\tau).$ (37)

The intensity multipoles are

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

Projections $(\kern-2.59995pt(l|\mu_{i})$ of the left multipole basis
$(\kern-2.59995pt(l|$ onto the right stream basis $|\mu_{i})$, and vice versa,
were given by Eqs. (84) and (85) of reference [1] in terms of Legendre
polynomials $P_{l}$, and weights $w_{i}$ of (13) as

$(\kern-2.59995pt(l|\mu_{i})=\frac{1}{2}P_{l}(\mu_{i}),$ (39)

and

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

Substituting (39) and (12) into (38), and noting from (37) that $I(\mu,\tau)$
can be written as a superpositon of the first $2n$ Legendre polynomials, we
can use the Gauss-Legendre quadrature [14] to write the $l$-th multipole
moment of the intensity as

$\displaystyle I_{l}(\tau)$ $\displaystyle=$
$\displaystyle\frac{1}{2}\sum_{i=1}^{2n}w_{i}P_{l}(\mu_{i})I(\mu_{i},\tau)$
(41) $\displaystyle\to$ $\displaystyle\frac{1}{2}\int_{-1}^{1}d\mu
P_{l}(\mu)I(\mu,\tau)\quad\hbox{as}\quad n\to\infty.$

In analogy to (19) the multipole basis vectors $|l^{\prime})$ and
$(\kern-2.59995pt(l|$ have been chosen to have the orthonormality property

$(\kern-2.59995pt(l|l^{\prime})=\delta_{ll^{\prime}}.$ (42)

In analogy to (20), they have the completeness property

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

Using (42) we can write the intensity vector as

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

The expansion coefficients $(\kern-2.59995pt(l|I(\tau)\\}=I_{l}(\tau)$ were
given by (41).

From (40) we see that the elements $(\kern-2.59995pt(\mu_{i}|0)$ of the right
monopole basis vector $|0)$ are the weights $w_{i}$ of (13)

$|0)=\sum_{i=1}^{2n}|\mu_{i})(\kern-2.59995pt(\mu_{i}|0)=\sum_{i=1}^{2n}|\mu_{i})w_{i}=\left[\begin{array}[]{c}w_{1}\\\
w_{2}\\\ \vdots\\\ w_{2n}\end{array}\right],$ (45)

From (39) we see that the elements $(\kern-2.59995pt(0|\mu_{i})$ of the left
monopole basis vector $(\kern-2.59995pt(0|$ are all equal to 1/2,

$(\kern-2.59995pt(0|=\sum_{i=1}^{2n}(\kern-2.59995pt(0|\mu_{i})(\kern-2.59995pt(\mu_{i}|=\frac{1}{2}\sum_{i=1}^{2n}(\kern-2.59995pt(\mu_{i}|=\bigg{[}\frac{1}{2}\quad\frac{1}{2}\quad\cdots\quad\frac{1}{2}\bigg{]}.$
(46)

An identity from Eq. (53) of reference [3] that will be useful subsequently is

$(\kern-2.59995pt(0|\hat{\mu}=(\kern-2.59995pt(1|\quad\hbox{and}\quad\hat{\mu}|0)=\frac{1}{3}|1).$
(47)

To facilitate subsequent discussions, we note the identity from Eq. (217) of
reference [3],

$(\kern-2.59995pt(0|\hat{\mu}_{\bf u}^{q}e^{-\hat{\varsigma}_{\bf
u}\tau}|0)=(-1)^{q}(\kern-2.59995pt(0|\hat{\mu}_{\bf
d}^{q}e^{\hat{\varsigma}_{\bf d}\tau}|0)=\frac{1}{2}E^{\\{n\\}}_{q+2}(\tau).$
(48)

Here $q$ is an integer, most often $q=0$ or $q=1$ in our work. The $n$-stream
exponential integral functions,

$E^{\\{n\\}}_{q}(\tau)=\sum_{k=n+1}^{2n}w_{k}\mu_{k}^{q-2}e^{-\tau/\mu_{k}},$
(49)

can be obtained by evaluating the exact exponential integral functions,

$E_{q}(\tau)=\int_{0}^{1}d\mu\,\mu^{q-2}e^{-\tau/\mu},$ (50)

with Gauss-Legendre quadratures. The exponential integral functions (50) are
discussed in Appendix I of Chandrasekhar’s book [12]. They account for the
contributions to radiation transfer of intensity propagating at various slant
angles with respect to the vertical. Graphical plots of the functions (49) and
(50) for $n\geq 5$ can hardly be distinguished on a linear scale, as shown by
Fig. 9 of reference [3]. Values of $E_{q}^{\\{n\\}}(0)$ are given in Table 1
of reference [3].

The first few multipole moments $I_{l}$ of the intensity have useful physical
interpretations. As shown by Eq. (19) of reference [1], the monopole moment
($l=0$) is proportional to the volume energy density $u$ of the radiation,

$(\kern-2.59995pt(0|I\\}=I_{0}=\frac{cu}{4\pi},$ (51)

where $c$ is the speed of light. As shown by Eq. (20) of reference [1], the
dipole moment ($l=1$) is proportional to the vertical energy flux $Z$ of the
radiation,

$(\kern-2.59995pt(1|I\\}=I_{1}=\frac{Z}{4\pi}.$ (52)

Under many conditions of practical importance one finds that the intensity is
nearly isotropic and its quadrupole and higher moments are very small compared
to the monopole moment

$(\kern-2.59995pt(l|I\\}=I_{l}\ll I_{0},\quad\hbox{for}\quad l=2,3,4,\ldots$
(53)

When radiation-transfer conditions are such that (53) is valid, it is common
to use the Eddington approximation

$(\kern-2.59995pt(l|I\\}=I_{l}=0\quad\hbox{for}\quad
l=2,3,4,\ldots\quad\hbox{(the Eddington approximation).}$ (54)

As shown in Fig. 1 and Fig. 2 of reference [2], the Eddington approximation
can be very good deep inside an optically thick cloud. But the approximation
is not good just inside the top and bottom surfaces. We do not use the
Eddington approximation in this paper.

###  2.3 Equation of radiative transfer

According to Eq. (62) of reference [1], for a $2n$-stream model of radiative
transfer the integro-differential equation of transfer (1) simplifies to the
first-order, linear differential equation for the intensity vector $|I\\}$,

$\frac{d}{d\tau}|I\\}=\hat{\kappa}\left(|B\\}-|I\\}\right).$ (55)

The Planck intensity vector $|B\\}=|B(\tau)\\}$ of (55) is

$|B\\}=|0)B.$ (56)

Here the right monopole basis vector $|0)$ was given by (45), and the scalar
Planck intensity $B=B(\tau)$ was given by (5). The exponentiation-rate matrix
$\hat{\kappa}$ of (55) was given by Eq. (63) of reference [1] as

$\hat{\kappa}=\hat{\varsigma}\hat{\eta},$ (57)

the product of the direction-secant matrix $\hat{\varsigma}$ of (22), and the
efficiency matrix $\hat{\eta}$. As discussed in Section 3.2.1 of reference
[1], the eigenvalues $\eta_{l}$ of the efficiency matrix are the fraction of
the $l$th multipole moment that remains after each generation of scattering.
As shown in (54) of reference [1], the efficiency matrix can be written as

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

Here the single-scattering albedo $\tilde{\omega}$, which we mentioned in
connection with (1), is the probability that a photon that collides with a
cloud particulate is scattered, rather than being absorbed and converted to
heat. Probabilities must be nonnegative and no larger than 1. So
$\tilde{\omega}$ must be bounded by

$0\leq\tilde{\omega}\leq 1.$ (59)

The continuous phase function $p(\mu,\mu^{\prime})$ of (3) is represented by
the scattering phase matrix $\hat{p}$ in (58). The matrix elements
$(\kern-2.59995pt(\mu_{i}|\hat{p}|\mu_{i^{\prime}})=w_{i}p(\mu_{i},\mu_{i^{\prime}})$
give the representation of $\hat{p}$ in $\mu$-space as a $2n\times 2n$ array
of numbers. The matrix $\hat{p}$ is defined such that a photon in the stream
$i^{\prime}$ that is not absorbed in a collision with a cloud particulate has
a probability $(\kern-2.59995pt(\mu_{i}|\hat{p}|\mu_{i^{\prime}})/2$ to be
scattered into the stream $i$. Therefore, in analogy to (4), we must have

$\sum_{i=1}^{2n}\frac{1}{2}(\kern-2.59995pt(\mu_{i}|\hat{p}|\mu_{i^{\prime}})=1\quad\hbox{and}\quad
0\leq\frac{1}{2}(\kern-2.59995pt(\mu_{i}|\hat{p}|\mu_{i^{\prime}})\leq 1.$
(60)

In accordance with Eq. (40) of reference [1], for a cloud of randomly oriented
scattering particulates or gas molecules the scattering matrix $\hat{p}$ can
be written in terms of the right and left multipole basis vectors, $|l)$ and
$(\kern-2.59995pt(l|$ of (39) and (40), as

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

where the multipole phase coefficients are $p_{l}$. The monopole coefficient
is always unity

$p_{0}=1.$ (62)

The coefficients $p_{l}$ of higher multipolarity, $1\leq l\leq 2n-1$, can be
used to represent any physically permissible (nonnegative) phase function for
a $2n$-stream model, for example: isotropic scattering, Rayleigh scattering,
strongly peaked forward or backward scattering, etc.

For isotropic scattering, the phase matrix is simply

$\displaystyle\hat{p}=2|0)(\kern-2.59995pt(0|.$ (63)

According to Eq. (133) of reference [1], the only non-zero multipole phase
coefficients $p_{l}$ of a Rayleigh-scattering phase function are $p_{0}=1$ and
$p_{2}=1/10$. So the phase matrix of (61) for Rayleigh scattering is

$\displaystyle\hat{p}=2|0)(\kern-2.59995pt(0|+\frac{1}{5}|2)(\kern-2.59995pt(2|.$
(64)

For a $2n$ model, the possible multipole indices are $l=0,1,2,\ldots,2n-1$.
Therefore, to represent the Rayleigh scattering operator (64), which includes
terms with $l=2$, we must have $2n-1\geq 2$ or $n\geq 3/2$. Since $n$ must be
an integer, Raleigh scattering phase operators (64) can only be represented in
models with $n\geq 2$.

As in Section 5.2 of reference [1], we let $p(\mu,1)=\varpi^{\\{p\\}}(\mu)$
denote the phase function $p(\mu,1)$ that can be constructed from the first
$2p$ Legendre polynomials $P_{0}(\mu),P_{1}(\mu),\ldots,P_{2p-1}(\mu)$, and
which gives the maximum possible forward scattering $p(1,1)$, subject to the
constraint that $p(\mu,1)\geq 0$ for any direction cosine, $-1\leq\mu\leq 1$.
The multipole coefficients $p_{l}$ of $\varpi^{\\{p\\}}(\mu)$ are denoted by
$p_{l}=\varpi^{\\{p\\}}_{l}$ and are listed in Table 1 of reference [1]. Then
the phase operator (61) for maximum forward scattering of radiation modeled
with $2n$ streams is

$\hat{p}=2\sum_{l=0}^{2n-1}\varpi^{\\{n\\}}_{l}|l)(\kern-2.59995pt(l|.$ (65)

For an illustative model with $2n=10$, the phase operator (65) for maximum
forward scattering becomes

$\displaystyle\hat{p}$ $\displaystyle=$ $\displaystyle
2\bigg{(}|0)(\kern-2.59995pt(0|+.8182\,|1)(\kern-2.59995pt(1|+.7273\,|2)(\kern-2.59995pt(2|+.5967\,|3)(\kern-2.59995pt(3|+.4988\,|4)(\kern-2.59995pt(4|$
(66)
$\displaystyle+.3869\,|5)(\kern-2.59995pt(5|+.2937\,|6)(\kern-2.59995pt(6|+.2016\,|7)(\kern-2.59995pt(7|+.1209\,|8)(\kern-2.59995pt(8|+.0573\,|9)(\kern-2.59995pt(9|\bigg{)}.$

The phase function for maximum backward scattering differs from (65) for
maximum forward scattering by having alternating signs for the multipole
expansion coefficients,

$\hat{p}=2\sum_{l=0}^{2n-1}(-1)^{l}\varpi^{\\{n\\}}_{l}|l)(\kern-2.59995pt(l|.$
(67)

As shown in Eq. (138) of reference [1], the phase function
$p(\mu,1)=\varpi^{\\{n\\}}(\mu)$ modeled by (66) is strongly peaked in the
forward direction, $\varpi^{\\{n\\}}(1)=n(n+1)$. More detailed discussions of
the scattering-phase matrix $\hat{p}$ and its multipole coefficients $p_{l}$
can be found in Section 5 of reference [1].

###  2.4 Vertical flux

For quantitative studies of vertical energy transfer, the vertical flux vector

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

given by Eq. (210) of reference [1] and the corresponding scalar flux

$\displaystyle Z$ $\displaystyle=$ $\displaystyle(\kern-2.59995pt(0|Z\\}$ (69)
$\displaystyle=$ $\displaystyle 4\pi(\kern-2.59995pt(0|\hat{\mu}|I\\}$
$\displaystyle=$ $\displaystyle 4\pi(1|I\\}$ $\displaystyle=$ $\displaystyle
4\pi I_{1},$

are more directly useful than the intensity vector $|I\\}$ of (15).

According to (12), the elements of the intensity vector are always nonnegative
$(\kern-2.59995pt(\mu_{i}|I\\}\geq 0$. But according to (10) and (11), the
elements $(\kern-2.59995pt(\mu_{i}|Z\\}$ of the vertical flux vector can have
either positive or negative signs, with

$\displaystyle(\kern-2.59995pt(\mu_{j}|Z\\}$ $\displaystyle=$ $\displaystyle
4\pi\mu_{j}(\kern-2.59995pt(\mu_{j}|I\\}\leq 0\quad\hbox{for}\quad
j=1,2,3,\ldots,n,$ (70) $\displaystyle(\kern-2.59995pt(\mu_{k}|Z\\}$
$\displaystyle=$ $\displaystyle 4\pi\mu_{k}(\kern-2.59995pt(\mu_{k}|I\\}\geq
0\quad\hbox{for}\quad k=n+1,n+2,n+3,\ldots,2n.$ (71)

We can use (20), (46), (37), (38), (18) and (12) to write the scalar vertical
flux (69) as

$\displaystyle Z$ $\displaystyle=$ $\displaystyle
4\pi\sum_{i=1}^{2n}(\kern-2.59995pt(0|\mu_{i})(\kern-2.59995pt(\mu_{i}|\hat{\mu}|I\\}$
(72) $\displaystyle=$ $\displaystyle
2\pi\sum_{i=1}^{2n}w_{i}\mu_{i}I(\mu_{i})$ $\displaystyle\to$ $\displaystyle
2\pi\int_{-1}^{1}d\mu\,\mu I(\mu)=\int_{4\pi}d\Omega\mu I(\mu).$

We noted that the second line of (72) is proportional to a Gauss-Legendre
quadrature of the function $2\pi\mu\,I(\mu)$, which converges to the
continuous integral of the third line when $n\to\infty$. The unit of solid
angle for axially symmetric radiation is $d\Omega=2\pi d\mu$.

Using (70) and (71) with the definitions (35) and (36) of downward and upward
parts of radiation vectors, we write the downward and upward parts of the
scalar flux as

$\displaystyle Z_{\bf d}$ $\displaystyle=$
$\displaystyle(\kern-2.59995pt(0|Z_{\bf
d}\\}=(\kern-2.59995pt(0|\mathcal{M}_{\bf d}|Z\\}\leq 0,$ (73) $\displaystyle
Z_{\bf u}$ $\displaystyle=$ $\displaystyle(\kern-2.59995pt(0|Z_{\bf
u}\\}=(\kern-2.59995pt(0|\mathcal{M}_{\bf u}|Z\\}\geq 0,$ (74)

with

$Z=Z_{\bf d}+Z_{\bf u}.$ (75)

###  2.5 Outgoing and incoming radiation

The radiation coming into and going out of a stack of $m$ clouds can be
characterized with the intensity vector $|I^{(m)}\\}$ just above the top of
the stack, and by the intensity vector $|I^{(0)}\\}$ just below the bottom.
Alternatively, one can characterize the intensity with the incoming intensity
$|I^{(\rm in)}\\}$ and outgoing intensity $|I^{(\rm out)}\\}$. As discussed in
Eq. (174) and (175) of reference [1], we can write the incoming intensity
vector as

$\displaystyle|I^{(\rm in)}\\}$ $\displaystyle=$ $\displaystyle|I^{(\rm
in)}_{\bf d}\\}+|I^{(\rm in)}_{\bf u}\\}$ (76) $\displaystyle=$
$\displaystyle|I^{(m)}_{\bf d}\\}+|I^{(0)}_{\bf u}\\}.$

where

$|I^{(\rm in)}_{\bf d}\\}=|I^{(m)}_{\bf d}\\}\quad\hbox{and}\quad|I^{(\rm
in)}_{\bf u}\\}=|I^{(0)}_{\bf u}\\}.$ (77)

The upward and downward parts of the intensity vectors of (76) and (77) are
defined in accordance with (36). In like manner, the outgoing intensity vector
can be written as

$\displaystyle|I^{(\rm out)}\\}$ $\displaystyle=$ $\displaystyle|I^{(\rm
out)}_{\bf d}\\}+|I^{(\rm out)}_{\bf u}\\}$ (78) $\displaystyle=$
$\displaystyle|I^{(0)}_{\bf d}\\}+|I^{(m)}_{\bf u}\\}.$

where

$|I^{(\rm out)}_{\bf d}\\}=|I^{(0)}_{\bf d}\\}\quad\hbox{and}\quad|I^{(\rm
out)}_{\bf u}\\}=|I^{(m)}_{\bf u}\\}.$ (79)

The inverses of (76) and (78) are

$\displaystyle|I^{(0)}\\}$ $\displaystyle=$ $\displaystyle|I^{(0)}_{\bf
d}\\}+|I^{(0)}_{\bf u}\\}$ (80) $\displaystyle=$ $\displaystyle|I^{(\rm
out)}_{\bf d}\\}+|I^{(\rm in)}_{\bf u}\\},$

and

$\displaystyle|I^{(m)}\\}$ $\displaystyle=$ $\displaystyle|I^{(m)}_{\bf
d}\\}+|I^{(m)}_{\bf u}\\}$ (81) $\displaystyle=$ $\displaystyle|I^{(\rm
in)}_{\bf d}\\}+|I^{(\rm out)}_{\bf u}\\}.$

In accordance with (68) we write the flux vectors above and below a cloud as

$\displaystyle|Z^{(0)}\\}$ $\displaystyle=$ $\displaystyle
4\pi\hat{\mu}|I^{(0)}\\},$ (82) $\displaystyle|Z^{(m)}\\}$ $\displaystyle=$
$\displaystyle 4\pi\hat{\mu}|I^{(m)}\\}.$ (83)

Substituting (80) into (82) we find

$|Z^{(0)}\\}=|Z^{(\rm in)}_{\bf u}\\}-|Z^{(\rm out)}_{\bf d}\\},$ (84)

where we define the upward part $|Z^{(\rm in)}_{\bf u}\\}$ of the incoming
flux and the downward part $|Z^{(\rm out)}_{\bf d}\\}$ of the outgoing flux by

$\displaystyle|Z^{(\rm in)}_{\bf u}\\}$ $\displaystyle=$ $\displaystyle
4\pi\mu_{\bf u}|I^{(\rm in)}\\},$ (85) $\displaystyle|Z^{(\rm out)}_{\bf
d}\\}$ $\displaystyle=$ $\displaystyle-4\pi\mu_{\bf d}|I^{(\rm out)}\\}.$ (86)

By convention, we use a negative sign on the right of (86) to ensure that
elements of the vector $|Z^{(\rm out)}_{\bf d}\\}$ are nonnegative in $\mu$
space. Multiplying (84) on the left by $(\kern-2.59995pt(0|$ we find the that
the scalar flux defined by (69) is

$Z^{(0)}=Z^{(\rm in)}_{\bf u}-Z^{(\rm out)}_{\bf d}.$ (87)

Substituting (81) into (83) we find

$|Z^{(m)}\\}=-|Z^{(\rm in)}_{\bf d}\\}+|Z^{(\rm out)}_{\bf u}\\},$ (88)

where

$\displaystyle|Z^{(\rm in)}_{\bf d}\\}$ $\displaystyle=$
$\displaystyle-4\pi\mu_{\bf d}|I^{(\rm in)}\\},$ (89) $\displaystyle|Z^{(\rm
out)}_{\bf u}\\}$ $\displaystyle=$ $\displaystyle 4\pi\mu_{\bf u}|I^{(\rm
out)}\\}.$ (90)

Multiplying (88) on the left by $(\kern-2.59995pt(0|$ we find the that the
scalar flux defined by (69) is

$Z^{(m)}=-Z^{(\rm in)}_{\bf d}+Z^{(\rm out)}_{\bf u}.$ (91)

We can use (84) – (90) to write the outgoing and incoming flux vectors as

$\displaystyle|Z^{(\rm out)}\\}$ $\displaystyle=$ $\displaystyle
4\pi(\hat{\mu}_{\bf u}-\hat{\mu}_{\bf d})|I^{(\rm out)}\\}=4\pi\sum_{\bf
q}s_{\bf q}\hat{\mu}_{\bf q}|I^{(\rm out)}\\}$ (92) $\displaystyle|Z^{(\rm
in)}\\}$ $\displaystyle=$ $\displaystyle 4\pi(\hat{\mu}_{\bf u}-\hat{\mu}_{\bf
d})|I^{(\rm in)}\\}=4\pi\sum_{\bf q}s_{\bf q}\hat{\mu}_{\bf q}|I^{(\rm
in)}\\}$ (93)

Here the sign factors are

$s_{\bf u}=1\quad\hbox{and}\quad s_{\bf d}=-1.$ (94)

The negative signs in the equations (92) and (93) ensure that the elements of
the outgoing and incoming flux vectors are nonnegative,

$(\kern-2.59995pt(\mu_{i}|Z^{(\rm out)}\\}\geq
0\quad\hbox{and}\quad(\kern-2.59995pt(\mu_{i}|Z^{(\rm in)}\\}\geq
0\quad\hbox{for}\quad i=1,2,3,\ldots,2n.$ (95)

From (46) and (95) we see that the scalar fluxes corresponding to (92) and
(93) are nonnegative

$\displaystyle Z^{(\rm out)}$ $\displaystyle=$
$\displaystyle(\kern-2.59995pt(0|Z^{(\rm
out)}\\}=\frac{1}{2}\sum_{i=1}^{2n}(\kern-2.59995pt(\mu_{i}|Z^{(\rm
out)}\\}\geq 0,$ (96) $\displaystyle Z^{(\rm in)}$ $\displaystyle=$
$\displaystyle(\kern-2.59995pt(0|Z^{(\rm
in)}\\}=\frac{1}{2}\sum_{i=1}^{2n}(\kern-2.59995pt(\mu_{i}|Z^{(\rm in)}\\}\geq
0.$ (97)

Subtracting (88) from (84) and using the identity (26) we find

$|Z^{(0)}\\}-|Z^{(m)}\\}=|Z^{(\rm in)}\\}-|Z^{(\rm out)}\\}.$ (98)

Multiplying (98) on the left by $(\kern-2.59995pt(0|$ we find

$Z^{(0)}-Z^{(m)}=Z^{(\rm in)}-Z^{(\rm out)}.$ (99)

According to (99), the difference between the scalar fluxes $Z^{(0)}$ and
$Z^{(m)}$, each of which can be positive, negative or zero, is equal to the
difference between the nonnegative flux $Z^{(\rm in)}\geq 0$ flowing into the
cloud through the top and bottom, and the nonnegative flux $Z^{(\rm out)}\geq
0$ flowing out. Positive values of the flux differences of (99) correspond to
radiative heating of the cloud; negative values correspond to radiative
cooling.

###  2.6 Incident and thermal radiation

As shown in Secton 2.6 of reference [3], the intensity at an optical depth
$\tau$ above the bottom of a cloud can be written as the sum of a part
$|\dot{I}\\}=|\dot{I}(\tau)\\}$ from thermal emission of particulates and gas
molecules inside the cloud and a part $|\ddot{I}\\}=|\ddot{I}(\tau)\\}$ from
the transmission, absorption and scattering of incoming (incident) radiation,

$|I\\}=|\dot{I}\\}+|\ddot{I}\\}.$ (100)

The single dots denote quantities originating from internally generated
thermal radiation. The double dots denote quantities originating from external
incoming radiation. We do not use single and double dot to represent first and
second time derivatives, a convention that goes back to Isaac Newton. Most of
the work of this paper is focussed on steady-state radiation transfer for
which there is no time dependence.

We write the intensity (76) that is incident on the top and bottom of a stack
of $m\geq 1$ clouds, or onto a single isolated cloud, with $m=1$ as

$\displaystyle|I^{(\rm in)}\\}$ $\displaystyle=$ $\displaystyle|\dot{I}^{(\rm
in)}\\}+|\ddot{I}^{(\rm in)}\\}$ (101) $\displaystyle=$
$\displaystyle|\ddot{I}^{(\rm in)}\\}$

Thermal emission of particulates and gas molecules in the clouds can generate
outgoing intensity, $|\dot{I}^{(\rm out)}\\}>\breve{0}$, but it cannot
generate incoming intensity. Therefore,

$|\dot{I}^{(\rm in)}\\}=\breve{0},\quad\hbox{and}\quad|\dot{Z}^{(\rm
in)}\\}=\breve{0}.$ (102)

The second equation of (102) follows from (93).

In accordance with (100) we write the outgoing radiation (78)

$|I^{(\rm out)}\\}=|\dot{I}^{(\rm out)}\\}+|\ddot{I}^{(\rm out)}\\},$ (103)

As shown by Eq. (111) of reference [3], the outgoing intensity vector
$|\ddot{I}^{\\{\rm out\\}}\\}$ is proportional to the incoming intensity
vector $|\ddot{I}^{\\{\rm in\\}}\\}$. The coefficient of proportionality is
the scattering matrix $\mathcal{S}$,

$|\ddot{I}^{(\rm out)}\\}=\mathcal{S}|\ddot{I}^{(\rm in)}\\}.$ (104)

We will frequently write the scattering matrix of a single cloud as the
$2\times 2$ block matrix

$\mathcal{S}=\left[\begin{array}[]{ll}\mathcal{S}_{\bf dd}&\mathcal{S}_{\bf
du}\\\ \mathcal{S}_{\bf ud}&\mathcal{S}_{\bf
uu}\end{array}\right]\quad\hbox{where}\quad\mathcal{S}_{\bf
qq^{\prime}}=\mathcal{M}_{\bf q}\mathcal{S}\mathcal{M}_{\bf q^{\prime}}.$
(105)

Possible values of the stream-direction indices ${\bf q}$ and ${\bf
q^{\prime}}$ of (105) are ${\bf d}$ and ${\bf u}$. In Section 2.8 we review
how to calculate the scattering matrix $\mathcal{S}$ for a homogeneous cloud.

Using (93) and (104) we write (92) as

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

In accordance with Eq. (216) of reference [1], the cloud albedo matrix
$\Omega$ of (106) is a similarity transformation of the scattering matrix,

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

The part of the outgoing intensity (103) that comes from thermal emission of
cloud particulates and gas molecules can be written as

$|\dot{I}^{(\rm
out)}\\}=|\dot{J}\\}=\int_{0}^{\tau}d\tau^{\prime}G(\tau^{\prime})|0)B(\tau^{\prime}).$
(108)

In (108) $\tau$ is the vertical optical thickness of the cloud,
$G(\tau^{\prime})$ is the $2n\times 2n$ continuous Green’s matrix for
thermally generated outgoing radiation, the monopole basis vector of (45) is
$|0)$, and $B(\tau^{\prime})$, given by (5), is the Planck intensity at the
source optical depth $\tau^{\prime}$ above the bottom of the cloud. We review
how to construct $G(\tau^{\prime})$ for homogeneous clouds in Section 2.9. The
analogous discrete Green’s matrix $G^{[mc\\}}$ of (261) depends on the
discrete index $c=1,2,3,\ldots,m$ of individual clouds in a stack of $m$
clouds, rather than on the continuous optical depth $\tau^{\prime}$ above the
bottom of a single cloud.

For the special case of a single isothermal cloud of constant Planck intensity
$B(\tau^{\prime})=B$, (108) simplifies to

$|\dot{I}^{(\rm out)}\\}=|\dot{J}\\}=\mathcal{E}|0)B.$ (109)

The emissivity matrix $\mathcal{E}$ of the isothermal cloud is the integral of
the continuous Green’s matrix over the entire optical depth $\tau$ of the
cloud,

$\mathcal{E}=\int_{0}^{\tau}d\tau^{\prime}G(\tau^{\prime}).$ (110)

According to Kirchhoff’s law of radiation, Eq. (279) of reference [1], the
emissivity matrix $\mathcal{E}$ of (110) is related to the scattering matrix
$\mathcal{S}$ by

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

An explicit proof of (111) for homogeneous clouds is given by (139) below. In
analogy to (105) we will frequently write the emissivity matrix of an
isothermal cloud as the $2\times 2$ block matrix

$\mathcal{E}=\left[\begin{array}[]{ll}\mathcal{E}_{\bf dd}&\mathcal{E}_{\bf
du}\\\ \mathcal{E}_{\bf ud}&\mathcal{E}_{\bf
uu}\end{array}\right]\quad\hbox{where}\quad\mathcal{E}_{\bf
qq^{\prime}}=\mathcal{M}_{\bf q}\mathcal{E}\mathcal{M}_{\bf q^{\prime}}.$
(112)

In summary, for a single cloud we can write (103) as

$\displaystyle|I^{(\rm out)}\\}=|\dot{J}\\}+\mathcal{S}|\ddot{I}^{(\rm
in)}\\}.$ (113)

The output intensity vector $|I^{(\rm out)}\\}$ of (113) is a linear
combination of intensity $|\dot{J}\\}$, thermally emitted by cloud
particulates and gas molecules, and transmitted and scattered incoming
intensity $\mathcal{S}|\ddot{I}^{(\rm in)}\\}$. For long wave thermal
radiation, both cloud particulates and greenhouse gas molecules absorb and
emit radiation, but only only cloud particulates contribute significantly to
scattering. Short wave sunlight is efficiently scattered by cloud
particulates. Especially for blue and ultraviolet sunlight, there is also
significant Rayleigh scattering by atmospheric gases. Both particulates and
gases absorb small fractions of sunlight, but Earth’s atmosphere is normally
too cool to emit short wave radiation.

###  2.7 Identities for scattering matrices $\mathcal{S}$

Some important fundamental properties of scattering matrices $\mathcal{S}$ of
(104) are summarized here. Formal proofs will be given in a subsequent paper.

For conservatively scattering clouds, with unit single-scattering albedo,
$\tilde{\omega}=1$, no thermal radiation can be emitted. That is, we must have

$|\dot{I}^{(\rm out)}\\}=\breve{0}\quad\hbox{when}\quad\tilde{\omega}=1.$
(114)

Comparing (114) with (109) we see that for a conservatively scattering cloud,
with $\tilde{\omega}=1$,

$\mathcal{E}|0)=\breve{0}.$ (115)

Using (115) with Kirchhoff’s law (111) we find for a conservatively scattering
cloud

$\mathcal{S}|0)=|0).$ (116)

Eq. (116) is the conservative-scattering limit of the more general inequality

$0\leq(\kern-2.59995pt(\mu_{i}|\mathcal{S}|0)\leq(\kern-2.59995pt(\mu_{i}|0)=w_{i}\quad\hbox{for}\quad
0\leq\tilde{\omega}\leq 1.$ (117)

The elements of any scattering matrix, for a homogeneous or inhomogeneous
cloud, must satisfy the Helmholtz-reciprocity symmetry

$\frac{(\kern-2.59995pt(\mu_{o}|\mathcal{S}|\mu_{i})}{w_{o}|\mu_{i}|}=\frac{(\kern-2.59995pt(\mu_{r(i)}|\mathcal{S}|\mu_{r(o)})}{w_{i}|\mu_{o}|}.$
(118)

The index reflection function $r(i)$ was defined by (9). The reciprocity
theorem (118) quantifies the plausible fact that for clouds described by (55),
the rate of scattering from an input stream with index $i$ to an output stream
with index $o$ is proportional to the rate of scattering of the time-reversed
output stream with index $r(o)$ into the time-reversed input stream with index
$r(i)$.

The elements of scattering matrices of homogeneous clouds also satisfy the
simpler reflection symmetry

$(\kern-2.59995pt(\mu_{o}|\mathcal{S}|\mu_{i})=(\kern-2.59995pt(\mu_{r(o)}|\mathcal{S}|\mu_{r(i)}).$
(119)

Scattering matrix elements of more general and realistic inhomogeneous clouds
do not satisfy the reflection symmetry (119) but they do satisfy the Helmholtz
reciprocity symmetry (118).

Helmholtz reciprocity symmetries were discussed by Chandrasekhar [12] in his
Sections §13 and §52. Chandrasekhar uses a normalization of the scattering
matrix $\mathcal{S}$ that eliminates the factors $|\mu_{i}|$ and $|\mu_{o}|$
from (118). We use a slightly different normalization to simplify the forms of
equations like (116), (117) and (119). So formulas involving the $\mathcal{S}$
matrix are slightly different in our work from those in Chandrasekhar’s book.

###  2.8 Scattering matrix $\mathcal{S}$ for a homogeneous cloud

According to Eq. (206) of reference [1], the scattering matrix $\mathcal{S}$
of a homogeneous cloud is given by the simple formula

$\mathcal{S}=\mathcal{O}\mathcal{I}^{-1}.$ (120)

Here we review how to calculate the outgoing matrix $\mathcal{O}$ and incoming
matrix $\mathcal{I}$ of (120). Let $(\kern-2.59995pt(\lambda_{i}|$ and
$|\lambda_{i})$ be the left and right eigenvectors, and let $\lambda_{i}$ be
the corresponding eigenvalue of the penetration-length matrix,

$\hat{\lambda}=\hat{\kappa}^{-1}=\sum_{i=1}^{2n}\lambda_{i}|\lambda_{i})(\kern-2.59995pt(\lambda_{i}|,$
(121)

the inverse of the exponentiation rate matrix $\hat{\kappa}$ of (57). As in
(7) the real eigenvalues or penetration lengths are ordered such that

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

As in (8) the penetration lengths for a homogenous cloud have the reflection
symmetry

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

The index reflection function $r(i)$ was defined by (9). The $\lambda$-space
basis vectors $(\kern-2.59995pt(\lambda_{j}|$ and $|\lambda_{i})$ are chosen
to have orthonormality and completeness relations analogous to (20) and (21),

$(\kern-2.59995pt(\lambda_{i}|\lambda_{j})=\delta_{ij},$ (124)

and

$\displaystyle\hat{1}$ $\displaystyle=$
$\displaystyle\sum_{i=1}^{2n}|\lambda_{i})(\kern-2.59995pt(\lambda_{i}|$ (125)
$\displaystyle=$ $\displaystyle\mathcal{L}_{\bf d}+\mathcal{L}_{\bf u}.$

For the limit of vanishing single scattering albedos, $\tilde{\omega}\to 0$,
$\lambda$-space quantities are chosen to approach the corresponding
$\mu$-space quantities.

$\left[\begin{array}[]{c}\lambda_{i}\\\ (\kern-2.59995pt(\lambda_{i}|\\\
|\lambda_{i})\end{array}\right]\to\left[\begin{array}[]{c}\mu_{i}\\\
(\kern-2.59995pt(\mu_{i}|\\\
|\mu_{i})\end{array}\right]\quad\hbox{as}\quad\tilde{\omega}\to 0.$ (126)

In analogy to (25) we define the downward and upward projection operators in
$\lambda$ space by

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

The exponentiation-rate operator (57) can be written as the sum of downward
and upward parts in $\lambda$ space

$\hat{\kappa}=\hat{\kappa}_{\bf d}+\hat{\kappa}_{\bf u}.$ (128)

where

$\hat{\kappa}_{\bf
d}=\sum_{j=1}^{n}\kappa_{j}|\lambda_{j})(\kern-2.59995pt(\lambda_{j}|\quad\hbox{and}\quad\hat{\kappa}_{\bf
u}=\sum_{k=n+1}^{2n}\kappa_{k}|\lambda_{k})(\kern-2.59995pt(\lambda_{k}|.$
(129)

The eigenvalues $\kappa_{i}$ of $\hat{\kappa}$ are the inverses of the
eigenvalues $\lambda_{i}$ of $\hat{\lambda}$, and the eigenvectors are the
same as those of $\hat{\lambda}$,

$\kappa_{i}=\frac{1}{\lambda_{i}},\quad(\kern-2.59995pt(\kappa_{i}|=(\kern-2.59995pt(\lambda_{i}|,\quad\hbox{and}\quad|\kappa_{i})=|\lambda_{i}).$
(130)

The overlap matrix $\mathcal{C}$ between $\mu$ space and $\lambda$ space is
defined by Eq. (162) of reference [1] as

$\mathcal{C}=\left[\begin{array}[]{ll}\mathcal{C}_{\bf dd}&\mathcal{C}_{\bf
du}\\\ \mathcal{C}_{\bf ud}&\mathcal{C}_{\bf
uu}\end{array}\right]=\left[\begin{array}[]{ll}\mathcal{M}_{\bf
d}\mathcal{L}_{\bf d}&\mathcal{M}_{\bf d}\mathcal{L}_{\bf u}\\\
\mathcal{M}_{\bf u}\mathcal{L}_{\bf d}&\mathcal{M}_{\bf u}\mathcal{L}_{\bf
u}\end{array}\right].$ (131)

Note that the left directional index, ${\bf q=d}$ or ${\bf u}$, of
$\mathcal{C}_{\bf qq^{\prime}}$ refers to $\mu$ space, and the right
directional index, ${\bf q^{\prime}=d}$ or ${\bf u}$, refers to $\lambda$
space. A useful identity for the overlap matrices was already given as Eq.
(163) of reference [1]

$\displaystyle\mathcal{C}_{\bf dd}+\mathcal{C}_{\bf ud}+\mathcal{C}_{\bf
du}+\mathcal{C}_{\bf uu}$ $\displaystyle=$ $\displaystyle\mathcal{M}_{\bf
d}(\mathcal{L}_{\bf d}+\mathcal{L}_{\bf u})+\mathcal{M}_{\bf
u}(\mathcal{L}_{\bf d}+\mathcal{L}_{\bf u})$ (132) $\displaystyle=$
$\displaystyle\mathcal{M}_{\bf d}+\mathcal{M}_{\bf u}$ $\displaystyle=$
$\displaystyle\hat{1}.$

Here we noted the sum rules (26) and (125) of the upward and downward
projection operators $\mathcal{M}_{\bf q}$ in $\mu$ space and
$\mathcal{L}_{\bf q}$ in $\lambda$ space. Each of the matrices
$\mathcal{C}_{\bf qq^{\prime}}$ can be represented as a $2n\times 2n$ matrix
so it does not matter if we write the $\mathcal{C}_{\bf qq^{\prime}}$ as an
element of a $2\times 2$ block matrix or simply add them as in (132).

A slight modification of the overlap matrix gives the incoming matrix, defined
by Eq. (199) of reference [1] as

$\mathcal{I}=\left[\begin{array}[]{ll}\mathcal{C}_{\bf dd}&\mathcal{C}_{\bf
du}e^{-\hat{\kappa}_{\bf u}\tau}\\\ \mathcal{C}_{\bf ud}e^{\hat{\kappa}_{\bf
d}\tau}&\mathcal{C}_{\bf uu}\end{array}\right].$ (133)

The outgoing matrix is defined by Eq. (201) of reference [1] by

$\mathcal{O}=\left[\begin{array}[]{ll}\mathcal{C}_{\bf dd}e^{\hat{\kappa}_{\bf
d}\tau}&\mathcal{C}_{\bf du}\\\ \mathcal{C}_{\bf ud}&\mathcal{C}_{\bf
uu}e^{-\hat{\kappa}_{\bf u}\tau}\end{array}\right].$ (134)

A simple limiting case of (120) is a purely absorbing cloud with
$\tilde{\omega}\to 0$ and negligible scattering, like a layer of clear air
with greenhouse-gas absorption. Then $\mathcal{L}_{\bf q}\to\mathcal{M}_{\bf
q}$, $\hat{\kappa}\to\hat{\varsigma}$ and

$\mathcal{I}\to\hat{1},\quad\mathcal{O}\to\left[\begin{array}[]{ll}e^{\hat{\varsigma}_{\bf
d}\tau}&\breve{0}\\\ \breve{0}&e^{-\hat{\varsigma}_{\bf
u}\tau}\end{array}\right],\quad\hbox{and}\quad\mathcal{S}=\mathcal{O}\mathcal{I}^{-1}\to\left[\begin{array}[]{ll}e^{\hat{\varsigma}_{\bf
d}\tau}&\breve{0}\\\ \breve{0}&e^{-\hat{\varsigma}_{\bf
u}\tau}\end{array}\right].$ (135)

###  2.9 Continuous Green’s matrix $G(\tau^{\prime})$

In reference [3] we used a continuous Green’s vector
$|G(\tau^{\prime})\\}=G(\tau^{\prime})|0)$ instead of the continuous Green’s
matrix $G(\tau^{\prime})$ of (110). In Eq. (151) of that reference we showed
that $G(\tau^{\prime})$ for a homogeneous cloud can be written as

$\displaystyle G(\tau^{\prime})$ $\displaystyle=$
$\displaystyle\left[\mathcal{Q}(\tau^{\prime})+\mathcal{S}\mathcal{R}(\tau^{\prime})\right]\hat{\kappa}$
(136) $\displaystyle=$
$\displaystyle\frac{\partial}{\partial\tau^{\prime}}\left[\mathcal{Q}(\tau^{\prime})+\mathcal{S}\mathcal{R}(\tau^{\prime})\right].$

In (136) the scattering matrix $S=S(\tau)$ depends on the total optical
thickness $\tau$ of the cloud, but is independent of the variable optical
thickness $\tau^{\prime}$ above the bottom of the cloud. The retro matrix
$\mathcal{R}(\tau^{\prime})$ was given by Eq. (127) of reference [3] as

$\displaystyle\mathcal{R}(\tau^{\prime})$ $\displaystyle=$
$\displaystyle\mathcal{C}_{\bf ud}e^{\hat{\kappa}_{\bf
d}\tau^{\prime}}-\mathcal{C}_{\bf du}e^{-\hat{\kappa}_{\bf
u}(\tau-\tau^{\prime})}.$ (137)

The matrix $\mathcal{Q}(\tau^{\prime})$ was given by Eq. (149) of reference
[3] as

$\displaystyle\mathcal{Q}(\tau^{\prime})$ $\displaystyle=$
$\displaystyle\mathcal{C}_{\bf uu}e^{-\hat{\kappa}_{\bf
u}(\tau-\tau^{\prime})}-\mathcal{C}_{\bf dd}e^{\hat{\kappa}_{\bf
d}\tau^{\prime}}.$ (138)

To verify (110) we integrate (136) to find

$\displaystyle\int_{0}^{\tau}d\tau^{\prime}G(\tau^{\prime})$ $\displaystyle=$
$\displaystyle\mathcal{Q}(\tau)-\mathcal{Q}(0)+\mathcal{S}\left[\mathcal{R}(\tau)-\mathcal{R}(0)\right]$
(139) $\displaystyle=$ $\displaystyle\mathcal{C}_{\bf uu}-\mathcal{C}_{\bf
dd}e^{\hat{\kappa}_{\bf d}\tau}-\mathcal{C}_{\bf uu}e^{-\hat{\kappa}_{\bf
u}\tau}+\mathcal{C}_{\bf dd}$ $\displaystyle+\mathcal{S}\left[\mathcal{C}_{\bf
ud}e^{\hat{\kappa}_{\bf d}\tau}-\mathcal{C}_{\bf du}-\mathcal{C}_{\bf
ud}+\mathcal{C}_{\bf du}e^{-\hat{\kappa}_{\bf u}\tau}\right]$ $\displaystyle=$
$\displaystyle\hat{1}-\mathcal{O}+\mathcal{S}\left[\mathcal{I}-\hat{1}\right]$
$\displaystyle=$ $\displaystyle\hat{1}-\mathcal{S}$ $\displaystyle=$
$\displaystyle\mathcal{E}.$

To derive the expression to the right of the third equal sign of (139) from
the previous line, we used (132), (134) and (133). To derive fourth line we
used (120). The last line follows from (111). We will call (139) the Kirchhoff
identity. A discretized version (139) for cloud stacks is given by (276).

###  2.10 Black clouds

A single black cloud absorbs all incident intensity and scatters none so the
scattering matrix $\mathcal{S}^{\\{1\\}}$ is the null matrix,

$\mathcal{S}=\breve{0}.$ (140)

A purely absorbing cloud with the scattering matrix (135) becomes black as the
optical thickness approaches infinity, $\tau\to\infty$. From Kirchhoff’s law
(111) and from (140) we see that the emissivity matrix of a black cloud is the
identity matrix

$\mathcal{E}=\hat{1}.$ (141)

Then according to (109) and (141) the intensity emitted by an isothermal black
cloud with the Planck intensity $B$ is blackbody radiation

$|\dot{I}^{(\rm out)}\\}=|\dot{J}\\}=|0)B.$ (142)

The output flux vector (92) corresponding to (142) is

$|\dot{Z}^{(\rm out)}\\}=4\pi(\hat{\mu}_{\bf u}-\hat{\mu}_{\bf d})|0)B.$ (143)

Using (143) we write the scalar output flux (69) of a black cloud as

$\displaystyle\dot{Z}^{(\rm out)}=(\kern-2.59995pt(0|\dot{Z}^{(\rm
out)}\\}=2\dot{Z}^{(\rm bb)}.$ (144)

Here the scalar blackbody flux $\dot{Z}^{(\rm bb)}$, emitted upward from the
top of the black cloud, or downward from the bottom, is

$\displaystyle\dot{Z}^{(\rm bb)}$ $\displaystyle=$ $\displaystyle
4\pi(\kern-2.59995pt(0|\hat{\mu}_{\bf u}|0)B$ (145) $\displaystyle=$
$\displaystyle-4\pi(\kern-2.59995pt(0|\hat{\mu}_{\bf d}|0)B$ $\displaystyle=$
$\displaystyle 2\pi E^{\\{n\\}}_{3}(0)\,B$ $\displaystyle\to$
$\displaystyle\pi B\quad\hbox{as}\quad n\to\infty.$

To write the third line of (145) we used the expresssion (48) for the
$n$-stream analog $E^{\\{n\\}}_{3}(\tau)$ of the exponential integral function
$E_{3}(\tau)$. To write the last line of (145) we recalled from (222) of
reference [3] that $E^{\\{n\\}}_{3}(0)$ converges to $1/2$ as $n\to\infty$.
The convergence is rapid. For $n=5$, Table 1 of reference [3] gives
$E^{\\{5\\}}_{3}(0)=0.5038$.

The frequency-integrated value of the flux (145) is

$\displaystyle\overline{\dot{Z}^{(\rm bb)}}$ $\displaystyle=$
$\displaystyle\int_{0}^{\infty}d\nu\,\dot{Z}^{(\rm bb)}$ (146)
$\displaystyle=$ $\displaystyle 2E^{\\{n\\}}_{3}(0)\,\sigma_{\rm SB}T^{4}$
$\displaystyle\to$ $\displaystyle\sigma_{\rm SB}T^{4}\quad\hbox{as}\quad
n\to\infty.$

Here we noted that the integral of the Planck intensity (5) over all upward
solid angle increments, $2\pi\mu d\mu$, and frequency increments, $d\nu$, is

$2\pi\int_{0}^{1}d\mu\,\mu\int_{0}^{\infty}d\nu B=\pi\int_{0}^{\infty}d\nu
B=\sigma_{\rm SB}T^{4}.$ (147)

The Stefan-Boltzmann constant is

$\displaystyle\sigma_{\rm SB}$ $\displaystyle=$ $\displaystyle\frac{2\pi
k_{\rm B}^{4}}{c^{2}h_{\rm P}^{3}}\int_{0}^{\infty}dx\frac{x^{3}}{e^{x}-1}$
(148) $\displaystyle=$ $\displaystyle\frac{2\pi^{5}k_{\rm
B}^{4}}{15c^{2}h_{\rm P}^{3}}.$

According to (144), an isothermal black cloud emits equal upward and downward
fluxes $\dot{Z}^{(\rm bb)}$ through its top and bottom surfaces. The
frequency-integrated fluxes (146) are very nearly $\overline{\dot{Z}^{(\rm
bb)}}=\sigma_{\rm SB}T^{4}$.

###  2.11 Planck emissivities

For a non-black, isothermal cloud with Planck intensity $B$, we can use (109)
with (92) to write the outgoing thermal flux vector as

$\displaystyle|\dot{Z}^{(\rm out)}\\}$ $\displaystyle=$ $\displaystyle
4\pi(\hat{\mu}_{\bf u}-\hat{\mu}_{\bf d})|\dot{I}^{(\rm out)}\\}$ (149)
$\displaystyle=$ $\displaystyle 4\pi(\hat{\mu}_{\bf u}-\hat{\mu}_{\bf
d})\mathcal{E}|0)B.$

In accordance with (82) the scalar outgoing flux that corresponds to (149) is

$\dot{Z}^{(\rm out)}=(\kern-2.59995pt(0|\dot{Z}^{(\rm out)}\\}=\dot{Z}^{(\rm
out)}_{\bf u}+\dot{Z}^{(\rm out)}_{\bf d}.$ (150)

The upward flux from the top of the cloud is

$\dot{Z}^{(\rm out)}_{\bf u}=4\pi(\kern-2.59995pt(0|\hat{\mu}_{\bf
u}\mathcal{E}|0)B=\varepsilon_{\bf u}\dot{Z}^{(\rm bb)},$ (151)

where the blackbody flux $\dot{Z}^{(\rm bb)}$ was given by (145). The Planck
emissivity $\varepsilon_{\bf u}$ for the cloud top is

$\varepsilon_{\bf u}=\frac{\dot{Z}^{(\rm out)}_{\bf u}}{\dot{Z}^{(\rm
bb)}}=\frac{(\kern-2.59995pt(0|\hat{\mu}_{\bf
u}\mathcal{E}|0)}{(\kern-2.59995pt(0|\hat{\mu}_{\bf
u}|0)}=1-\frac{(\kern-2.59995pt(0|\hat{\mu}_{\bf
u}\mathcal{S}|0)}{(\kern-2.59995pt(0|\hat{\mu}_{\bf u}|0)}.$ (152)

In like manner, we write the downward flux from the bottom of the cloud as

$\dot{Z}^{(\rm out)}_{\bf d}=-4\pi(\kern-2.59995pt(0|\hat{\mu}_{\bf
d}\mathcal{E}|0)B^{\\{1\\}}=\varepsilon_{\bf d}\dot{Z}^{(\rm bb)}.$ (153)

The Planck emissivity $\varepsilon_{\bf d}$ for the cloud bottom is

$\varepsilon_{\bf d}=\frac{\dot{Z}^{(\rm out)}_{\bf d}}{\dot{Z}^{(\rm
bb)}}=\frac{(\kern-2.59995pt(0|\hat{\mu}_{\bf
d}\mathcal{E}|0)}{(\kern-2.59995pt(0|\hat{\mu}_{\bf
d}|0)}=1-\frac{(\kern-2.59995pt(0|\hat{\mu}_{\bf
d}\mathcal{S}|0)}{(\kern-2.59995pt(0|\hat{\mu}_{\bf d}|0)}.$ (154)

The total flux from the top and bottom of an isothermal cloud is

$\dot{Z}^{(\rm out)}=(\varepsilon_{\bf u}+\varepsilon_{\bf d})\dot{Z}^{(\rm
bb)},$ (155)

From (117), (152) and (154) we see that the emissivities are bounded by

$0\leq\varepsilon_{\bf u}\leq 1\quad\hbox{and}\quad 0\leq\varepsilon_{\bf
d}\leq 1$ (156)

By symmetry, a homogeneous isothermal cloud emits equal upward and downward
scalar fluxes, so $\varepsilon_{\bf u}=\varepsilon_{\bf d}$. For an
inhomogeneous cloud, $\varepsilon_{\bf u}$ and $\varepsilon_{\bf d}$ may not
be equal because particulates near the top of the cloud may have different
single-scattering albedos $\tilde{\omega}$ and scattering-phase matrices
$\hat{p}$ from those near the bottom.

Figure 3: Planck emissivities $\varepsilon_{\bf u}$ of (152) for homogeneous,
optically thick clouds with the scattering-phase matrices $\hat{p}$ of (63) –
(67) as functions of the single-scattering albedos $\tilde{\omega}$. See the
text for more detail.

Examples of Planck emissivities $\varepsilon_{\bf u}(\tau)$ for homogeneous
clouds of finite thickness $\tau$ were shown in Fig. 10 of reference [3],
where the expression (152) for the emissivity was given as Eq. (229). The
emissivities increase with increasing optical depths $\tau$ of the clouds. In
the limit $\tau\to\infty$, the emissivity saturates at a value,
$\varepsilon_{\bf u}(\infty)\leq 1$, that depends on the single-scattering
albedo $\tilde{\omega}$ and the scattering-phase matrix $\hat{p}$ of (61). For
purely absorbing clouds with $\tilde{\omega}=0$, the optically-thick
emissivity is $\varepsilon_{\bf u}(\infty)=1$.

Some representative examples of the Planck emissivities $\varepsilon_{\bf u}$
of homogeneous, optically thick clouds, calculated with (152) as functions of
the single-scattering albedo $\tilde{\omega}$ are shown in Fig. 3. The curves
are for fixed values of the four scattering-phase matrices $\hat{p}$ of (63)
–(67). Independent of $\hat{p}$, for black clouds with no scattering or
transmission and only absorption ($\tilde{\omega}=0$) we have
$\varepsilon_{\bf u}=1$. For clouds with 100% scattering or transmission and
no absorption ($\tilde{\omega}=1$) we have $\varepsilon_{\bf u}=0$. Only for
intermediate values of the single scattering albedo, $0<\tilde{\omega}<1$,
does the emissivity depend on the scattering-phase matrix $\hat{p}$ of (61).

###  2.12 Planck albedos

Suppose that the cloud discussed in the preceding section is illuminated from
above and below with isotropic radiation of Planck intensity $B$, so the
incoming intensity is

$|\ddot{I}^{(\rm in)}\\}=|0)B$ (157)

According to (104), the cloud will scatter the incoming intensity into
outgoing intensity

$|\ddot{I}^{(\rm out)}\\}=\mathcal{S}|\ddot{I}^{(\rm in)}\\}=\mathcal{S}|0)B.$
(158)

The scattering matrix $\mathcal{S}$ of the cloud is related to the emissivity
matrix $\mathcal{E}$ in accordance with Kirchhoff’s law (111). We use (82) to
write the scattered scalar flux as

$\ddot{Z}^{(\rm out)}=4\pi(\kern-2.59995pt(0|(\hat{\mu}_{\bf u}-\hat{\mu}_{\bf
d})|\ddot{I}^{(\rm out)}\\}=\ddot{Z}^{(\rm out)}_{\bf u}+\ddot{Z}^{(\rm
out)}_{\bf d}.$ (159)

We can write the upward part of (159) as

$\ddot{Z}^{(\rm out)}_{\bf u}=4\pi(\kern-2.59995pt(0|\hat{\mu}_{\bf
u}\mathcal{S}|0)B=\omega_{\bf u}\dot{Z}^{(\rm bb)}.$ (160)

where we can use the expression (145) for the blackbody flux $\dot{Z}^{(\rm
bb)}$, and Kirchhoff’s law (111) to write the upward part of the Planck albedo
$\omega_{\bf u}$ as

$\omega_{\bf u}=\frac{\ddot{Z}^{(\rm out)}_{\bf u}}{\dot{Z}^{(\rm
bb)}}=\frac{(\kern-2.59995pt(0|\hat{\mu}_{\bf
u}\mathcal{S}|0)}{(\kern-2.59995pt(0|\hat{\mu}_{\bf u}|0)}=1-\varepsilon_{\bf
u}.$ (161)

In like manner we can write the downward part of (159) as

$\ddot{Z}^{(\rm out)}_{\bf d}=-4\pi(\kern-2.59995pt(0|\hat{\mu}_{\bf
d}\mathcal{S}|0)B=\omega_{\bf d}\dot{Z}^{(\rm bb)}.$ (162)

where

$\omega_{\bf d}=\frac{\ddot{Z}^{(\rm out)}_{\bf d}}{\dot{Z}^{(\rm
bb)}}=\frac{(\kern-2.59995pt(0|\hat{\mu}_{\bf
d}\mathcal{S}|0)}{(\kern-2.59995pt(0|\hat{\mu}_{\bf u}|0)}=1-\varepsilon_{\bf
d}.$ (163)

The Planck albedos $\omega_{\bf u}$ and $\omega_{\bf d}$ of (161) and (163)
should not be confused with the single-scattering albedo $\tilde{\omega}$ of a
cloud particulate or gas molecule, which we distinguish with a tilde. A cloud
with purely absorbing particulates and gas molecules, and therefore with a
vanishing single-scattering albedo, $\tilde{\omega}=0$, can have a Planck
albedo $\omega_{\bf q}$ close to 1 if the cloud is optically thin enough that
most of the incoming radiation can be transmitted through the cloud without
absorption. As one would intuitively expect, and as follows formally from
(117) the Planck albedos are bounded by

$0\leq\omega_{\bf u}\leq 1\quad\hbox{and}\quad 0\leq\omega_{\bf d}\leq 1$
(164)

The isotropic input intensity (157) corresponds to the input flux

$\displaystyle\ddot{Z}^{(\rm in)}$ $\displaystyle=$ $\displaystyle
4\pi(\kern-2.59995pt(0|(\hat{\mu}_{\bf u}-\hat{\mu}_{\bf d})|\ddot{I}^{(\rm
in)}\\}$ (165) $\displaystyle=$ $\displaystyle
4\pi(\kern-2.59995pt(0|(\hat{\mu}_{\bf u}-\hat{\mu}_{\bf d})|0)B$
$\displaystyle=$ $\displaystyle 2\dot{Z}^{(\rm bb)},$

where the blackbody flux $\dot{Z}^{(\rm bb)}$ was given by (145). The ratio of
the output flux $\ddot{Z}^{(\rm out)}$ of (159) to the input flux
$\ddot{Z}^{(\rm in)}$ of (165) is therefore

$\frac{\ddot{Z}^{(\rm out)}}{\ddot{Z}^{(\rm in)}}=\frac{\omega_{\bf
u}+\omega_{\bf d}}{2}.$ (166)

The mean Planck albedo, $(\omega_{\bf u}+\omega_{\bf d})/2$ of (166), is much
like the Bond albedo of a planet [16], the fraction of the intercepted solar
energy flux that the planet transmits or scatters back to space without
absorption. But the mean Planck albedo $(\omega_{\bf u}+\omega_{\bf d})/2$ is
the fraction of isotropic, monochromatic incoming radiation per unit area of
the top and bottom of a cloud that is reflected or transmitted, rather than
being absorbed and converted to heat. The Bond albedo for planets in our solar
system is defined for nearly collimated illumination of the entire planet by
the full frequency spectrum of the Sun.

Using (87) and (91) we write the scalar flux $Z(\tau^{\prime})$ as
$\dot{Z}^{(0)}=\dot{Z}(0)$ at the bottom of the cloud, where $\tau^{\prime}=0$
and $\dot{Z}^{(1)}=\dot{Z}(\tau)$ at the top, where $\tau^{\prime}=\tau$.

$\displaystyle\dot{Z}^{(0)}$ $\displaystyle=$ $\displaystyle\dot{Z}^{(\rm
in)}_{\bf u}-\dot{Z}^{(\rm out)}_{\bf d}=-\dot{Z}^{(\rm out)}_{\bf
d}=4\pi(\kern-2.59995pt(0|\hat{\mu}_{\bf d}\mathcal{E}|0)B$ (167)
$\displaystyle\dot{Z}^{(1)}$ $\displaystyle=$ $\displaystyle-\dot{Z}^{(\rm
in)}_{\bf d}+\dot{Z}^{(\rm out)}_{\bf u}=\dot{Z}^{(\rm out)}_{\bf
u}=4\pi(\kern-2.59995pt(0|\hat{\mu}_{\bf u}\mathcal{E}|0)B.$ (168)

Here we noted that a single cloud has no incoming thermal flux,
$|\dot{Z}^{(\rm in)}\\}=\breve{0}$.

To facilitate subsequent discussions of stacks of more than one cloud, we
write (167) and (168) as the single vector equation

$|\dot{Z})=4\pi\dot{M}|B]\quad\hbox{or}\quad
Z^{(g)}=4\pi\sum_{c=1}^{1}\dot{M}^{(gc\\}}B^{\\{c\\}},$ (169)

where $|B]=B^{\\{c\\}}=B^{\\{1\\}}$ is the Planck intensity of cloud $c=1$ in
a stack of $m=1$ clouds, and $Z^{(g)}$ is the flux in the space $g=c$ above
the cloud $c$ or in the space $g=c-1$ below the cloud c. We will call the
spaces below or above clouds gaps. For a stack of $m$ clouds the cloud indices
$c$ and gap indices $g$ can take on the values

$c=1,2,3,\ldots,m\quad\hbox{and}\quad g=0,1,2,3,\ldots m.$ (170)

On the left of (169) the abstract vector and matrix symbols mean

$|\dot{Z})=\left[\begin{array}[]{c}\dot{Z}^{(0)}\\\
\dot{Z}^{(1)}\end{array}\right]\quad\hbox{and}\quad\dot{M}=\left[\begin{array}[]{c}\dot{M}^{(01\\}}\\\
\dot{M}^{(11\\}}\end{array}\right]=\left[\begin{array}[]{c}(\kern-2.59995pt(0|\hat{\mu}_{\bf
d}\mathcal{E}|0)\\\ (\kern-2.59995pt(0|\hat{\mu}_{\bf
u}\mathcal{E}|0)\end{array}\right].$ (171)

Since the abstract matrix $M$ converts a cloud property $|B]$ to a gap
quantity $|\dot{Z})$ we write its elements as $M^{(gc\\}}$ with a left
parenthesis and a right curly bracket as delimiters for the gap and cloud
indices.

###  2.13 Half isotropic incoming intensity

To fully specify the intensity vector $|\ddot{I}^{(\rm in)}\\}$ of incoming
radiation onto the bottom and top of a cloud stack, one needs $2n$ numbers,
for example, the stream amplitudes $(\kern-2.59995pt(\mu_{i}|\ddot{I}^{(\rm
in)}\\}$ for $i=1,2,3,\ldots,2n$. But for this expository work it is
convenient to model input radiation as half-isotropic Planck radiation, which
would be generated by an external black cloud with a Planck intensity
$B^{\\{\rm in\\}}_{\bf d}$, located above the cloud stack, and a second
external black cloud of Planck intensity $B^{\\{\rm in\\}}_{\bf u}$, located
below the cloud stack. The two nonnegative numbers, $B^{\\{\rm in\\}}_{\bf d}$
and $B^{\\{\rm in\\}}_{\bf u}$, are sufficient to specify half-isotropic
incoming radiation with the incoming intensity vector

$|\ddot{I}^{(\rm in)}\\}=\sum_{\bf q=d}^{\bf u}\mathcal{M}_{\bf q}|0)B^{\\{\rm
in\\}}_{\bf q}.$ (172)

The input flux vector (93) corresponding to (172) is

$|\ddot{Z}^{(\rm in)}\\}=4\pi\sum_{\bf q=d}^{\bf u}s_{\bf q}\hat{\mu}_{\bf
q}B^{\\{\rm in\\}}_{\bf q}|0).$ (173)

The sign factors $s_{\bf q}$ were defined by (94). The upward and downward
scalar incoming fluxes corresponding to (173) are

$\displaystyle\ddot{Z}^{(\rm in)}_{\bf u}$ $\displaystyle=$
$\displaystyle(\kern-2.59995pt(0|\mathcal{M}_{\bf u}|\ddot{Z}^{(\rm
in)}\\}=4\pi(\kern-2.59995pt(0|\hat{\mu}_{\bf u}|0)B^{\\{\rm in\\}}_{\bf u},$
(174) $\displaystyle\ddot{Z}^{(\rm in)}_{\bf d}$ $\displaystyle=$
$\displaystyle(\kern-2.59995pt(0|\mathcal{M}_{\bf d}|\ddot{Z}^{(\rm
in)}\\}=-4\pi(\kern-2.59995pt(0|\hat{\mu}_{\bf d}|0)B^{\\{\rm in\\}}_{\bf d}.$
(175)

From (106) we see that the outgoing flux corresponding to the incoming flux
(173) is

$|\ddot{Z}^{(\rm out)}\\}=4\pi\sum_{\bf q=d}^{\bf u}\Omega\,s_{\bf
q}\,\hat{\mu}_{\bf q}|0)B^{\\{\rm in\\}}_{\bf q}.$ (176)

The cloud albedo matrix $\Omega$ in (176) was given by (107). The upward and
downward scalar outgoing fluxes corresponding to the vector flux (176) are

$\displaystyle\ddot{Z}^{(\rm out)}_{\bf u}$ $\displaystyle=$
$\displaystyle(\kern-2.59995pt(0|\mathcal{M}_{\bf u}|\ddot{Z}^{(\rm
out)}\\}=4\pi\sum_{\bf q}(\kern-2.59995pt(0|\mathcal{M}_{\bf u}\Omega\,s_{\bf
q}\,\hat{\mu}_{\bf q}|0)B^{\\{\rm in\\}}_{\bf q},$ (177)
$\displaystyle\ddot{Z}^{(\rm out)}_{\bf d}$ $\displaystyle=$
$\displaystyle(\kern-2.59995pt(0|\mathcal{M}_{\bf d}|\ddot{Z}^{(\rm
out)}\\}=4\pi\sum_{\bf q}(\kern-2.59995pt(0|\mathcal{M}_{\bf d}\Omega\,s_{\bf
q}\,\hat{\mu}_{\bf q}|0)B^{\\{\rm in\\}}_{\bf q},$ (178)

We can use (174), (175), (177) and (178) with (87) to write the scalar flux
$\ddot{Z}^{(0)}$ below the cloud as a linear combination of the input Planck
intensities $B^{\\{\rm in\\}}_{\bf d}$ and $B^{\\{\rm in\\}}_{\bf u}$,

$\displaystyle\ddot{Z}^{(0)}$ $\displaystyle=$ $\displaystyle\ddot{Z}^{(\rm
in)}_{\bf u}-\ddot{Z}^{(\rm out)}_{\bf d}$ (179) $\displaystyle=$
$\displaystyle 4\pi(\kern-2.59995pt(0|\hat{\mu}_{\bf u}-\mathcal{M}_{\bf
d}\Omega\hat{\mu}_{\bf u}|0)B^{\\{\rm in\\}}_{\bf
u}+4\pi(\kern-2.59995pt(0|\mathcal{M}_{\bf d}\Omega\hat{\mu}_{\bf
d}|0)B^{\\{\rm in\\}}_{\bf d}$ $\displaystyle=$
$\displaystyle\ddot{N}^{(0)}_{\bf u}B^{\\{\rm in\\}}_{\bf
u}+\ddot{N}^{(0)}_{\bf d}B^{\\{\rm in\\}}_{\bf d}.$

In like manner, can use (174), (175), (177) and (178) with (91) to write the
scalar flux $\ddot{Z}^{(1)}$ above the cloud as

$\displaystyle\ddot{Z}^{(1)}$ $\displaystyle=$ $\displaystyle-\ddot{Z}^{(\rm
in)}_{\bf d}+\ddot{Z}^{(\rm out)}_{\bf u}$ (180) $\displaystyle=$
$\displaystyle 4\pi(\kern-2.59995pt(0|\hat{\mu}_{\bf d}-\mathcal{M}_{\bf
u}\Omega\hat{\mu}_{\bf d}|0)B^{\\{\rm in\\}}_{\bf
d}+4\pi(\kern-2.59995pt(0|\mathcal{M}_{\bf u}\Omega\hat{\mu}_{\bf
u}|0)B^{\\{\rm in\\}}_{\bf u}$ $\displaystyle=$
$\displaystyle\ddot{N}^{(1)}_{\bf d}B^{\\{\rm in\\}}_{\bf
d}+\ddot{N}^{(1)}_{\bf u}B^{\\{\rm in\\}}_{\bf u}.$

To facilitate subsequent discussions of stacks of more than one cloud, we
write (179) and (180) as the single vector equation

$|\ddot{Z})=4\pi\ddot{N}|B^{\\{\rm
in\\}}\rangle\quad\hbox{or}\quad\ddot{Z}^{(g)}=4\pi\sum_{\bf q=d}^{\bf
u}\ddot{N}^{(g)}_{\bf q}B^{\\{\rm in\\}}_{\bf q},$ (181)

where

$|\ddot{Z})=\left[\begin{array}[]{c}\ddot{Z}^{(0)}\\\
\ddot{Z}^{(1)}\end{array}\right],$ (182)

and

$|B^{\\{\rm in\\}}\rangle=\left[\begin{array}[]{c}B^{\\{\rm in\\}}_{\bf d}\\\
B^{\\{\rm in\\}}_{\bf u}\end{array}\right].$ (183)

The elements of the $(m+1)\times 2=2\times 2$ matrix $\ddot{N}$ of (181)
follow from (179) and (180) and are

$\ddot{N}=\left[\begin{array}[]{cc}\ddot{N}^{(0)}_{\bf d}&\ddot{N}^{(0)}_{\bf
u}\\\ \ddot{N}^{(1)}_{\bf d}&\ddot{N}^{(1)}_{\bf
u}\end{array}\right]=\left[\begin{array}[]{lr}(\kern-2.59995pt(0|\mathcal{M}_{\bf
d}\Omega\hat{\mu}_{\bf d}|0)&(\kern-2.59995pt(0|\hat{\mu}_{\bf
u}-\mathcal{M}_{\bf d}\Omega\hat{\mu}_{\bf u}|0)\\\
(\kern-2.59995pt(0|\hat{\mu}_{\bf d}-\mathcal{M}_{\bf u}\Omega\hat{\mu}_{\bf
d}|0)&(\kern-2.59995pt(0|\mathcal{M}_{\bf u}\Omega\hat{\mu}_{\bf
u}|0)\end{array}\right].$ (184)

Summing the flux $|\dot{Z})$ of (169) from thermal emission and the flux
$|\ddot{Z})$ of (181) from scattering of incoming radiation we find the total
fluxes

$\displaystyle|Z)$ $\displaystyle=$
$\displaystyle\left[\begin{array}[]{c}Z^{(0)}\\\
Z^{(1)}\end{array}\right]=|\dot{Z})+|\ddot{Z})=4\pi\dot{M}|B]+4\pi\ddot{N}|B^{\\{\rm
in\\}}\rangle.$ (187)

###  2.14 Radiative heating and cooling of an isolated cloud

It is instructive to consider a 1-cloud stack consisting of a single isolated
cloud. The number of clouds in the stack is $m=1$ and the cloud index can only
be $c=1$. Using (187) we write the net radiative absorption rate $R^{\\{c\\}}$
per unit area of the isolated cloud as the difference between the vertical
flux $Z^{(0)}$ at the bottom of the cloud and the vertical flux $Z^{(1)}$ at
the top,

$\displaystyle R^{\\{c\\}}$ $\displaystyle=$ $\displaystyle Z^{(0)}-Z^{(1)}$
(188) $\displaystyle=$ $\displaystyle-\sum_{g=0}^{1}\Delta^{\\{cg)}Z^{(g)}.$

The elements of the $m\times(m+1)=1\times 2$ differencing matrix $\Delta$ are

$\displaystyle\Delta$ $\displaystyle=$
$\displaystyle\left[\begin{array}[]{rr}\Delta^{\\{10)}&\Delta^{\\{11)}\end{array}\right]$
(190) $\displaystyle=$
$\displaystyle\left[\begin{array}[]{rr}-1&1\end{array}\right].$ (192)

In preparation for discussions of stacks of more than one cloud we have
written (188) as the element of an abstract vector equation

$\displaystyle|R]$ $\displaystyle=$ $\displaystyle-\Delta|Z)$ (193)
$\displaystyle=$
$\displaystyle-4\pi\Delta\dot{M}|B]-4\pi\Delta\ddot{N}|B^{\\{\rm
in\\}}\rangle.$

The rate $R^{\\{c\\}}$ of (188) is the diabatic heating rate (or cooling rate
if $R^{\\{c\\}}<0$) due to absorption and emission of radiation by cloud
particulates and gas molecules. One definition of diabatic heating is: A
process that occurs with the addition or loss of heat. The opposite of
adiabatic. Meteorological examples include air parcels warming due to the
absorption of infrared radiation or release of latent heat [17]. This
definition refers to the enthalpy the non-condensible (non-water) molecules of
an air parcel. If one were to include the enthalpy of water and water vapor in
the total enthalpy of the air parcel, condensation or evaporation would
generate no diabatic heating or cooling. For example, the heat (enthalpy)
added by condensation of water vapor to the dry air of an expanding air parcel
would be equal and opposite to the heat (enthalpy) lost from the condensing
vapor.

We can write the net heating rate of (188) as the difference between the
heating rate $\ddot{H}^{\\{c\\}}$ due to absorption of external radiation and
the cooling rate $\dot{C}^{\\{c\\}}$ due to thermal emission by cloud
particulates and gas molecules.

$\displaystyle R^{\\{c\\}}=\ddot{H}^{\\{c\\}}-\dot{C}^{\\{c\\}}.$ (194)

From inspection of (188) we see that

$\displaystyle\ddot{H}^{\\{c\\}}$ $\displaystyle=$
$\displaystyle\ddot{Z}^{(0)}-\ddot{Z}^{(1)}$ (195) $\displaystyle=$
$\displaystyle\ddot{Z}^{(\rm in)}-\ddot{Z}^{(\rm out)}$ $\displaystyle=$
$\displaystyle(\kern-2.59995pt(0|\ddot{Z}^{(\rm
in)}\\}-(\kern-2.59995pt(0|\ddot{Z}^{(\rm out)}\\}$ $\displaystyle=$
$\displaystyle(\kern-2.59995pt(0|(\hat{1}-\Omega^{\\{c\\}})|\ddot{Z}^{(\rm
in)}\\}$ $\displaystyle=$
$\displaystyle(\kern-2.59995pt(0|\mathcal{A}^{\\{c\\}}|\ddot{Z}^{(\rm
in)}\\}.$

We used (99) to write the second line of (195) and we used (106) to write the
fourth line. In the last line we have introduced the absorptivity matrix
$\mathcal{A}^{\\{c\\}}$ which we define as the complement of the albedo matrix
$\Omega^{\\{c\\}}$,

$\displaystyle\mathcal{A}^{\\{c\\}}$ $\displaystyle=$
$\displaystyle\hat{1}-\Omega^{\\{c\\}}$ (196) $\displaystyle=$
$\displaystyle(\hat{\mu}_{\bf u}-\hat{\mu}_{\bf
d})\mathcal{E}^{\\{c\\}}(\hat{\varsigma}_{\bf u}-\hat{\varsigma}_{\bf d}).$

The second line of (196) comes from (107) and (111). From (196) we see that
the same similarity transformation (107) that converts the scattering matrix
$\mathcal{S}^{\\{c\\}}$ to the albedo matrix $\Omega^{\\{c\\}}$, converts the
emissivity matrix $\mathcal{E}^{\\{c\\}}$ to the absorptivity matrix
$\mathcal{A}^{\\{c\\}}$. For future reference, we use (173) to write the last
line of (195) as

$\displaystyle\ddot{H}^{\\{c\\}}$ $\displaystyle=$
$\displaystyle(\kern-2.59995pt(0|\mathcal{A}^{\\{c\\}}|\ddot{Z}^{(\rm in)}\\}$
(197) $\displaystyle=$
$\displaystyle(\kern-2.59995pt(0|\mathcal{A}^{\\{c\\}}\left(\hat{\mu}_{\bf
u}B^{(\rm in)}_{\bf u}-\hat{\mu}_{\bf d}B^{(\rm in)}_{\bf d}\right)|0)$

From inspection of (188) we see that the cooling rate due to thermal emission
by particulates and gas molecules of a single isothermal cloud of Planck
intensity $B^{\\{c\\}}$ is

$\displaystyle\dot{C}^{\\{c\\}}$ $\displaystyle=$
$\displaystyle\dot{Z}^{(1)}-\dot{Z}^{(0)}$ (198) $\displaystyle=$
$\displaystyle\Delta|\dot{Z})$ $\displaystyle=$ $\displaystyle
4\pi\Delta\dot{M}B^{\\{c\\}}$ $\displaystyle=$ $\displaystyle
4\pi(\kern-2.59995pt(0|(\hat{\mu}_{\bf u}-\hat{\mu}_{\bf
d})\mathcal{E}^{\\{c\\}}|0)B^{\\{c\\}}$ $\displaystyle=$
$\displaystyle(\varepsilon_{\bf u}+\varepsilon_{\bf d})\dot{Z}^{\\{\rm
bb\\}}.$

We used the $1\times 2$ matrix $\Delta$ of (192) and the $2\times 1$ vector
$|\dot{Z})$ of (171) to write the second line of (198); we used the expression
(192) for $\Delta$ and the expression (169) for $|\dot{Z})$ to write the third
line; we used the matrix elements of (171) to write the fourth line; for the
fifth line we used the definitions (152) and (154) of the Planck emissivities
$\varepsilon_{\bf u}$ and $\varepsilon_{\bf d}$, along with the definition
(145) of the blackbody flux $\dot{Z}^{\\{\rm bb\\}}$. Not surprisingly, the
cooling rate (198) of a single, isolated, isothermal cloud is the sum of the
thermal flux $\varepsilon_{\bf u}\dot{Z}^{\\{\rm bb\\}}$ emitted upward from
the cloud top and the thermal flux $\varepsilon_{\bf d}\dot{Z}^{\\{\rm bb\\}}$
emitted downward from the cloud bottom.

We can use the second line of (196) to write the fourth line of (198) as the
isolated cloud cooling rate

$\dot{C}^{\\{c\\}}=4\pi(\kern-2.59995pt(0|\mathcal{A}^{\\{c\\}}(\hat{\mu}_{\bf
u}-\hat{\mu}_{\bf d})|0)B^{\\{c\\}}.$ (199)

Substituting (197) and (199) into (194) we find that the net heating rate of
the cloud is

$\displaystyle R^{\\{c\\}}$ $\displaystyle=$
$\displaystyle\ddot{H}^{\\{c\\}}-\dot{C}^{\\{c\\}}$ (200) $\displaystyle=$
$\displaystyle 4\pi(\kern-2.59995pt(0|\mathcal{A}^{\\{c\\}}([B^{\\{\rm
in\\}}_{\bf u}-B^{\\{c\\}}]\hat{\mu}_{\bf u}-[B^{\\{\rm in\\}}_{\bf
d}-B^{\\{c\\}}]\hat{\mu}_{\bf d})|0).$

###  2.15 Radiative equilibrium

The net radiative heating rate for an isolated cloud will vanish,
$R^{\\{c\\}}=0$. Then the cloud thermally emits just as much radiation out of
its top and bottom surfaces as it absorbs from incoming radiation. Suppose
that convection and other heat transfer mechanism are negligibly small
compared to radiative heat transfer. Then for fixed half-isotropic incoming
radiation, described by the $2\times 1$ vector of upward and downward Planck
intensities, $|B^{(\rm in)}\rangle$ of (183), the temperature $T^{\\{c\\}}$ of
the cloud will rise if $R^{\\{c\\}}>0$ or drop if $R^{\\{c\\}}<0$ until the
Planck intensity $B^{\\{c\\}}$ of the isothermal cloud makes the expression
(200) equal to zero. Setting $R^{\\{c\\}}=0$ in (200) we see that for given
values, $B^{\\{\rm in\\}}_{\bf d}$ and $B^{\\{\rm in\\}}_{\bf u}$, of half-
isotropic incoming radiation the heating rate will vanish if the Planck
intensity of a isolated, isothermal cloud is

$B^{\\{c\\}}=\frac{(\kern-2.59995pt(0|\mathcal{A}^{\\{c\\}}(B^{\\{\rm
in\\}}_{\bf u}\hat{\mu}_{\bf u}-B^{\\{\rm in\\}}_{\bf d}\hat{\mu}_{\bf
d})|0)}{(\kern-2.59995pt(0|\mathcal{A}^{\\{c\\}}(\hat{\mu}_{\bf
u}-\hat{\mu}_{\bf d})|0)}.$ (201)

The cloud Planck intensity $B^{\\{c\\}}$ for radiative equilibrium will be
somewhere between $B^{\\{\rm in\\}}_{\bf d}$ and $B^{\\{\rm in\\}}_{\bf u}$,
and the cloud temperature $T^{\\{c\\}}$ will be somewhere between the
temperatures $T^{\\{\rm in\\}}_{\bf d}$ and $T^{\\{\rm in\\}}_{\bf u}$ of the
downward and upward incoming radiation.

The fundamental symmetries of a homogeneous cloud, with the same single-
scattering albedo $\tilde{\omega}$ and the same scattering-phase matrix
$\hat{p}$ of (61) from top to bottom, are such that

$(\kern-2.59995pt(0|\mathcal{A}^{\\{c\\}}\hat{\mu}_{\bf
u}|0)=-(\kern-2.59995pt(0|\mathcal{A}^{\\{c\\}}\hat{\mu}_{\bf d}|0).$ (202)

Then (201) simplifies to

$B^{\\{1\\}}=\frac{B^{\\{\rm in\\}}_{\bf u}+B^{\\{\rm in\\}}_{\bf d}}{2}.$
(203)

In radiative equilibrium the Planck intensity $B^{\\{1\\}}$ of a homogeneous
cloud is the average of the Planck intensities $B^{\\{\rm in\\}}_{\bf u}$ and
$B^{\\{\rm in\\}}_{\bf d}$ of the half isotropic incoming radiation.

###  2.16 1-cloud examples

To graph numerical results it is convenient to introduce a reference Planck
intensity $B_{0}$ and a reference flux $Z_{0}$. In accordance with (145),
these are related as the flux of a blackbody is related to its Planck
intensity

$\displaystyle Z_{0}$ $\displaystyle=$ $\displaystyle
4\pi(\kern-2.59995pt(0|\hat{\mu}_{\bf u}|0)B_{0}$ (204) $\displaystyle=$
$\displaystyle 2\pi E^{\\{n\\}}_{3}(0)B_{0}$ $\displaystyle\to$
$\displaystyle\pi B_{0}\quad\hbox{as}\quad n\to\infty.$

We use caligraphic fonts to denote fluxes measured in units of $Z_{0}$ or
Planck intensities measured in units of $B_{0}$. For example, for half-
isotropic incoming radiation incident on a homogeneous cloud we can use (204)
and (173) to write

$\displaystyle\ddot{\mathcal{Z}}^{(\rm in)}_{\bf q}=\frac{\ddot{Z}^{(\rm
in)}_{\bf q}}{Z_{0}}=\frac{B^{\\{\rm in\\}}_{\bf
q}}{B_{0}}=\mathcal{B}^{\\{\rm in\\}}_{\bf q}.$ (205)

Figure 4:  A Rayleigh scattering cloud of optical depth $\tau=1$ with the same
temperature as incoming radiation from above and below. The downward and
upward relative fluxes of incoming radiation are denoted by
$\ddot{\mathcal{Z}}^{(\rm in)}_{\bf d}$ and $\ddot{\mathcal{Z}}^{(\rm
in)}_{\bf u}$, and are shown as the dotted cyan lines. The downward and upward
parts of transmitted and reflected radiation are denoted by
$\ddot{\mathcal{Z}}^{(\rm out)}_{\bf d}$ and $\ddot{\mathcal{Z}}^{(\rm
out)}_{\bf u}$, and are shown as the dashed cyan lines. The net fluxes above
and below the cloud from the incoming radiation are denoted by
$\ddot{\mathcal{Z}}^{(1)}$ and $\ddot{\mathcal{Z}}^{(0)}$, and are shown as
the continuous cyan lines. For thermal emission by cloud particulates and gas
molecules, the fluxes below and above the clouds are denoted by
$\dot{\mathcal{Z}}^{(0)}$ and $\dot{\mathcal{Z}}^{(1)}$, and are shown as the
continuous red lines. The net fluxes are denoted by $\mathcal{Z}^{(0)}$ and
$\mathcal{Z}^{(1)}$, and are shown as the dashed green lines. Fluxes with only
upward or only downward radiation are distinguished with arrowheads from
fluxes which are the net of both upward and downward radiation. See the text
for more detailed discussions.  Figure 5: The same cloud as for Fig. 4, but
with only upward incoming radiation incident on the bottom of the cloud. See
the text for more detailed discussions.  Figure 6: The same incoming radiation
as for Fig. 5, but after the cloud has cooled to radiative equilibrium. See
the text for more detailed discussions.

In Fig. 4 – 6 we show numerical examples of radiation transport by a single
homogeneous cloud of optical thickness

$\displaystyle\tau=1.$ (206)

For the examples considered here, we assume that both the particulates and gas
molecules have a Rayleigh scattering-phase matrix $\hat{p}$ of (64) with a
single-scattering albedo

$\displaystyle\tilde{\omega}=\frac{1}{2}.$ (207)

We can use the value (207) for the single-scattering albedo $\tilde{\omega}$,
together with the Rayleigh-scattering-phase matrix $\hat{p}$ of (64) to
evaluate the efficiency matrix $\hat{\eta}$ of (58). Multiplying $\hat{\eta}$
on the left with the secant matrix $\hat{\varsigma}$ of (32) gives the
exponentiation matrix $\hat{\kappa}$ of (57). As summarized in connection with
(105), $\hat{\kappa}$, together with the optical depth $\tau$ of (206), can be
used to evaluate the scattering matrix $\mathcal{S}$ of the cloud. The closely
related albedo matrix $\Omega$ follows from (107), the emissivity matrix
$\mathcal{E}$ follows from (111) and the absorptivity matrix $\mathcal{A}$
follows from (196). We use $2n=10$ streams so $\mathcal{S}$, $\mathcal{E}$,
$\Omega$ and $\mathcal{A}$ are $10\times 10$ matrices, too large for
convenient display in this paper, but readily used by modern mathematical
software packages like Matlab. Using $\mathcal{E}$ to evaluate the upward
Planck emissivity $\varepsilon_{\bf u}$ of (152), and using $\mathcal{S}$ to
evaluate the upward Planck albedo $\omega_{\bf u}$ of (161), we find

$\displaystyle\varepsilon_{\bf u}=0.5568\quad\hbox{and}\quad\omega_{\bf
u}=0.4432.$ (208)

Doubling the number of streams to $2n=20$ decreases the Planck albedo by about
$0.15\%$, from $\omega_{\bf u}=0.4432$ to $\omega_{\bf u}=0.4426$, too small
to be displayed in Fig. 4. So $2n=10$ streams gives good modeling accuracy.

Using $\mathcal{E}$ to evaluate the matrix $\dot{M}$ of (171) we find

$4\pi\dot{M}=\left[\begin{array}[]{r}-0.5568\\\
0.5568\end{array}\right]\frac{Z_{0}}{B_{0}}.$ (209)

Here $Z_{0}$ is the reference flux and $B_{0}$ is the reference Planck
intensity of (204). Using $\Omega$ to evaluate the matrix $\ddot{N}$ of (184)
we find

$4\pi\ddot{N}=\left[\begin{array}[]{rr}-0.3058&0.8627\\\
-0.8627&0.3058\end{array}\right]\frac{Z_{0}}{B_{0}}.$ (210)

Let the Planck intensity $B^{\\{c\\}}$ of the cloud and the Planck intensities
$B^{\\{\rm in\\}}_{\bf d}$ and $B^{\\{\rm in\\}}_{\bf u}$ of the half-
isotropic incoming radiation be

$B^{\\{c\\}}=B_{0}\quad\hbox{and}\quad|B^{\\{\rm
in\\}}\rangle=\left[\begin{array}[]{c}B^{\\{\rm in\\}}_{\bf d}\\\ B^{\\{\rm
in\\}}_{\bf u}\end{array}\right]=\left[\begin{array}[]{c}1\\\
1\end{array}\right]B_{0}.$ (211)

Evaluating (169) with $B^{\\{c\\}}$ from (211) and $\dot{M}$ from (209) we
find

$|\dot{Z})=\left[\begin{array}[]{c}\dot{Z}^{(0)}\\\
\dot{Z}^{(1)}\end{array}\right]=4\pi\dot{M}B^{\\{c\\}}=\left[\begin{array}[]{r}-0.5568\\\
0.5568\end{array}\right]Z_{0}.$ (212)

Evaluating (181) with $|B^{\\{\rm in\\}}\rangle$ from (211) and $\ddot{N}$
from (210) we find (after adjusting the last numerical digit for roundoff
error)

$|\ddot{Z})=\left[\begin{array}[]{c}\ddot{Z}^{(0)}\\\
\ddot{Z}^{(1)}\end{array}\right]=4\pi\ddot{N}|B^{\\{\rm
in\\}}\rangle=\left[\begin{array}[]{r}0.5568\\\
-0.5568\end{array}\right]Z_{0}.$ (213)

We see that the flux (212) due to thermal emission of cloud particulates and
gas molecules is equal and opposite the flux (213) due to transmission and
scattering of incoming radiation. Summing (212) and (213) we find that the
total flux vanishes,

$|Z)=\left[\begin{array}[]{c}Z^{(0)}\\\
Z^{(1)}\end{array}\right]=|\dot{Z})+|\ddot{Z})=\left[\begin{array}[]{r}0\\\
0\end{array}\right].$ (214)

The cloud heating rate (188) also vanishes

$R^{\\{c\\}}=Z^{(0)}-Z^{(1)}=0.$ (215)

The cloud of Fig. 4 is in full thermal equilibium with the incoming radiation,
and it neither heats nor cools.

In Fig. 5 we show what happens if the downward incoming radiation is removed
but everything else remains the same. The Planck intensity $|B^{\\{\rm
in\\}}\rangle$ of (211) is changed to

$|B^{\\{\rm in\\}}\rangle=\left[\begin{array}[]{c}B^{\\{\rm in\\}}_{\bf d}\\\
B^{\\{\rm in\\}}_{\bf u}\end{array}\right]=\left[\begin{array}[]{c}0\\\
1\end{array}\right]B_{0}.$ (216)

Since the cloud retains the same Planck intensity $B^{\\{1\\}}=B_{0}$ as for
(211), the thermally emitted flux $|\dot{Z})$ stays the same as for (212).
Evaluating (181) with $|B^{\\{\rm in\\}}\rangle$ from (216) and $\ddot{N}$
from (210) we find values different from (213)

$|\ddot{Z})=\left[\begin{array}[]{c}\ddot{Z}^{(0)}\\\
\ddot{Z}^{(1)}\end{array}\right]=\left[\begin{array}[]{r}0.8627\\\
0.3058\end{array}\right]Z_{0}.$ (217)

Summing (212) and (217) we find that the net flux is positive both below and
above the cloud

$|Z)=\left[\begin{array}[]{c}Z^{(0)}\\\
Z^{(1)}\end{array}\right]=|\dot{Z})+|\ddot{Z})=\left[\begin{array}[]{r}0.3058\\\
0.8627\end{array}\right]Z_{0}.$ (218)

The cloud heating rate (188) is negative, that is, the cloud cools by
releasing more vertical flux $Z^{(1)}=0.8627\,Z_{0}$ out of the top than the
flux $Z^{(0)}=0.3058\,Z_{0}$ that comes into the bottom

$R^{\\{c\\}}=Z^{(0)}-Z^{(1)}=-0.5568Z_{0}.$ (219)

In Fig. 6 we show what happens to the cloud of Fig. 5 if it is allowed to cool
to radiative equilibrium. The incoming radiation $|B^{\\{\rm in\\}}\rangle$
retains the value of (216). In accordance with (203), the Planck intensity of
the cloud is halved from the value of (211).

$B^{\\{c\\}}=\frac{1}{2}B_{0}.$ (220)

Evaluating (169) with $B^{\\{c\\}}$ from (220) and $\dot{M}$ from (209) we
find that the thermally emitted fluxes $|\dot{Z})$ are half as large as those
of (212)

$|\dot{Z})=\left[\begin{array}[]{c}\dot{Z}^{(0)}\\\
\dot{Z}^{(1)}\end{array}\right]=\left[\begin{array}[]{r}-0.2784\\\
0.2784\end{array}\right]Z_{0}.$ (221)

The flux $|\ddot{Z})$ from scattered incoming radiation remains the same as
(217). Summing the thermally emitted flux (221) and scattered flux
$|\ddot{Z})$ of (217) we find that the total flux for radiative equilibrium is

$|Z)=\left[\begin{array}[]{c}Z^{(0)}\\\
Z^{(1)}\end{array}\right]=|\dot{Z})+|\ddot{Z})=\left[\begin{array}[]{r}0.5842\\\
0.5842\end{array}\right]Z_{0}.$ (222)

The cloud heating rate (188) is zero, that is, the cloud thermally emits the
same amount of radiation as it absorbs from incoming radiation

$R^{\\{c\\}}=Z^{(0)}-Z^{(1)}=0.$ (223)

But the cloud is colder than the source of external radiation. The Planck
intensity $B^{\\{c\\}}=B_{0}/2$ of the cloud is only half of the Planck
intensity of $B^{\\{\rm in\\}}_{\bf u}=B_{0}$ of the radiation coming up from
below.

Figure 7:  A stack of $m$ clouds, labeled by the indices $c=1,2,3,\dots,m$. As
described in the text, the figure illustrates the fundamental equations (229)
– (236) for radiation transfer in a cloud stack.

##  3 Cloud Stacks

Having discussed transmission, scattering and absorption of radiation, along
with radiative heating and cooling of single clouds, we turn to stacks of $m$
clouds, like those shown schematically by Fig. 7. Each cloud is identified by
a superscript $\\{c\\}$ where $c=1,2,\ldots,m$. The cloud $c$ will have a
scattering matrix $\mathcal{S}^{\\{c\\}}$. The thermal radiation emitted by
cloud $c$ will be denoted by $|\dot{J}^{\\{c\\}}\\}$. The dotted red arrows of
Fig. 7 indicate the upward and downward parts, $|\dot{J}^{\\{c\\}}_{\bf u}\\}$
and $|\dot{J}^{\\{c\\}}_{\bf d}\\}$. For a nonisothermal cloud,
$|\dot{J}^{\\{c\\}}\\}$ is given by the (110) as

$|\dot{J}^{\\{c\\}}\\}=\int_{0}^{\tau^{\\{c\\}}}d\tau^{\prime}G^{\\{c\\}}(\tau^{\prime})|0)B^{\\{c\\}}(\tau^{\prime}).$
(224)

Here $B^{\\{c\\}}(\tau^{\prime})$ is the Planck intensity at the optical depth
$\tau^{\prime}$ above the bottom of the cloud, and
$|G^{\\{c\\}}(\tau^{\prime})\\}$ is the Green’s function vector for a cloud
with a variable internal temperature, $T(\tau^{\prime})$ at $\tau^{\prime}$.
The vertical optical depth between the bottom and top of the cloud is
$\tau^{\\{c\\}}$. For an isothermal cloud of constant Planck intensity
$B^{\\{c\\}}$ we can use the simpler expression (109),

$|\dot{J}^{\\{c\\}}\\}=\mathcal{E}^{\\{c\\}}|0)B^{\\{c\\}},$ (225)

where the emissivity matrix $\mathcal{E}^{\\{c\\}}$ of the cloud follows from
(110) and is

$\mathcal{E}^{\\{c\\}}=\int_{0}^{\tau^{\\{c\\}}}d\tau^{\prime}G^{\\{c\\}}(\tau^{\prime}).$
(226)

In accordance with Kirchhoff’s law (111), the sum of the emissivity matrix
$\mathcal{E}^{\\{c\\}}$ and the scattering matrix $\mathcal{S}^{\\{c\\}}$ is
the $2n\times 2n$ identity matrix $\hat{1}$.

$\mathcal{E}^{\\{c\\}}+\mathcal{S}^{\\{c\\}}=\hat{1}.$ (227)

We characterize the incoming radiation by the intensity vector
$|\ddot{I}^{(\rm in)}\\}=|\ddot{I}^{(m)}_{\bf d}\\}+|\ddot{I}^{(0)}_{\bf
u}\\}$ of (76), the sum of the downward part $|I^{(m)}_{\bf d}\\}$ of the
intensity vector in the $m$th “gap,” above the top cloud ($c=m$), and the
upward part $|I^{(0)}_{\bf u}\\}$ of the intensity vector in the $0$th “gap”
below the bottom cloud ($c=1$). We recall from (101) that $|\dot{I}^{(\rm
in)}\\}=\breve{0}$ and $|I^{(\rm in)}\\}=|\ddot{I}^{(\rm in)}\\}$. There is no
internally generated thermal radiation coming into the bottom and top of the
stack. In Fig. 7 the incoming intensity is shown as the dotted cyan arrows.

Just as for a single cloud there is an outgoing intensity vector $|I^{(\rm
out)}\\}$ like that of (78), the sum of the downward part $|I^{(\rm out)}_{\bf
d}\\}=|I^{(0)}_{\bf d}\\}$ of the intensity vector in the $0$th gap below the
cloud stack and the upward part $|I^{(\rm out)}_{\bf u}\\}=|I^{(m)}_{\bf
u}\\}$ of the intensity vector in the $m$th gap above the cloud stack. The
downward and upward parts of the outgoing intensity, $|I^{(\rm out)}_{\bf
d}\\}$ and $|I^{(\rm out)}_{\bf u}\\}$ are denoted by the green arrows above
and below the cloud stack of Fig. 7.

The most convenient descriptors of radiation transfer in cloud stacks are the
gap intensities

$\displaystyle|I^{(g)}\\}$ $\displaystyle=$
$\displaystyle|\dot{I}^{(g)}\\}+|\ddot{I}^{(g)}\\}$ (228) $\displaystyle=$
$\displaystyle|I^{(g)}_{\bf d}\\}+|I^{(g)}_{\bf u}\\},$

The gap index can take on the values $g=0,1,2,\ldots,m$. We already defined
the gap intensity $|I^{(0)}\\}$ in the bottom gap, that is, the intensity
below the stack in (80). The gap intensity $|I^{(m)}\\}$, that is, the
intensity above the stack, was defined by (81). In addition, a stack of $m>1$
clouds has $m-1$ internal gaps. For internal gaps $g=1,2,3,\ldots,m-1$ we
denote the intensity in the gap between cloud $c=g$ and the next higher cloud
$c=g+1$ by $|I^{(g)}\\}$ . The downward and upward parts of the internal gap
intensities, $|I^{(g)}_{\bf d}\\}$ and $|I^{(g)}_{\bf u}\\}$ for
$g=1,2,3,\ldots,m-1$ are denoted by the green arrows between adjacent clouds
of Fig. 7.

We see that the radiation in a stack of clouds is determined by $2n(m+1)$
known, nonnegative numbers, the $2n$ stream projections
$(\kern-2.59995pt(\mu_{i}|\ddot{I}^{(\rm in)}\\}$ of the incoming intensity,
and the $2nm$ stream projections
$(\kern-2.59995pt(\mu_{i}|\dot{J}^{\\{c\\}}\\}$ of the thermal source vectors.
These known variables of the problem determine the values of $2nm$ unknown
variables, the $2n$ nonnegative stream projections
$(\kern-2.59995pt(\mu_{i}|I^{(\rm out)}\\}$ of the outgoing radiation, and the
$2n(m-1)$ nonnegative stream projections $(\kern-2.59995pt(\mu_{i}|I^{(g)}\\}$
of the internal gap intensities.

From inspection of Fig. 7 we see that the $2nm$ unknown intensities can be
determined by $2nm$ linear equations, which we abbreviate as $2m$ equations of
$n\times 1$ vectors

$\displaystyle|I^{(1)}_{\bf d}\\}$ $\displaystyle=$
$\displaystyle\mathcal{S}^{\\{2\\}}_{\bf du}|I^{(1)}_{\bf
u}\\}+\mathcal{S}^{\\{2\\}}_{\bf dd}|I^{(2)}_{\bf
d}\\}+|\dot{J}^{\\{2\\}}_{\bf d}\\}.$ (229) $\displaystyle|I^{(1)}_{\bf u}\\}$
$\displaystyle=$ $\displaystyle\mathcal{S}^{\\{1\\}}_{\bf uu}|\ddot{I}^{(\rm
in)}_{\bf u}\\}+\mathcal{S}^{\\{1\\}}_{\bf ud}|I^{(1)}_{\bf
d}\\}+|\dot{J}^{\\{1\\}}_{\bf u}\\}.$ (230) $\displaystyle|I^{(2)}_{\bf d}\\}$
$\displaystyle=$ $\displaystyle\mathcal{S}^{\\{3\\}}_{\bf du}|I^{(2)}_{\bf
u}\\}+\mathcal{S}^{\\{3\\}}_{\bf dd}|I^{(3)}_{\bf
d}\\}+|\dot{J}^{\\{3\\}}_{\bf d}\\}.$ (231) $\displaystyle|I^{(2)}_{\bf u}\\}$
$\displaystyle=$ $\displaystyle\mathcal{S}^{\\{2\\}}_{\bf uu}|I^{(1)}_{\bf
u}\\}+\mathcal{S}^{\\{2\\}}_{\bf ud}|I^{(2)}_{\bf
d}\\}+|\dot{J}^{\\{2\\}}_{\bf u}\\}.$ (232) $\displaystyle\vdots$
$\displaystyle=$ $\displaystyle\vdots$ $\displaystyle|I^{(m-1)}_{\bf d}\\}$
$\displaystyle=$ $\displaystyle\mathcal{S}^{\\{m\\}}_{\bf du}|I^{(m-1)}_{\bf
u}\\}+\mathcal{S}^{\\{m\\}}_{\bf dd}|\ddot{I}^{(\rm in)}_{\bf
d}\\}+|\dot{J}^{\\{m\\}}_{\bf d}\\}.$ (233) $\displaystyle|I^{(m-1)}_{\bf
u}\\}$ $\displaystyle=$ $\displaystyle\mathcal{S}^{\\{m-1\\}}_{\bf
uu}|I^{(m-2)}_{\bf u}\\}+\mathcal{S}^{\\{m-1\\}}_{\bf ud}|I^{(m-1)}_{\bf
d}\\}+|\dot{J}^{\\{m-1\\}}_{\bf u}\\}.$ (234) $\displaystyle|I^{(\rm
out)}_{\bf d}\\}$ $\displaystyle=$ $\displaystyle\mathcal{S}^{\\{1\\}}_{\bf
du}|\ddot{I}^{(\rm in)}_{\bf u}\\}+\mathcal{S}^{\\{1\\}}_{\bf dd}|I^{(1)}_{\bf
d}\\}+|\dot{J}^{\\{1\\}}_{\bf d}\\}.$ (235) $\displaystyle|I^{(\rm out)}_{\bf
u}\\}$ $\displaystyle=$ $\displaystyle\mathcal{S}^{\\{m\\}}_{\bf
uu}|I^{(m-1)}_{\bf u}\\}+\mathcal{S}^{\\{m\\}}_{\bf ud}|\ddot{I}^{(\rm
in)}_{\bf d}\\}+|\dot{J}^{\\{m\\}}_{\bf u}\\}.$ (236)

We see that the downward radiation $|I^{(g)}_{\bf d}\\}$ in the gap $g$ is the
sum of three parts

  1. 1.

$\mathcal{S}^{\\{g+1\\}}_{\bf du}|I^{(g)}_{\bf u}\\}$, the downward reflection
of upward radiation $|I^{(g)}_{\bf u}\\}$ in gap $g$ by the next higher cloud
$c=g+1$.

  2. 2.

$\mathcal{S}^{\\{g+1\\}}_{\bf dd}|I^{(g+1)}_{\bf d}\\}$, the downward
transmission of downward radiation $|I^{(g+1)}_{\bf d}\\}$ in gap $g+1$ by the
cloud $c=g+1$.

  3. 3.

The downward thermal radiation $|\dot{J}^{\\{g+1\\}}_{\bf d}\\}$ emitted by
the cloud $c=g+1$.

The upward radiation $|I^{(g)}_{\bf u}\\}$ in the gap $g$ is the sum of three
analogous parts.

To account for intensity below and above the cloud stack, there are four
exceptional cases in (229) – (236). In (230) and (235) the upward gap
radiation on the right side of the equations is replaced by $|\ddot{I}^{(\rm
in)}_{\bf u}\\}$ of (76) the upward incoming intensity incident on the bottom
of the stack. In (233) and (236) the downward gap radiation on the right side
of the equations is replaced by $|\ddot{I}^{(\rm in)}_{\bf d}\\}$, the
downward incoming intensity incident on the top of the stack. The left sides
of (235) and (236) show how the downward and upward parts of the unknown
intensity $|I^{(\rm out)}\\}$ are related to the other unknown intensities and
to the thermal and external sources.

###  3.1 Stack analysis

To find the unknown intensities on the left sides of equations (229) – (236)
it is helpful to write the equations as the formally simpler stack equation

$|U]=X|U]+P|\dot{J}]+Y|\ddot{I}^{\\{\rm in\\}}\\}.$ (237)

The left side of (237) is the unknown intensity stack vector, which we define
as the concatenation

$|U]=\left[\begin{array}[]{l}|U^{[1]}\\}\\\ |U^{[2]}\\}\\\ \vdots\\\
|U^{[m-1]}\\}\\\
|U^{[m]}\\}\end{array}\right]=\left[\begin{array}[]{l}|U^{[1]}_{\bf d}\\}\\\
|U^{[1]}_{\bf u}\\}\\\ |U^{[2]}_{\bf d}\\}\\\ |U^{[2]}_{\bf u}\\}\\\ \vdots\\\
|U^{[m-1]}_{\bf d}\\}\\\ |U^{[m-1]}_{\bf u}\\}\\\ |U^{[m]}_{\bf d}\\}\\\
|U^{[m]}_{\bf u}\\}\end{array}\right].$ (238)

For $x=1,2,3,\ldots,m-1$ the unknown intensity vectors $|U^{[x]}\\}$ are
identical to the gap intensity vectors $|I^{(x)}\\}$, that is

$|U^{[x]}\\}=|I^{(x)}\\}\quad\hbox{for}\quad x=1,2,3,\ldots,m-1.$ (239)

For the special case of $x=m$, the unknown intensity vector is the same as the
outgoing intensity (78)

$|U^{[m]}\\}=|I^{(\rm out)}\\}=|I^{(0)}_{\bf d}\\}+|I^{(m)}_{\bf u}\\}.$ (240)

We use square brackets $[$ and $]$ as delimiters for the superscript indices
$x$ of the unknown intensity vectors $|U^{[x]}\\}$.

The first term on the right of the stack equation (237) is the product of the
generation matrix $X$ with the unknown intensity stack vector $|U]$ of (238),

$X|U]=\sum_{xx^{\prime}{\bf qq}^{\prime}}X^{[xx^{\prime}]}_{\bf
qq^{\prime}}|U^{[x^{\prime}]}_{\bf q^{\prime}}\\}.$ (241)

The unknown intensity indices, $x$ and $x^{\prime}$, can have the values
$1,2,3,\ldots,m$. The directional indices ${\bf q}$ and ${\bf q}^{\prime}$ can
have the values ${\bf d}$ and ${\bf u}$.

From inspection of (229) – (236) and (237) we see that if $m$ is not so small
that the pattern is distorted by the four exceptions, (230), (233), (235) and
(236), the $n\times n$ block matrix elements in the upper left corner of $X$
are

$X=\left[\begin{array}[]{cccccccccc}\breve{0}&\mathcal{S}^{\\{2\\}}_{\bf
du}&\mathcal{S}^{\\{2\\}}_{\bf
dd}&\breve{0}&\breve{0}&\breve{0}&\breve{0}&\cdots\\\
\mathcal{S}^{\\{1\\}}_{\bf
ud}&\breve{0}&\breve{0}&\breve{0}&\breve{0}&\breve{0}&\breve{0}&\cdots\\\
\breve{0}&\breve{0}&\breve{0}&\mathcal{S}^{\\{3\\}}_{\bf
du}&\mathcal{S}^{\\{3\\}}_{\bf dd}&\breve{0}&\breve{0}&\cdots\\\
\breve{0}&\mathcal{S}^{\\{2\\}}_{\bf uu}&\mathcal{S}^{\\{2\\}}_{\bf
ud}&\breve{0}&\breve{0}&\breve{0}&\breve{0}&\cdots\\\
\breve{0}&\breve{0}&\breve{0}&\breve{0}&\breve{0}&\mathcal{S}^{\\{4\\}}_{\bf
du}&\mathcal{S}^{\\{4\\}}_{\bf dd}&\cdots\\\
\breve{0}&\breve{0}&\breve{0}&\mathcal{S}^{\\{3\\}}_{\bf
uu}&\mathcal{S}^{\\{3\\}}_{\bf ud}&\breve{0}&\breve{0}&\cdots\\\
\vdots&\vdots&\vdots&\vdots&\vdots&\vdots&\vdots&\ddots\\\
\end{array}\right].$ (242)

The block matrix elements in the bottom right corner of $X$ are

$X=\left[\begin{array}[]{cccccccccc}\vdots&\ddots&\vdots&\vdots&\vdots&\vdots&\vdots&\vdots&\vdots\\\
\breve{0}&\cdots&\breve{0}&\breve{0}&\mathcal{S}^{\\{m-1\\}}_{\bf
du}&\mathcal{S}^{\\{m-1\\}}_{\bf dd}&\breve{0}&\breve{0}&\breve{0}\\\
\breve{0}&\cdots&\mathcal{S}^{\\{m-2\\}}_{\bf uu}&\mathcal{S}^{\\{m-2\\}}_{\bf
ud}&\breve{0}&\breve{0}&\breve{0}&\breve{0}&\breve{0}\\\
\breve{0}&\cdots&\breve{0}&\breve{0}&\breve{0}&\breve{0}&\mathcal{S}^{\\{m\\}}_{\bf
du}&\breve{0}&\breve{0}\\\
\breve{0}&\cdots&\breve{0}&\breve{0}&\mathcal{S}^{\\{m-1\\}}_{\bf
uu}&\mathcal{S}^{\\{m-1\\}}_{\bf ud}&\breve{0}&\breve{0}&\breve{0}\\\
\mathcal{S}^{\\{1\\}}_{\bf
dd}&\cdots&\breve{0}&\breve{0}&\breve{0}&\breve{0}&\breve{0}&\breve{0}&\breve{0}\\\
\breve{0}&\cdots&\breve{0}&\breve{0}&\breve{0}&\breve{0}&\mathcal{S}^{\\{m\\}}_{\bf
uu}&\breve{0}&\breve{0}\\\ \end{array}\right].$ (243)

The last two columns of $X$ have only $n\times n$ null matrices, $\breve{0}$,
as elements. For the simplest case of $m=2$ clouds, $X$ is constructed
entirely from the four exceptional equations, (230), (233), (235) and (236),
and can be written as

$X=\left[\begin{array}[]{llll}X^{[11]}_{\bf dd}&X^{[11]}_{\bf
du}&X^{[12]}_{\bf dd}&X^{[12]}_{\bf du}\\\ X^{[11]}_{\bf ud}&X^{[11]}_{\bf
uu}&X^{[12]}_{\bf ud}&X^{[12]}_{\bf uu}\\\ X^{[21]}_{\bf dd}&X^{[21]}_{\bf
du}&X^{[22]}_{\bf dd}&X^{[22]}_{\bf du}\\\ X^{[21]}_{\bf ud}&X^{[21]}_{\bf
uu}&X^{[22]}_{\bf ud}&X^{[22]}_{\bf
uu}\end{array}\right]=\left[\begin{array}[]{llll}\breve{0}&\mathcal{S}^{\\{2\\}}_{\bf
du}&\breve{0}&\breve{0}\\\ \mathcal{S}^{\\{1\\}}_{\bf
ud}&\breve{0}&\breve{0}&\breve{0}\\\ \mathcal{S}^{\\{1\\}}_{\bf
dd}&\breve{0}&\breve{0}&\breve{0}\\\ \breve{0}&\mathcal{S}^{\\{2\\}}_{\bf
uu}&\breve{0}&\breve{0}\\\ \end{array}\right].$ (244)

The second term on the right of (237)

$P|\dot{J}]=\sum_{xc^{\prime}{\bf qq}^{\prime}}P^{[xc^{\prime}\\}}_{\bf
qq^{\prime}}|J^{\\{c^{\prime}\\}}_{\bf q^{\prime}}\\}.$ (245)

is the product of the permutation matrix $P$ with the known thermal source
stack matrix $|\dot{J}]$, which we write as

$|\dot{J}]=\left[\begin{array}[]{l}|\dot{J}^{\\{1\\}}\\}\\\
|\dot{J}^{\\{2\\}}\\}\\\ \vdots\\\ |\dot{J}^{\\{m-1\\}}\\}\\\
|\dot{J}^{\\{m\\}}\\}\end{array}\right]=\left[\begin{array}[]{l}|\dot{J}^{\\{1\\}}_{\bf
d}\\}\\\ |\dot{J}^{\\{1\\}}_{\bf u}\\}\\\ |\dot{J}^{\\{2\\}}_{\bf d}\\}\\\
|\dot{J}^{\\{2\\}}_{\bf u}\\}\\\ \vdots\\\ |\dot{J}^{\\{m-1\\}}_{\bf d}\\}\\\
|\dot{J}^{\\{m-1\\}}_{\bf u}\\}\\\ |\dot{J}^{\\{m\\}}_{\bf d}\\}\\\
|\dot{J}^{\\{m\\}}_{\bf u}\\}\end{array}\right].$ (246)

Here the thermal source vectors $|\dot{J}^{\\{c\\}}\\}$ are given by (224) for
a cloud $c$ with a variable internal temperature, or by (225) for an
isothermal cloud. The downward and upward parts of the source vectors are
$|\dot{J}^{\\{c\\}}_{\bf d}\\}=\mathcal{M}_{\bf d}|\dot{J}^{\\{c\\}}\\}$ and
$|\dot{J}^{\\{c\\}}_{\bf u}\\}=\mathcal{M}_{\bf d}|\dot{J}^{\\{c\\}}\\}$.

From inspection of (229) – (236) and (237) we see that the matrix $P$ is a
block version of a permutation matrix [18]

$P=\left[\begin{array}[]{cccccccccc}\breve{0}&\breve{0}&\hat{1}&\breve{0}&\breve{0}&\cdots&\breve{0}&\breve{0}&\breve{0}&\breve{0}\\\
\breve{0}&\hat{1}&\breve{0}&\breve{0}&\breve{0}&\cdots&\breve{0}&\breve{0}&\breve{0}&\breve{0}\\\
\breve{0}&\breve{0}&\breve{0}&\breve{0}&\hat{1}&\cdots&\breve{0}&\breve{0}&\breve{0}&\breve{0}\\\
\breve{0}&\breve{0}&\breve{0}&\hat{1}&\breve{0}&\cdots&\breve{0}&\breve{0}&\breve{0}&\breve{0}\\\
\vdots&\vdots&\vdots&\vdots&\vdots&\ddots&\vdots&\vdots&\vdots&\vdots\\\
\breve{0}&\breve{0}&\breve{0}&\breve{0}&\breve{0}&\cdots&\breve{0}&\breve{0}&\hat{1}&\breve{0}\\\
\breve{0}&\breve{0}&\breve{0}&\breve{0}&\breve{0}&\cdots&\breve{0}&\hat{1}&\breve{0}&\breve{0}\\\
\hat{1}&\breve{0}&\breve{0}&\breve{0}&\breve{0}&\cdots&\breve{0}&\breve{0}&\breve{0}&\breve{0}\\\
\breve{0}&\breve{0}&\breve{0}&\breve{0}&\breve{0}&\cdots&\breve{0}&\breve{0}&\breve{0}&\hat{1}\end{array}\right].$
(247)

Here $\hat{1}$ denotes an $n\times n$ identity matrix and $\breve{0}$ is an
$n\times n$ null matrix. For the simplest case of $m=2$, the matrix $P$
becomes

$P=\left[\begin{array}[]{llll}P^{[11\\}}_{\bf dd}&P^{[11\\}}_{\bf
du}&P^{[12\\}}_{\bf dd}&P^{[12\\}}_{\bf du}\\\ P^{[11\\}}_{\bf
ud}&P^{[11\\}}_{\bf uu}&P^{[12\\}}_{\bf ud}&P^{[12\\}}_{\bf uu}\\\
P^{[21\\}}_{\bf dd}&P^{[21\\}}_{\bf du}&P^{[22\\}}_{\bf dd}&P^{[22\\}}_{\bf
du}\\\ P^{[21\\}}_{\bf ud}&P^{[21\\}}_{\bf uu}&P^{[22\\}}_{\bf
ud}&P^{[22\\}}_{\bf uu}\\\ \end{array}\right]\\\
=\left[\begin{array}[]{cccc}\breve{0}&\breve{0}&\hat{1}&\breve{0}\\\
\breve{0}&\hat{1}&\breve{0}&\breve{0}\\\
\hat{1}&\breve{0}&\breve{0}&\breve{0}\\\
\breve{0}&\breve{0}&\breve{0}&\hat{1}\\\ \end{array}\right].$ (248)

The last term on the right of (237)

$Y|\ddot{I}^{(\rm in)}\\}=\sum_{x{\bf qq^{\prime}}}Y^{[x]}_{\bf
qq^{\prime}}|\ddot{I}^{(\rm in)}_{\bf q^{\prime}}\\}$ (249)

is the product of the insertion matrix $Y$ with the incoming intensity vector
$|\ddot{I}^{(\rm in)}\\}$ of (76). From inspection of (229) – (236) and (237)
we see that $Y$ has the general form

$Y=\left[\begin{array}[]{ll}Y^{[1]}_{\bf dd}&Y^{[1]}_{\bf du}\\\ Y^{[1]}_{\bf
ud}&Y^{[1]}_{\bf uu}\\\ \vdots&\vdots\\\ Y^{[m-1]}_{\bf dd}&Y^{[m-1]}_{\bf
du}\\\ Y^{[m-1]}_{\bf ud}&Y^{[m-1]}_{\bf uu}\\\ Y^{[m]}_{\bf dd}&Y^{[m]}_{\bf
du}\\\ Y^{[m]}_{\bf ud}&Y^{[m]}_{\bf uu}\\\
\end{array}\right]=\left[\begin{array}[]{ll}\breve{0}&\breve{0}\\\
\breve{0}&\mathcal{S}^{\\{1\\}}_{\bf uu}\\\ \vdots&\vdots\\\
\mathcal{S}^{\\{m\\}}_{\bf dd}&\breve{0}\\\ \breve{0}&\breve{0}\\\
\breve{0}&\mathcal{S}^{\\{1\\}}_{\bf du}\\\ \mathcal{S}^{\\{m\\}}_{\bf
ud}&\breve{0}\\\ \end{array}\right].$ (250)

The four non-zero elements, $Y^{[1]}_{\bf uu}$, $Y^{[m-1]}_{\bf dd}$,
$Y^{[m]}_{\bf du}$, and $Y^{[m]}_{\bf ud}$, are parts of the scattering
matrices $\mathcal{S}^{\\{1\\}}$ and $\mathcal{S}^{\\{m\\}}$ of the bottom and
top clouds of the stack. From inspection of (250) we see that (249)
determines:

  1. 1.

how much the upward incoming intensity $|\ddot{I}^{(\rm in)}_{\bf u}\\}$
contributes to the upward unknown intensity $|\ddot{U}^{[1]}_{\bf u}\\}$ in
the gap above the bottom cloud with $c=1$

  2. 2.

how much the downward incoming intensity $|\ddot{I}^{(\rm in)}_{\bf d}\\}$
contributes to the downward unknown intensity $|\ddot{U}^{[m-1]}\\}$ in the
gap below the top cloud with $c=m-1$

  3. 3.

how much the upward incoming intensity $|\ddot{I}^{(\rm in)}_{\bf u}\\}$
contributes to the downward unknown intensity element $|\ddot{U}^{[m]}_{\bf
d}\\}$. According to (240), this is equal to the downward outgoing intensity
$|\ddot{I}^{\\{\rm out\\}}_{\bf d}\\}$ from the cloud stack,
$|\ddot{U}^{[m]}_{\bf d}\\}=|\ddot{I}^{\\{\rm out\\}}_{\bf d}\\}$

  4. 4.

how much the downward incoming intensity $|\ddot{I}^{(\rm in)}_{\bf d}\\}$
contributes to the upward unknown intensity element $|\ddot{U}^{[m]}_{\bf
u}\\}$. According to (240), this is equal to the upward outgoing intensity
$|\ddot{I}^{\\{\rm out\\}}_{\bf u}\\}$ from the cloud stack,
$|\ddot{U}^{[m]}_{\bf u}\\}=|\ddot{I}^{\\{\rm out\\}}_{\bf u}\\}$

For the simplest case of $m=2$ we see from inspection of (229) – (236) that
the matrix $Y$ simpifies to

$Y=\left[\begin{array}[]{ll}Y^{[1]}_{\bf dd}&Y^{[1]}_{\bf du}\\\ Y^{[1]}_{\bf
ud}&Y^{[1]}_{\bf uu}\\\ Y^{[2]}_{\bf dd}&Y^{[2]}_{\bf du}\\\ Y^{[2]}_{\bf
ud}&Y^{[2]}_{\bf uu}\\\
\end{array}\right]=\left[\begin{array}[]{ll}\mathcal{S}^{\\{2\\}}_{\bf
dd}&\breve{0}\\\ \breve{0}&\mathcal{S}^{\\{1\\}}_{\bf uu}\\\
\breve{0}&\mathcal{S}^{\\{1\\}}_{\bf du}\\\ \mathcal{S}^{\\{2\\}}_{\bf
ud}&\breve{0}\\\ \end{array}\right].$ (251)

###  3.2 Thermal equilibrium

Suppose that all clouds, as well as the external radiation that is incident on
them, have the same temperature and therefore the same Planck intensity $B$.
Then the gap intensities must be given by

$|U]=|0]B,$ (252)

The angular part of the thermal-equilibrium intensity (252) is a stack of $m$
copies of the monopole basis vector $|0)$ of (45)

$|0]=\left[\begin{array}[]{c}|0)\\\ |0)\\\ \vdots\\\ |0)\end{array}\right].$
(253)

For isothermal clouds, each with the same Planck intensity $B$, we can use
(225) to write the thermal source stack vector of (246) as

$|\dot{J}]=[\mathcal{E}]|0]B.$ (254)

The stack matrix for emissivity is block diagonal and given by

$[\mathcal{E}]=\left[\begin{array}[]{llllll}\mathcal{E}^{\\{1\\}}&\breve{0}&\breve{0}&\cdots&\breve{0}&\breve{0}\\\
\breve{0}&\mathcal{E}^{\\{2\\}}&\breve{0}&\cdots&\breve{0}&\breve{0}\\\
\breve{0}&\breve{0}&\mathcal{E}^{\\{3\\}}&\cdots&\breve{0}&\breve{0}\\\
\vdots&\vdots&\vdots&\ddots&\vdots&\vdots\\\
\breve{0}&\breve{0}&\breve{0}&\cdots&\mathcal{E}^{\\{m-1\\}}&\breve{0}\\\
\breve{0}&\breve{0}&\breve{0}&\cdots&\breve{0}&\mathcal{E}^{\\{m\\}}\\\
\end{array}\right].$ (255)

The emissivity matrices $\mathcal{E}^{\\{c\\}}$ of individual clouds were
given by (111).

External radiation in thermal equilibrium is simply Planck radiation, or
blackbody radiation, like that of (142),

$|\ddot{I}^{(\rm in)}\\}=|0)B.$ (256)

Substituting (252), (254) and (256) into (237) we find

$\left(\hat{1}-X-P[\mathcal{E}]\right)|0]=Y|0).$ (257)

The identity (257) can be used as a consistency check of the numerical values
of the matrices $X$, $P$, $[\mathcal{E}]$ and $Y$.

###  3.3 Formal solutions

Rewriting (237) as

$(1-X)|U]=P|\dot{J}]+Y|\ddot{I}^{(\rm in)}\\},$ (258)

we see that the formal solution is

$\displaystyle|U]$ $\displaystyle=$
$\displaystyle(\hat{1}-X)^{-1}P|\dot{J}]+(\hat{1}-X)^{-1}Y|\ddot{I}^{(\rm
in)}\\}$ (259) $\displaystyle=$ $\displaystyle G|\dot{J}]+W|\ddot{I}^{(\rm
in)}\\}$ $\displaystyle=$ $\displaystyle|\dot{U}]+|\ddot{U}].$

According to the last line of (259), the unknown intensity stack vector $|U]$
is a linear combination of the known thermal source stack vector $|\dot{J}]$
of (246) and the known incoming intensity vector $|\ddot{I}^{(\rm in)}\\}$ of
(76). We will discuss the proportionality matrices $G$ and $W$ below.

The contribution to (259) of thermal emission by cloud particulates and gas
molecules is

$|\dot{U}]=G|\dot{J}]\quad\hbox{or}\quad|\dot{U}^{[x]}\\}=\sum_{c=1}^{m}G^{[xc\\}}|\dot{J}^{\\{c\\}}\\}.$
(260)

Here the stack matrix $G$ of (260) can be written as the $m\times m$ block
matrix

$G=(\hat{1}-X)^{-1}P=\left[\begin{array}[]{ccccccc}G^{[11\\}}&G^{[12\\}}&\cdots&G^{[1,m-1\\}}&G^{[1,m\\}}\\\
G^{[21\\}}&G^{[22\\}}&\cdots&G^{[2,m-1\\}}&G^{[2,m\\}}\\\
\vdots&\vdots&\ddots&\vdots&\vdots\\\
G^{[m-1,1\\}}&G^{[m-1,2\\}}&\cdots&G^{[m-1,m-1\\}}&G^{[m-1,m\\}}\\\
G^{[m,1\\}}&G^{[m,2\\}}&\cdots&G^{[m,m-1\\}}&G^{[m,m\\}}\\\
\end{array}\right].$ (261)

The $2n\times 2n$ discrete Green’s matrices $G^{[xc\\}}$, that are the
elements of the block matrix (261), are themselves $2\times 2$ block matrices,

$G^{[xc\\}}=\left[\begin{array}[]{ccccccc}G^{[xc\\}}_{\bf dd}&G^{[xc\\}}_{\bf
du}\\\ G^{[xc\\}}_{\bf ud}&G^{[xc\\}}_{\bf
uu}\end{array}\right]\quad\hbox{where}\quad\mathcal{M}_{\bf
q}G^{[xc\\}}\mathcal{M}_{\bf q^{\prime}}.$ (262)

Both the row index $x$, which labels unknown intensities, and the column index
$c$ of $G^{[xc\\}}$, which labels clouds, can have the values
$1,2,3,\ldots,m$.

We can use the last block element of the matrix equation (260) to write the
equivalent thermal source vector $|\dot{J}^{\\{\rm ev\\}}\\}$ of the entire
stack as a linear combination of the thermal source vectors of the individual
clouds,

$\displaystyle|\dot{J}^{\\{\rm ev\\}}\\}$ $\displaystyle=$
$\displaystyle|\dot{I}^{(\rm out)}\\}=|\dot{U}^{[m]}\\}$ (263)
$\displaystyle=$ $\displaystyle\sum_{c=1}^{m}G^{[mc\\}}|\dot{J}^{\\{c\\}}\\}.$

For the special case that each cloud is isothermal with a Planck intensity
$B^{\\{c\\}}$ and a corresponding thermal source vector
$|\dot{J}^{\\{c\\}}\\}=\mathcal{E}^{\\{c\\}}|0)B^{\\{c\\}}$ from (225), we can
simplify (263) to

$\displaystyle|\dot{J}^{\\{\rm ev\\}}\\}$ $\displaystyle=$
$\displaystyle|\dot{I}^{(\rm out)}\\}=|\dot{U}^{[m]}\\}$ (264)
$\displaystyle=$
$\displaystyle\sum_{c=1}^{m}G^{[mc\\}}\mathcal{E}^{\\{c\\}}|0)B^{\\{c\\}}.$

If all clouds have the same temperature and the same Planck intensities
$B^{\\{c\\}}=B$, (264) can be further simplified to

$\displaystyle|\dot{J}^{\\{\rm ev\\}}\\}$ $\displaystyle=$
$\displaystyle|\dot{I}^{(\rm out)}\\}=|\dot{U}^{[m]}\\}$ (265)
$\displaystyle=$ $\displaystyle\mathcal{E}^{\\{\rm ev\\}}|0)B.$

The emissivity of the isothermal cloud stack is a weighted sum of the
individual cloud emissivities,

$\mathcal{E}^{\\{\rm ev\\}}=\sum_{c=1}^{m}G^{[mc\\}}\mathcal{E}^{\\{c\\}}.$
(266)

From inspection (266) and (139) we see that in the limit of clouds $c$ of
infinitesimal optical thickness, $\tau^{\\{c\\}}\to d\tau^{\prime}\to 0$, the
discrete Green’s matrix $G^{[mc\\}}$ and the emissivity
$\mathcal{E}^{\\{c\\}}$ approach the limits

$\frac{G^{[mc\\}}}{\tau^{\\{c\\}}}\to
G(\tau^{\prime})\quad\hbox{and}\quad\mathcal{E}^{\\{c\\}}\tau^{\\{c\\}}\to\hat{1}\,d\tau^{\prime}\quad\hbox{as}\quad\tau^{\\{c\\}}\to
d\tau^{\prime}\to 0.$ (267)

Here $G^{[mc\\}}$ is the discrete Green’s matrix of (276), $G(\tau^{\prime})$
is the continuous Green’s matrix of (139), and $\tau^{\prime}$ is the total
optical depth from the bottom of the cloud stack to the bottom of the cloud
$c$.

The contribution to (259) of incoming radiation is

$|\ddot{U}]=W|\ddot{I}^{(\rm
in)}\\}\quad\hbox{or}\quad|\ddot{U}^{[x]}\\}=W^{[x]}|\ddot{I}^{(\rm in)}\\}.$
(268)

The scattering coefficients $W^{[x]}$ are generalized scattering matrices that
give the contribution to the unknown gap intensity $|U^{[x]}\\}$ from the
incoming intensity $|I^{\\{\rm in\\}}\\}$ after transmissions and reflections
by all of the clouds of the stack. The $m\times 1$ block matrix $W$ of
scattering coefficients is

$W=(\hat{1}-X)^{-1}Y=\left[\begin{array}[]{cc}W^{[1]}\\\ W^{[2]}\\\ \vdots\\\
W^{[m-1]}\\\ W^{[m]}\\\ \end{array}\right].$ (269)

The elements $W^{[x]}$ of the block matrix (269) are themselves $2\times 2$
block matrices, with $n\times n$ elements $W^{[x]}_{\bf qq^{\prime}}$,

$W^{[x]}=\left[\begin{array}[]{ccccccc}W^{[x]}_{\bf dd}&W^{[x]}_{\bf du}\\\
W^{[x]}_{\bf ud}&W^{[x]}_{\bf
uu}\end{array}\right]\quad\hbox{where}\quad\mathcal{M}_{\bf
q}W^{[x]}\mathcal{M}_{\bf q^{\prime}}.$ (270)

We write the last element of the matrix equation (268) as

$|\ddot{U}^{[m]}\\}=W^{[m]}|\ddot{I}^{(\rm in)}\\}.$ (271)

or since $|\ddot{U}^{[m]}\\}=|\ddot{I}^{(\rm out)}\\}$ according to (240),

$|\ddot{I}^{(\rm out)}\\}=\mathcal{S}^{\\{\rm ev\\}}|\ddot{I}^{(\rm in)}\\}.$
(272)

From inspection of (271) and (272) we see that the last scattering coefficient
$W^{[m]}$ is the equivalent scattering matrix of the entire cloud stack,

$\mathcal{S}^{\\{\rm ev\\}}=W^{[m]}.$ (273)

As we will discuss in more detail in Section 4, the equivalent scattering
matrix $\mathcal{S}^{\\{\rm ev\\}}$ of a cloud stack is the Redheffer star
product [19] of the scattering matrices of the individual clouds,

$\mathcal{S}^{\\{\rm
ev\\}}=W^{[m]}=\mathcal{S}^{\\{m\\}}\raisebox{-1.0pt}{*}\mathcal{S}^{\\{m-1\\}}\raisebox{-1.0pt}{*}\cdots\raisebox{-1.0pt}{*}\mathcal{S}^{\\{2\\}}\raisebox{-1.0pt}{*}\mathcal{S}^{\\{1\\}}.$
(274)

An explicit expression for the star product of two clouds is given in Section
4.3 as Eq. (382). The general expression (274) for a stack of $m\geq 2$ clouds
follows by mathematical induction.

In summary, we can use (263) and (271) to write the total output intensity,
the sum of a part due to thermal emission of cloud particulates and gas
molecules, and a part due to transmission, scattering and absorption of
incoming radiation, as

$|I^{(\rm out)}\\}=|\dot{J}^{\\{\rm ev\\}}\\}+\mathcal{S}^{\\{\rm
ev\\}}|\ddot{I}^{(\rm in)}\\}.$ (275)

The thermal source vector $|\dot{J}^{\\{\rm ev\\}}\\}$ of the equivalent cloud
was given by (263) for non-isothermal clouds and by the simpler expression
(264) for isothermal clouds. For a stack of isothermal clouds, we can
substitute the emissivity matrix (266), and the scattering matrix (273) into
Kirchhoff’s law (111), to find the identity,

$\sum_{c=1}^{m}G^{[mc\\}}\mathcal{E}^{\\{c\\}}+\mathcal{S}^{\\{m\\}}\raisebox{-1.0pt}{*}\mathcal{S}^{\\{m-1\\}}\raisebox{-1.0pt}{*}\cdots\raisebox{-1.0pt}{*}\mathcal{S}^{\\{2\\}}\raisebox{-1.0pt}{*}\mathcal{S}^{\\{1\\}}=\hat{1}.$
(276)

Eq. (276), a discretized version of the Kirchhoff identity (139), can provide
a useful check of numerical calculations.

Using the stack Green’s matrix $G$ to find $|\dot{U}]=G|\dot{J}]$ in (260),
and using the scattering coefficient matrix $W$ to find
$|\ddot{U}]=W|\ddot{I}^{(\rm in)}\\}$ in (268), have the advantage of physical
clarity. But this is not the most efficient and accurate way to numerically
solve (259). For extensive numerical calculations, it is better to avoid
calculating the inverse matrix $(1-X)^{-1}$ needed to evaluate $G$ and $W$,
and to evaluate the contribution $|\dot{U}]$ of internal thermal emission and
$|\ddot{U}]$ of incoming radiation with backslash fractions,

$|\dot{U}]=(\hat{1}-X)\backslash P|\dot{J}],$ (277)

and

$|\ddot{U}]=(\hat{1}-X)\backslash Y|\ddot{I}^{\\{\rm in\\}}\\}.$ (278)

The algorithms used by modern mathematical software like Matlab [20] to
evaluate backslash fractions like (277) and (278) are more efficient and much
more accurate than those used to evaluate the inverse matrix
$(\hat{1}-X)^{-1}$.

###  3.4 Gap intensities and fluxes

Having solved (259) for the stack vector $|U]$ of unknown intensities, we can
use the results to evaluate the gap intensity stack vector

$|I)=\left[\begin{array}[]{l}|I^{(0)}\\}\\\ |I^{(1)}\\}\\\ |I^{(2)}\\}\\\
\vdots\\\ |I^{(m-1)}\\}\\\ |I^{(m)}\\}\\\ \end{array}\right].$ (279)

The gap intensity stack vector $|I)$ of (279) and the unknown intensity vector
$|U]$ of (238) are almost the same, but the small differences are important.
The gap intensity vector $|I)$ of (279) has $m+1$ block elements,
$|I^{(g)}\\}$ for $g=0,1,2,\ldots,m$, in contrast to the $m$ block elements
$|U^{[x]}\\}$ for $x=1,2,3,\ldots,m$ of the unknown intensity vector $|U]$ of
(238). According to (239) the $m-1$ middle block elements
$|I^{(g)}\\}=|U^{[g]}\\}$ for the internal gap indices $g=1,2,3,\ldots,m-1$.
The first block element $|I^{(0)}\\}$ of the gap intensity vector $|I)$ of
(279) was given by (80), which we rewrite as

$\displaystyle|I^{(0)}\\}$ $\displaystyle=$ $\displaystyle|I^{(\rm out)}_{\bf
d}\\}+|I^{(\rm in)}_{\bf u}\\}$ (280) $\displaystyle=$
$\displaystyle|U^{[m]}_{\bf d}\\}+|\ddot{I}^{(\rm in)}_{\bf u}\\}.$

The second line of (280) comes from (101) and (240). The last block element
$|I^{(m)}\\}$ of the gap intensity vector $|I)$ of (279) was given by (81),
which we rewrite as

$\displaystyle|I^{(m)}\\}$ $\displaystyle=$ $\displaystyle|I^{(\rm in)}_{\bf
d}\\}+|I^{(\rm out)}_{\bf u}\\}$ (281) $\displaystyle=$
$\displaystyle|\ddot{I}^{(\rm in)}_{\bf d}\\}+|U^{[m]}_{\bf u}\\}.$

The second line of (281) comes from (101) and (240).

The gap intensity stack vector $|I)$ of (279) is fully characterized by
$2n(m+1)$ nonnegative numbers $(\kern-2.59995pt(\mu_{i}|I^{(g)}\\}$. One can
also fully characterize the intensity with the $2n(m+1)$ multipole moments

$I^{(g)}_{l}=(\kern-2.59995pt(l|I^{(g)}\\},$ (282)

where the possible values of the multipole index are $l=0,1,2,...,2n-1$ and
the possible values of the gap index are $g=0,1,2,\ldots,m$. Expressions for
the multipole amplitudes $I^{(g)}_{l}$ were given in Section 2.2. Of
particular importance are the dipole moments $I^{(g)}_{1}$. According to (52)
we can write the scalar gap fluxes $Z^{(g)}$ as

$Z^{(g)}=4\pi I^{(g)}_{1}=4\pi(\kern-2.59995pt(1|I^{(g)}\\}.$ (283)

We will denote the $(m+1)\times 1$ array of scalar gap fluxes by

$|Z)=\left[\begin{array}[]{l}Z^{(0)}\\\ Z^{(1)}\\\ Z^{(2)}\\\ \vdots\\\
Z^{(m-1)}\\\ Z^{(m)}\end{array}\right].$ (284)

###  3.5 Model cloud stacks

To gain insight into radiation transfer for cloud stacks we consider the
simplified model where each cloud is isothermal and has a Planck intensity
$B^{\\{c\\}}$, and where the incoming radiation $|\ddot{I}^{(\rm in)}\\}$ is
the half-isotropic radiation of (172). We assume that the scattering matrices,
$\mathcal{S}^{\\{c\\}}$, and therefore the emissivity matrices,
$\mathcal{E}^{\\{c\\}}=\hat{1}-\mathcal{S}^{\\{c\\}}$, are known for all of
the clouds. Then the radiation transfer is determined by only $m+2$ initial
conditions: the $m$ Planck intensities $B^{\\{c\\}}$ of the isothermal clouds
and the upward and downward Planck intensities, $B^{\\{\rm in\\}}_{\bf u}$ and
$B^{\\{\rm in\\}}_{\bf d}$ of the half-isotropic incoming radiation.

We first consider the scalar fluxes of internal gaps $g=1,2,3,\ldots,m-1$,
when $|I^{(g)}\\}=|U^{[g]}\\}$ in accordance with (239). Using (260) with
$|\dot{J}^{\\{c\\}}\\}=\mathcal{E}^{\\{c\\}}|0)B^{\\{c\\}}$ in accordance with
(225), we write the thermally generated part of (283) as

$\dot{Z}^{(g)}=4\pi\sum_{c=1}^{m}(\kern-2.59995pt(1|G^{[gc\\}}\mathcal{E}^{\\{c\\}}|0)B^{\\{c\\}}.$
(285)

We can use (268) and (172) to write the contribution of incoming radiation to
the internal gap flux as

$\ddot{Z}^{(g)}=4\pi\sum_{\bf q=d}^{\bf
u}(\kern-2.59995pt(1|W^{[g]}\mathcal{M}_{\bf q}|0)B^{\\{\rm in\\}}_{\bf q}.$
(286)

Summing (285) and (286), we find that the total scalar flux for the internal
gap $g$ is

$\displaystyle Z^{(g)}$ $\displaystyle=$
$\displaystyle\dot{Z}^{(g)}+\ddot{Z}^{(g)}$ (287) $\displaystyle=$
$\displaystyle
4\pi\sum_{c=1}^{m}(\kern-2.59995pt(1|G^{[gc\\}}\mathcal{E}^{\\{c\\}}|0)B^{\\{c\\}}+4\pi\sum_{\bf
q=d}^{\bf u}(\kern-2.59995pt(1|W^{[g]}\mathcal{M}_{\bf q}|0)B^{\\{\rm
in\\}}_{\bf q}$ $\displaystyle=$ $\displaystyle
4\pi\sum_{c=1}^{m}\dot{M}^{(gc\\}}B^{\\{c\\}}+4\pi\sum_{\bf q=d}^{\bf
u}\ddot{N}^{(g)}_{\bf q}B^{\\{\rm in\\}}_{\bf q}.$

For the internal gaps with $g=1,2,3,\ldots,m-1$, the coupling matrices are

$\displaystyle\dot{M}^{(gc\\}}=(\kern-2.59995pt(1|G^{[gc\\}}\mathcal{E}^{\\{c\\}}|0)\quad\hbox{for}\quad
c=1,2,3,\ldots,m,\quad\hbox{and}\quad\ddot{N}^{(g)}_{\bf
q}=(\kern-2.59995pt(1|W^{[g]}\mathcal{M}_{\bf q}|0).$ (288)

The gap fluxes $Z^{(0)}$ below the stack, and $Z^{(m)}$ above, require special
attention. Using (264), (268) and (273) we write (280) as

$\displaystyle|I^{(0)}\\}$ $\displaystyle=$ $\displaystyle|U^{[m]}_{\bf
d}\\}+|\ddot{I}^{(\rm in)}_{\bf u}\\}$ (289) $\displaystyle=$
$\displaystyle\mathcal{M}_{\bf
d}\left(\sum_{c=1}^{m}G^{[mc\\}}\mathcal{E}^{\\{c\\}}|0)B^{\\{c\\}}+W^{[m]}|\ddot{I}^{(\rm
in)}\\}\right)+\mathcal{M}_{\bf u}|\ddot{I}^{(\rm in)}\\}$ $\displaystyle=$
$\displaystyle\sum_{c=1}^{m}\mathcal{M}_{\bf
d}G^{[mc\\}}\mathcal{E}^{\\{c\\}}|0)B^{\\{c\\}}+\left(\mathcal{M}_{\bf
u}+\mathcal{M}_{\bf d}W^{[m]}\right)|\ddot{I}^{(\rm in)}\\}.$

Multiplying the last line of (289) on the left by $4\pi(\kern-2.59995pt(1|$,
and using (172) for $|\ddot{I}^{(\rm in)}\\}$ we write the scalar flux (283)
for the gap $g=0$ as

$\displaystyle Z^{(0)}$ $\displaystyle=$ $\displaystyle
4\pi\sum_{c=1}^{m}(\kern-2.59995pt(1|\mathcal{M}_{\bf
d}G^{[mc\\}}\mathcal{E}^{\\{c\\}}|0)B^{\\{c\\}}+4\pi\sum_{\bf
q}(\kern-2.59995pt(1|(\mathcal{M}_{\bf u}+\mathcal{M}_{\bf
d}W^{[m]})\mathcal{M}_{\bf q}|0)B^{\\{\rm in\\}}_{\bf q}$ (290)
$\displaystyle=$ $\displaystyle
4\pi\sum_{c=1}^{m}\dot{M}^{(0c\\}}B^{\\{c\\}}+4\pi\sum_{\bf q=d}^{\bf
u}\ddot{N}^{(0)}_{\bf q}B^{\\{\rm in\\}}_{\bf q}.$

The coupling matrices are

$\displaystyle\dot{M}^{(0c\\}}=(\kern-2.59995pt(1|\mathcal{M}_{\bf
d}G^{[mc\\}}\mathcal{E}^{\\{c\\}}|0)\quad\hbox{and}\quad\ddot{N}^{(0)}_{\bf
q}=(\kern-2.59995pt(1|(\mathcal{M}_{\bf u}+\mathcal{M}_{\bf
d}W^{[m]})\mathcal{M}_{\bf q}|0).$ (291)

Using (264), (268) and (273) we write (281) as

$\displaystyle|I^{(m)}\\}$ $\displaystyle=$ $\displaystyle|U^{[m]}_{\bf
u}\\}+|\ddot{I}^{(\rm in)}_{\bf d}\\}$ (292) $\displaystyle=$
$\displaystyle\mathcal{M}_{\bf
u}\left(\sum_{c=1}^{m}G^{[mc\\}}\mathcal{E}^{\\{c\\}}|0)B^{\\{c\\}}+W^{[m]}|\ddot{I}^{(\rm
in)}\\}\right)+\mathcal{M}_{\bf d}|\ddot{I}^{(\rm in)}\\}$ $\displaystyle=$
$\displaystyle\sum_{c=1}^{m}\mathcal{M}_{\bf
u}G^{[mc\\}}\mathcal{E}^{\\{c\\}}|0)B^{\\{c\\}}+\left(\mathcal{M}_{\bf