On the fundamental solution of the radiative transfer equation

JOURNAL OF GEOPHYSICAL RESEARCH, VOL. 93, NO. D9, PAGES 11,037-11,050, SEPTEMBER 20, 988 On the Fundamental Solution of the Radiative Transfer Equation PIOTR J. FLATAU AND GRAEME L. STEPHENS Department of Atmospheric Science, Colorado State University Fort Collins, Colorado This paper outlines the general solution of the one-dimensional, azimuthally averaged reliative transfer equation in terms of a matrix exponential. The link between this fundxnental solution and those more commonly used in radiative transfer is established. The formulation is developed for a general vertically inhomogeneous atmosphere with sources. Several new concepts, based on properties of the matrix exponentials, are described in the context of radiative transfer, including the use of the commutator an.d product integrals. It is also demonstrated how the matrix exponen- tial formulation provides for new insights, not only into improvements of the numerical efficiency and stability of the solution, but also into the understanding of raliative transfer through a lay- ered atmosphere. The various concepts introduced in this paper are illustrated throughout by the two-stream simplification of the general radiative transfer equation. 1. INTRODUCTION The transfer of solar and infrared radiation through the terrestrial atmosphere is one of a class of transport phe- nomena that are important to many studies in the engineer- ing and physical sciences (e.g., see Duderstradt and Mar- tin [1979] for a general review). Radiative transport theory, which began in its present form with the work of Schuster in 1905, has been widely used in many diverse problems, such as in the calculation of the radiative energy budgets of stars, planetary atmosph_.eres, crops, cattle, and sheep; in atmo- spheric and underwater visibility studies; in the determina- tion of the reflecting properties of snow, ice sheets, paint sur- faces, photo emulsions, and of blood tissues and viruses; in the estimation of the insulating properties of fiberglass ma- terials; and in the research on the thermal structures of blast furnaces and other high-temperature machinery, among sev- eral others. The interested reader is referred to the Journal of Quantitative Spectroscopy and Radiative Transfer, 36 (1), 1986, for examples of many different applications of radia- tive transfer. Because of the diversity of topics that use radiative trans- fer theory, a wide range of solution methodologies of the radiative transfer equation (RTE) can be found in the lit- erature (as exemplified in the monograph of Lenoble [1985] and in the review article of Stamnes [1986] ). While there seems to be an abundance of different methods of solution, most start by replacing the continuous dependence of radi- ation intensity on direction by a discrete set of directions. The discretized transfer equations then take the form of a coupled system of first-order differential equations. An anal- ogous set of coupled equations also results when the spatial Fourier transform method is employed on the full three- dimensional transfer equation [e.g., Stephens, 1986, 1988]. As it turns out, such systems can be represented in a classi- cal mathematical way in terms of a matrix-exponential op- eration, in which the array of coupling coefficients appears in the exponent. This procedure has been applied to various problems in a limited way; notably in studies involving elec- tromagnetic wave propagation [Keller, 1960], in seismology [Frer and Frazer, 1984], and in radiative transfer [Water- Copyright 1988 by the American Geophysical Union. Paper number 8D0281. 0148-0227/88/008D--0281505.00 man, 1981]. Even though many methods exist for solving the radiative transfer equation, there is still a need to develop more efficient numerical techniques. This need is prticu- larly apparent in those studies of radiative transfer set on two- and three-dimensional geometries [e.g., Stephens, 1988]. The object of this paper is to formulate the solution of the radiative transfer equation in terms of the matrix exponen- tial and to demonstrate the link between this fundamental solution and those more commonly encountered in radia- tive transfer literature. This formulation is developed for a general vertically inhomogeneous atmosphere with sources. Several concepts, such as the commutator and commuta- tor identities, product integral, extinction transformation, and two-point boundary value problems, are introduced in this paper and are relatively new to the context of radiative transfer but are actively used in other physical disciplines. -- The formulation described in this paper presents the framework for the unification of the many different exist- ing methods of solution and not only provides new insights into the possible improvements in numerical efficiency and stability, but also provides important new insight into ra- diative transfer through a vertically inhomogeneous atmo- sphere. Examples can be found in the discussion of small ' solutions, the discussion of stability, and the discussion of the commutator. These aspects, to our knowledge, are new to radiative transfer and warrant even further consideration than is given here. The formulations described in this paper use the two-stream representation of the RTE throughout to illustrate the different concepts introduced. As such, the examples provide new and meaningful insight on the math- ematical structure of the two-stream model. The outline of the paper is as follows. Section 2 briefly in- troduces the radiative transfer equation in its matrix form. The fundamental solution of this matrix equation is then represented in terms of the propagator function and the boundary conditions. Section 3 describes methods for evalu- ating the matrix exponential exp(A') particularly for small '. The interaction principle and the extinction operator are introduced in section 4, and the connection between the ex- tinction operator and the Riccati equations is demonstrated. Section 5 describes the eigenvector method for evaluating exp(._-) in terms of similarity transformations and demon- strates that the use of a Riccati transform is, in fact, a sim- ilarity transformation. The solution for a vertically inho- mogeneous atmosphere is presented in section 6, in which 11,037

11,038 ILATAU AND STEPlIENS: SOLUTION OF THE RADIATIVE TRANSFER EQUATION adding, doubling, initialization, stability, and the commuta- tor are discussed. Examples of the latter two concepts are presented in the context of a two-stream model. 2. BASIC CONCEPTS It will be supposed throughout this paper that the absorbing-scattering medium is horizontany homogeneous. The radiative transfer equation describing the transfer of monochromatic radiation through such a medium is then d (, , ) =  if, , ) _    &  (..;'.') ß I (%t',qb') dp'dqb' q- S (%p, qb) (1) where I (% p, b) is the specific intensity at level ' along the direction represented by the angle pair p, b, an where - is the optical depth, b the azimuth angle, and p the cosine of the polar angle, and the wavelength dependence on all variables is taken to be understood. We consider, for convenience, that the intensity is az- imuthally independent, in which case (1) simplifies to dI o f_+x x   (' ) = z (' ) - T  if' ' ') () ß  (, ') ' + s (, ) where I, p, and S are now taken to be azimuthally averaged quantities. Discretization of RTE Consider the medium as shown in Figure 1 and represent the continuous intensity field as a set of 2n discrete directions (as shown). If we choose these directions according to some quadrature scheme, then the integral (i.e., path radiance term of equation (2)) becomes  (,,') (, ') ' ,(,,,)(,,) 1 i=l +  , (.. .-,)  (..-,) With the definition of the n-dimensional specific intensity vector I as (a similar definition applies for source vector S) (,-) z(.,+) z(.,-) z(., +) I, = . I. = . (4) (,-.) (,+.) and with the definition of the fonowing matrices o o p = Tp(, ) i,y, = ,...,  (5) the discretized, matrix form of (2) becomes dIa - rn- = Ia - plwlt - pbwl. + St (6) dl. rn-- = I. - pwl - plwl. + $. (7) n streams ._ v K .- Fig. 1. Schematic representation of the radiation flow arid dis- cretition of the medium. The convention used for the optical depth scale to the left is for the layered medium, while the  convention to the right applies to the continuous case. where we use the convection p and p! for the symmetric phase matrices associated with back and forward scattering, respectively. Subscripts d and tt refer to the sense of the flow in the downward and upward directions, respectively. Equations (6) and (7) can be more succinctly written as dl : AI + $ ($) where now the 2n x 2n matrix A is composed of the n x  matrices ax and a2 which are defined by al = m-lp.w- m -I (10) a2 = m-Xpaw (11) and I = I. ' The matrices ax and a2 depend on local optical properties of the atmosphere (such as the single scattering albedo o and the asymmetry parameter !/). Since we use two-stream models as a means of illustrating various of the concepts in- troduced throughout this paper, we now present the appro- priate form of these coeffidents for the Eddington approx- imation. From Meador and Weaver [1980], for example, it readily follows 1 ax = - [7 -- &0 (4 + 3!/)] (13)

FLATAU AND STEPHENS' SOLUTION OF THE uADTI q.a. NSFER EQUATION 11,039 - a0 (4 - where al and a2 are now scalar quantities (note also that the two-stream version of equation (6) follows with n = 1). Slightly different definitions of these coefficients which corre- spond to different types of two-stream models, as discussed by Meador and Weaver [1980], could also have been used. The Eddington approximation is used only as a matter of convenience. The generality of the results that follows will not be altered with different forms of ax and a2. Fundamental Solution of the Matrix RTE It can be shown (e.g., refer to mathematical texts on or- dinary differential equations (ODE), such as that of Cole [1968], and Gantmacher [1959]) that for 2(,r) continuous on ,rt _ ,r _ ,rb, there exists a unique solution to the vector differential equation dl = l(.r)l (15) d,r with (one-point) boundary conditions l(-rt). This solution may be written in the form I (r) = Pal (rt) (16) where Pa is a 2n x 2n fundamental solution matrix which is the unique matrix satisfying the matrix differential equation = A (,r) ra (17) with Pa (.rt) = 1. (Here the subscript d reminds us that \"propagation\" occurs from the top of the medium downward to level The propagator for a homogeneous atmosphere. Consider the solution of the homogeneous part of (8) for constant A; then (16) may be written as I (.r) = eA('-\")l (.r) (18) for a layer of thickness .r-.ft. As discussed later, the matrix exponential P(% .r,)=e A('-\") (19) contains most of the information about the solution of the RTE. It is called the fundamental solution or the matriser in the theory of ordinary differential equations [Cole, 1968; Coddington and Levinson, 1955; Gantmacher, 1959], the transfer matrix solution [Aronson and Yarmush, 1966] in neutron physics, and the propagator [Kennett and Kerry, 1979; Gilbert and Backus, 1966; Fryer and Frazer, 1984] in seisinology. In oceanography, Preisendorfer [1988] uses the term fundamental matrix, and in optimal control, Denman [1986] calls P the transition matrix. We have chosen to call this factor the propagator, in keeping with its interpreta- tive meaning. For example, given the boundary conditions at one point, say I (.r) in (16), then the solution deeper in the medium can be recovered (propagated) down from this boundary by applying this propagator. However, it has to be stressed in this context that the radiative transfer equation can not be treated as an initial value problem. The propaga- tor matrix as such is not a solution to the RTE. Only after the boundary conditions are applied are the solutions re- covered. The above-mentioned \"boundary conditions at one point\" mean both sets of upward and downward intensities, and only one set is known in advance. This by no means reduces the usefulness of propagator, as demonstrated more fully later. Inhomogeneous atmosphere: Definition of the product in- tegral. It can be shown that the propagator Pa(-r,-rt) be- tween the top of the atmosphere .rt and the internal point .r can be written as a product of individual fundamental solutions of the sublayers [Gilbert and Backus, 1966; Gant- macher, 1959; Karp et ed., 1980] in the form Pa (,r, ,r) = Pa (,r, ,r,-) Pa (,r,-x, ,r,-.) x'\" x P(,r, ,r) = H Pa (.r, .r_) (20) Notice that the continuous interval (.r, -rt) has been divided here into k layers and that -rt = -r0 and .r = -rk (see Fig- ure 1). Equation (20) has simple intuitive appeal: starting with the intensities at the top of the medium (which are as- sumed to be known) we get intensities at the internal point -rx by applying the downward propagator Pa (-rx, -rt). Op- erator Pa (-r2, -rx) is then applied to these intensities to get the intensities even further down into the layer, this time at point -r2. Continuing this procedure, solutions at each inter- nal point are obtained, and Pa (.r, -rt) is constructed. For a vertically contiguous medium composed of k homogeneous layers, (20) can be expressed, with the aid of (19), as x... x e A(\"'\")xn = H eA('\"\"-')x\" (21) In the limiting case of infinitesimally small sublayers for which assumption about homogeneity is satisfied, we can define the product integral T ra (r, rt) = rleA(') a, (22) Because we use the subscript d in definitions (20, (21), and (22), it is logical to assume that the operator Pa does not depend on the order of its arguments. Propagation of the solution from the bottom to the top of the atmosphere will be indicated by subscript . The appendix contains more details on properties of product integrals, which will prove to be useful in further analysis of (21). Boundary Conditions In atmospheric radiative transfer problems one usually en- counters a two-point (TP) boundary value problem (BVP) where for solar radiation, say, the downward intensities at cloud top and the upward intensities at cloud base reflected from an underlying surface, are specified. The emerging intensities from cloud top and base are then the required quantities. Thus the set of 2n linear ordinary differential equations givenby (8) with A defined by (9) is then sup- plemented by 2n two-point boundary conditions at the top, .r = -rt = 0, and bottom, .r = ,rb, of the atmosphere. These boundary conditions can be expressed as b = B,I (.rt) + BI (.r) (23) where Bt and B, are 2n x 2n matrices and b and I are 2n coltann vectors. Usually, boundary conditions are separable,

11,040 FLATAU AND STEPHENS: SOLUTION OF THE I:[ADIATIVE TRANSFER EQUATION which leads to Be and Bb of the form Bt-- 0 0 and (24) 0 0] (2) B,-- 0 1 Here 1 is an n x n unit matrix and 0 is an n x n matrix whose entries are all zeros. Notice that vector b contains information on boundary conditions at both points. Using (16), we can write (23) as DI (r,) = b (26) where the characteristic matrix for the boundary problem D is defined as D = B,Pa (pt, pt) + B,Pt (r,, pt) (27) and where (26) gives upward and downward intensities I (r) at one point, in this case the top of the atmosphere, in terms of boundary conditions b specified at two points. Note also that Pa(rt, r) = 1. More explicitly, inverting (26) yields I (r,) = D-b (28) wtfich provides the full solution to the external response problem (without sources S). Combining (16) and (28) we also get intensities at every internal point of the atmospheric layer between 't _< ' _< ', (again for a sourceless problem). The analogy between the present formulation of the internal and external response problems in terms of the propagator and the response problems discussed in terms of doubling and adding methods will be made more evident later. Example of boundareal conditions for the two-stream model. Consider the sourceless two-stream model equations, which dr I,, -a. -a I,, The boundary conditions for I(r,) nd I(n) can be written in the notation employed so far as 0 .() ] (a0) The characteristic matrix for this boundary value problem is given by (27), with B, and B, defined by the 2 x 2 analog of (24) and (25), respectively. Using the propagator defined later by (100), we get D=n,+--n x+ +---n x- and obtain the solution to the internal problem as () () Sozt'ce$ Knowing the fundamental matrix of the homogeneous part of RTE (8), one can construct a particular solution for the given source function S [e.g., Cole, 1968; Bellman, 1970; Gantrnacher, 1959] as r () = r () r () S () d (2) Thus the intensities in an inhomogeneous atmosphere with arbitrary sources are obtained as I (r) = Pa (r) D-b + P (r) P (s) S (s) ds (33) or more simply as I(r) = ra(r)I(r,) + F(r) (34) Using the Green's function notation, we define the following matrix function as [Russell, 1985] G(t, s) = r(r)D-B(r,)r-(s) G(r, s) = -Pa(r)[1 - D -B(r,)]Pa(s) -x () and write (33) as I(r) = Pa(') D-b+ G(% s)S(s) ds (36) Two examples of source functions in the two-stream model. Radiative transfer problems often encounter two different types of source functions; one that represents the source of diffuse intensity as a result of the single scatter of a collimated beam of radiation and a second that repre- sents the source of intensity that arises from thermal emis- sion within the atmosphere. The latter is encountered in infrared transfer studies and is often approximated by S, = ( -- o) (to + r) So (37) which assumes the emission to be a linear function of '. The former type of source function is usually encountered in problems dealing with solar radiative transfer and is ex- pressed as S,ola,. = aae-\"/øSo (38) where the constant vector So is defined as For a homogeneous atmosphere there is, at least concep- tually, no difficulty in evaluating (32). The propagator is exponential in - and source terms such as (37) or (38) are relatively simple functions of r. Thus in some instances the integral of (32) can be evaluated analytically. Let us again provide the following example, based on two-stream theory, for an illustration. Consider a homogeneous atmo- sphere with constant A, for which we have from (32) and ()  - ,x('-')s()a (40) Substitution of (100) into (40) yields  X)Soh( ' , ,) +  ( _ x  = ( + ¾  ¾)s0(-, , ,) (4) where the function h(k, , ') in (41)is h(, , ,) = '('-') i() a, (42) and where f(') is given by (37) or (38) without the So factor. In both cases, So is defined by (39). The integral in (42) can be evaluated for a homoge- neons atmosphere, and the two-stream approximation to

iLATAU AND STEPHENS: SOLUTION OF THE IIATIVE TRANSFER EQUATION 11,041 the RTE's particular solution is given analytically for two important cases of source functions: ha(k, r, O) ---- rBx (Bx + kBo) (e k, _ 1) (43) k } k 2 ( hoo(k, r, O) = 1 + kilo e  - e- For the more complex case of the n-stream problem and for a more general form of source function, F can be obtained numerically by techniques to be described later. 3. EVALUATION OF THE MATRIX EXPONENTIAL We now describe some methods used to evaluate the prop- agator. The review article by Moler and van Loan [1978] provides a much more extensive discussion of the subject of calculating exp(i') and is recommended to the interested reader. We will first consider those methods that provide an approximation to exp(ir) for small r, as these will be utilized in later discussions. We also describe the similarity transformation, which provides the basis for the eigenvector approach often used in radiative transfer literature. Taylor Series Decomposition of Matrix Exponent The formal Taylor series definition of (19) is A A e =X+Ar+..-+ n! +'\" (45) where I is a 2n x 2n unit matrix. This series exists for all A and ' and is convergent [Cole, 1968; Gantrnacher, 1959; Campbell, 1980]. Suppose now that ' is small and consider omy ,o,d-od tm  , the,, umg (4) d (, w obtain e A*= l+Ar+- +'\" 0 1 al. a2 T 2 al. a2  --a2 --ax -- --a2 --ax (46) It is shown later that (46) constitutes the existing methods of initialization for the doubling and adding schemes. The assumption behind (46) is that Ar is somehow small. A more quantitative representation of this criterion can be IIi II < 1 w denote the matrix norm of A by II A II (e Co,, [1960] and Moler and van Loan [1970] for definitions of various matrix norms and their properties). Rational Approximation and Ward 's Algorithm Here we describe a method based on the work of Ward [1977], which uses the diagonal Pad table rational approxi- mation to compute the exponential matrix. Following Ward [1977], the pth diagonal Pad4 approximation to the expo- nential of the matrix i is given by Q- (-A)Q,(A) (47) where Q,(X) is a polynomial in X of degree p defined by Q,(X) ---- E c'X' (48) I--0 where e x [02I ' '- eX ' ' ' ' ' [ ...... 3rd order Pad' I ...... 3rd order Taylor exponsion 'øøto.% x Fig. 2. Compaxisons of various approximations to the exponent for small m. (2p- k)!p! (49) c -- (2p)!k!(p- The error in approximating the matrix exponent by diagonal Pad entries increases as the norm of the matrix increases. Therefore one of the objectives of Ward's [1977] algorithm is to decrease the matrix norm to some desirable prescribed level. The algorithm uses three techniques to achieve such a decrease in matrix norm. The first technique uses trans- formations based on the estimated value of the eigenvalues. The second technique attempts to minimize the norm of A over all possible diagonal similarity transformations. The third technique exploits the identity exp(i) -- (exp(i/m) \" (property 2, equation (153)). The general idea of this tech- nique is to choose ra as a power of 2 starting with exp which has been determined reliably and efficiently, say, using (46) and then to form the matrix [exp (A/m)]'\" by repeated squaring. Thus it is required to choose m as the smallest of 2 which II II/m _< 1, This mhod of squ- ing (doubling) has been used in mathematical literature to evaluate other functions which involve matrix exponentials (as, for example, in the work of van Loan [1978]). This same technique is also equivalent to the doubling scheme which is a method frequently used in radiative transfer literature. We now provide an example of the diagonal Pad approxi- mation to e'. With some manipulation, it is straightforward to show that ,. = + + + (50) 1 -- x/2 + x2/10-- xa/120 Figure 2 presents a comparison of e ', third-order power se- ries of exp(x) and Paa. It can be seen that the rational function gives considerably better approximation to exp(x) than the Taylor series, and it is to be anticipated that (50) will provide a better and more efficient approximation to r (i) th th,t gi,, by (4). Similarity Transformations We now pursue the matrix decomposition technique which is closely related to the discrete ordinate method, but the description offered here is more general. All of the matrix decompositions are based on similarity transformations of the form [e.g., Moler and van Loan, 1978]

11,042 FLATAU AND STEPHENS: SOLUTION OF THE RADIATIVE TRANSFER EQUATION A = TCT - (51) and after some algebra involving the Taylor series definition (45), can be shown that (51) implies A, TeC, T- (52) The idea now is to find the similarity transformation T for which eC is easy to compute. In particular, if C is diagonal of the form C ---- diag (kx, .. ', k2,) (53) where ki are eigenvalues associated with A, then the expo- nent of C is trivial to calculate, namely, C. (e k\" , ?\"') (54) e -- diag ,... In the degenerate case of repetitive eigenvalues the Jordan canonical form could be used [see Karp et al., 1980] , and there is no difficulty in finding exponent of this Jordan form [Moler and van Loan, 1978] . Splitting The familar additive law for scalar exponentials xp(,O = + doesn't hold when one deals with matrices, as the following short calculation shows. Using the Taylor series expansion, we get -A e-A, e A,+As _ (I- A2 + A e  2 .... ) + -y .... + + + 2 Multiplying terms on the right-hand side, we obtain e-Ae-A'e A+A = I-- [A,A2] +'\" where we defined the commutator as [A, A:] = AA: - A:A (58) Only if Ax and A2 commute, i.e., [A1,A:] = 0, does the additive law hold. Using (157), we can finally write the following formula, which is valid through quadratic terms in Ax and A2 e A'+A' = eA'eA'e-[A\"A']/ .-- (59) The exponential identity (59) can serve as a method to find matrix exponent when the matrix A can be split into matrices Ax, Am, for which their exponents and the expo- nents of their commutators can be calculated efficiently. It is also sometimes possible to express (59) by a finite number of factors. This can be achieved if the decomposition of A into basis matrices Ai forms solvable Lie algebra [ Wei and Norman, 1964; Magnus, 1954; Sauter, 1979]. In such a case it is possible to use Baker-Hausdorf-Carnpbell exponential identities of the type (59) to develop a dosed-form solution. (ro by 4. THE INTERACTION PRINCIPLE AND THE EXTINCTION OPERATOR The fundamental solution, characterized by the propa- gator, does not depend on boundary conditions in any way and is determined solely by the internal properties of a layer. These are the single scattering properties and their distribu- tion within the layer. The reflectance and transmittance ma- trices which are defined in the following discussion similarly depend only on internal characteristics of a layer and, once known, provide a solution for any arbitrary incident radia- tion field. Thus knowledge of reflectance and transmittance contains all information about the internal characteristics of a layer, as does the propagator. It is therefore reasonable to expect that these two concepts are closely related. It is the purpose of this section to illustrate the connection between these inherent properties of the layer. The significance of the interaction principle in establishing this link becomes evident through the definition of the extinction operator. Interaction Principle, Extinction Matrix, and Extinction Transformations Let us define, in the usual way, the reflectance and trans- mittance matrices as describing the amount of radiation re- flected from the layer and transmitted through the layer. These matrices are composed of bidirectional reflectances and transmittances which are, in turn, characterized by in- cident and response directions of the flow. In general, radi- ation coming from above will be reflected and transmitted differently than radiation incident from below. As above we use the subscripts d and u to indicate downward and upward propagating radiation. Given the incident intensities Ia(',) on the top of the layer, we get the following responses Ia (r) = taI (r,) (60) I. (r) =rtlt (r) (61) and for the case of upward intensities incident on the bottom of the layer, we get a similar set of responses; namely, I, = I. (r) = t.I. (r) (63) In general, radiation comes from both directions and inter- acts linearly with the medium. Thus (60) - (63) can be combined and, with the additional contributions by inter- nal sources, we obtain the matrix form of the interaction prindple I.(%) = R I.(r,) + G (64) where we define the extinction matrix R as and the source-related vector ( as Our goal now is to find a relationship between the 2n X 2n extinction matrix and the 2n x 2n propagator matrix. To this end, we define, after Saastamoinen [1980], the following matrix operations. Consider an arbitrary 2n x 2n matrix M partitioned into n x n blocks according to m2! m22

iLATAU AND STEPHENS: SOLUTION OF THE RADIATIVE TRANSFER EQUATION 11,043 We define the downward operator Sa as $a(M)= [mxx-mx2mm2 m2m] --! --rn m m d the upwd operator S by [ m' --m'mx ] $(M) = m2m' m:2 - m:m'm: (68) (69) It is simple to check that $. [$a(M)] = M, $ [$(M)] = M (70) and that the inverse of the extinction operator is the same as the operator itself, i.e., S(M) = si  (M) () and &(M) = SC(M) (7) Thus with matrix multiplication we find that Sa(M)S,,(M) = S,,(M)Sa(M) = X (73) IL = S [Pa(, 'a)] (82) and IL = S,, [P,,('a, 're)] (83) Let us remark about the inhomogeneous term. Again, com- paring (75) to (34), one has F = Pax G (84) which provides the transformation rule for inhomogeneous terms in the propagator and the extinction formulations of the RTE. Riccati Equations It is possible to derive a set of differential equations for the extinction matrix. Differentiating (80) with respect to using (17) and applying the $ transformation to P, we get aS(t) = S(t) (85) from which four Riccati-type matrix equations for t., and ta are determined as Relation Between Extinction and Fundamental Matrices We now utilize the matrix definitions given earlier. Con- sider the interaction principle (64) for intensities emerging from the top and at the bottom of the layer. With some rearrangement and with regrouping of the intensities, (64) can be written as r.I.(%! Ia(',, [2:](74, [0 ,.] [:: _0,] which is equivalent to the following form of the fundamental solution where we identify that [ [ ] 1 --r 1 --rt-! Pax= 0 --t = 0 --t -x (76) and Comparing (75) with (16), we see that PaO't, 'b) is Pa(, 't): Pax (, ', )Pa (, ',) (78) and with simple multiplication we Pa( ,) = [ ta- r'txra r'tx] (79) ' _tXra t Th (79) gives the relatiop between the fdental solution d the extction matr (65) [ Pa(,) = Sa(R) = ta--rtXra r.t x _t;Xra t; x (80) ' rat x t. -- ratXr Since the extinction operator is its own inverse, we get from () a (SO) dr dT : a2 -alr -rual -ra2r (86) drd = --a -b axr -- raax -- rar (87) d dta d7' = --axt -- arat (88) dt -- ait. + ar.t. (89) d' The extinction transformation of the two-stream propaga- tor. Consider a homogeneous layer and assume that r,, = ra and t = ta. It will be shown later that the propagator for the two-stream model is given by (100) which, in expanded form, is  = (,. ,. ' (o) 2k -- e- ) e-- --% + + = ( +=/), = = ( - /), a m where a, transformation S, applied to (90), produces 1 R = $a(e A'r) = (e,r/2)( 1 _ az/k) q- (e-'r/2)(1 -I- ax/k) [ 1 ß (\"2/2]½) (e 'r -- e-'') I (91) Thus with the definition of P, according to (65), we readily identify the transmission as 1 ---- e,r/2(l_ax/k) Jr_e_,r/2(1Jr_az/k) (92) and the reflectance as a212k(e\" --e -'v)  = ($) e'r/2(1 -- az/k) -I- e-\"/2(1 --[- az/k) which correspond to the expressions given by Meador and Weaver [1980] for their model 2. The etinction transforwaation of the n-stream model. For the n-stream model we use the eigenvalue decomposition,

1,044 17LATAU AND STEPHENS: SOLUTION OF THE lqADIATIVE 'INSFER EQUATION which will be discussed in more detail later. Applying the S transformation to (105), we get ( ekz+) -x t ---- u_e-k'z_ + u+ (94) where u+, zñ are cconls to Notice that (94) and (95) can be manipulated such that only negative exponents are present tAror, on and Yavmush, 1966; Aronson, 1970; Saastamoinen, 1980; Fryer and Frazer, 1984]. We will return to this aspect of (94) and (95) in our discussion on stability. 5. EIGENVECTOR ANALYSIS: APPLICATIONS TO THE RTE The Method of Discrete Ordinates One approach to evaluate the matrix exponential exp (A') is to find the eigenvalues ki and associated elgen- vectors ui of the matrix A.. This is the essence of the discrete ordinate method of Chandrasekhar [1960] , as discussed by jima and Tanaka [1986] in recent literature. The method discussed in these papers is based on the factorization or decomposition of the matrix A. (The terminology \"discrete ordinates\" is in fact misleading, as the discreteness of the method which results by replacing the continuous angular dependence by a set of discrete directions is not unique to this method. We propose that this method be more appro- priately referred to as the eigenvector method.) The traditional derivation of the discrete ordinate method avoids finding an explicit expression for the propagator by imposing boundary conditions and solving for matching co- efflclents arising from the requirements of continuity of in- tensity across each of the boundaries. The boundary condi- tions are therefore an integral part of the solution scheme, and a new set of calculations are necessary for derent set of boundary conditions. By contrast, the propagator is de- fined by the internal structure of the atmosphere and is not dependent in any way on boundary conditions. Therefore the fundamental solution is obtained first, and the bound- ary conditions are applied next. The Propagator of a Two-Stream Model We first consider the two-stream model as a relatively simple example of the eigenvector method. From the Caley- Hamilton theorem [e.g., Graham, 1979] we know that some function of a 2n x 2n matrix A can be represented by a polynomial r (A) in A of the order 2n- 1. For the case of a two-stream model, n --- 1, and we can represent the propagator as A e = o1 -[- A (96) To get a and/3, we first find the eigenvalues of Ax = kx (97) which gives k 2 = a- a, and since we require both e A and r (A) to possess the same eigenvalues [Pennisi, 1972; Graham, 1979], it follows that (9s) and -' f] (99) After simple algebra we finally arrive at e = e-k (1 -- 2 for the two-stream propagator where 1 is the 2 x 2 identity matrix. The Eigenvalues and Eigenvectors of œn Flow Model Eigenvalues of the 2n x 2n matrix A defined by (9) have been studied extensively by Aronson [1970], Aronson and Yammush [1966] and, more recently, by Waterman [1981], Stamnes and Swanson [1981], Nakajima and Tanaka [1986], and Preisendorfer [1988]. Aronson [1971], studied rela- tionship between the transfer matrix and Case's [Case and Zweifel, 1967] method of eigenfunctions. The structure of the eigenvectors is important when considering stability of numerical schemes. They occur in 4-k pairs and are all real [Stamnes and Swanson, 1981; Aronson, 1970]. The eigen- values are also invariant with respect to similarity transfor- mations. Let us define the 2n x 2n transformation matrix and note that T -x  T. Thus ----TAT-X---- d 0 where s = ax - a2 and d -- ax -F a. and (101) (102) = 0 ds Because matrix sd is ds transposed, it is enough to find eigenvectors of one of them. After Aronson [1970], let the diagonalization of sd be given by k  = y-x (sd)y (104) Thus it can be shown that e =  u_ u+ 0 e k z_ z+ 00) where u+ = (y 4-x) (106) ..+ = + and x -- dyk -x -- s-Xyk (108) For the 2n flow problem under consideration, u+, u_, z+, z_, s, d, a, a2, x, y, and k are n x n dimensional matrices. Invariant Imbedding The more traditional discussion of invariant imbedding usually focuses on the Riccati equations, such as derived earlier. The earlier works of Bellman and Kalaba [1959], Redheifer [1962], and Preisendorfer [1958] and the review of

:LATAU AND STEPHENS: SOLUTION OF THE I:[ADIATIV rIIqANSFER EQUATION 11,045 Scott [1975] provide ample illustration of the method. The method is a powerful tool for analyzing linear systems [e.g., Denman, 1979; 1986] and is frequently employed in current research on two-point boundary value problem for ODE's [19s4] mong While the method of invariant imbedding has been de- scribed in many different ways, we demonstrate here that the approach is a form of similarity transformation of matrix .&, much in the same way as the eigenvector analysis described earner. This connection is perhaps not surprising given that invariant imbedding techniques have been employed on cer- tain eigenvalue problems (e.g., Scott, [1975] and Preisendor- far, [1973] among others). Under this transformation, the real virtue of invariant imbedding and the reason perhaps why it is popular in two-point boundary value problems is made apparent as we demonstrate that this transformation yields the matrix of block diagonal form and thus transforms the two-point boundary value problem to a more tractable and separable pair of initial value problems. We start with a transformation T, which might be viewed as more general than that described by (103), as it is contin- uous in '. We now apply this transformation to the prop- agator defined by (22), with the end result given by the similarity rule of (161), namely, IX [X()d4 T(n) = H exp T - ()A(s)T(s) - T - (') dT(s) ds ds (09) and desire that the resulting transformation matrix be block diagonal to provide a stable solution to (8). Let us define the transformation matrix T(?') as 1 [1 --r.] (110) T(')-- r, -1 and impose the condition that r and r satisfy the matrix Riccati equations (85) and (86), then it can be shown that T - T dT() _ [ al- t2r& o ] --ax q- a2r,, (111) and, by virtue of the block diagonal form of (111), the prob- lem has been transformed into two uncoupled initial value problems, which are stable with respect to the direction of the solution. Dichotomy, Stability, and Decoupling We have already remarked that our general matrix A has a spectrum such that for every eigenvalue k/ there is an eigenvalue -ki. The solution, in the form of the propagator exp(A') is thus a combination of both growing and de- caying exponentials as indicated by (105). For large optical thicknesses, growing exponentials will dominate the solution given in this form, and it follows that the naive applica- tion of the propagator will lead to numerical instability. In such circumstances the true solution is lost and cannot be recovered. The behavior just described is referred to as di- chotomic and is the subject of vigorous ongoing research in mathematical literature [Russell, 1985; vanLoon, 1985; de- Hook and Mattheij, 1985; Mattheij, 1985; Smith, 1987] and [Evans, 1986]. In the case of the two-stream model, this dichotomy means that there is one exponentially increas- ing and one exponentially decreasing fundamental solution component [Russell, 1985]. In the case of a 2n flow model, n growing and n decaying modes can be expected. The issue of numerical stability has been realized in the radiative transfer literature for some time. To illustrate the nature of this in- stability and the remedies introduced to overcome them, we consider the two-layer atmosphere. We confine our interests to the case of a sourceless atmosphere for brevity and desire to determine the fluxes at each of the boundaries including at the internal layer. We employ the two-stream formula- tions in this illustration, although the same arguments apply to the more general 2n stream model. The method com- monly employed is exemplified in the two-stream work of ShetHe and Wainman [1970] and Pollack and McKee [1985] and in the 2n stream work of Liou [1975] and Starehas [1986] among others. In these studies the solutions for each indi- vidual layer are determined and the boundary conditions at the interface boundaries are matched, resulting in a system of simultaneous equations. In the present formulation the solutions for a two-layer atmosphere can be expressed as (112) where obviously, matrix multiplication serves as a matching condition between boundary of two layers. In Figure 3 we show each of the components of the propagator (denoted as Pij) as a function of increasing optical thickness. One of the components (namely, P22) becomes unbounded, while the others stay bounded. Thus the solution either in the form of the propagator or in the form of matching boundary con- ditions, as mentioned earlier will be numerically unstable. This can be overcome by employing the extinction formu- lation in the form of reflection and transmission functions, since they are well behaved and bounded. 6. THE CASE OF A VERTICALLY INHOMOGENEOUS ATMO SPHERE We have so far described different ways of evaluating exp (A.') using small ' approximations and similarity trans- 5.0 4.0 5.0 2.0 .0 0 i o:0.8 g: 0.9 ./ / o / ß / - ß / ß / ß _ / /' \"P2,2 ' P,I 0 1.0 2'.0 3.0 4.0 Optical Thickness r Fig. 3. Components of two-stream propagator as a function of increasing 'r.

11,046 LATAU AND STEPHENS: SOLUTION OF THE PtADIATIV, TRANSFElt F, QUATION formations. We will now show, by a study of the properties of products of exponentials of a matrix, that the propagator can be readily used to construct solutions of RTE for an in- homogeneous atmosphere and arrive at the adding scheme as a result. The virtue of the adding scheme was actually made apparent in the earlier discussion of stability and is of- ten used in radiative transfer studies. The relationship of the adding scheme to the propagator has already been realized in the atmospheric sciences literature by Waterman [1981], Karp et al. [1980], and Stamnes [1986]. Peebles and Ples- set [1951] applied a similar scheme in gamma ray transfer. Aronson and Yarmush [1966] studied the transfer matrix so- lution in connection with neutron transfer. Scattering and transfer matrix methods are used in optimal control [e.g., Denman, 1986], while Kennett [1974] showed how this tech- nique can be applied to seismic waves. The Adding Method for a Sourceless Medium We now treat a vertically contiguous atmosphere com- posed of a succession of overlying homogeneous layers and obtain composite expressions for reflectance and transmit- tance. We assume, for the moment, that there is no internal source, so that F = 0 and G = 0. Using (20) and (83) the extinction matrix for the layer ('t, ) is R(rt, rb) = Sa [Pa(r, rb)] = Sa[ fi Pa(ri+x, ri)] (113) Since each component in the product of (113) can be ex- pressed by the extinction matrix (65) through the extinction transformation, (113) can be written in the form Equation (114) actually represents the well-known adding scheme [e.g., Grant and Hunt, 1968 a, b, 1969a, b]. It allows us to construct reflectances and transmittances of combined layers from our knowledge of reflectances and transmittances of individual layers, and it is equivalent to the star product composition rule of Redheifer [1962]. Consider, for illustra- tion, an atmosphere composed of two layers extending from ('b, ') and (% -). The extinction matrix for the two com- bined layers follows from (114) and (113) as n) = S IS S n)]] (115) With the notation R(, ') = R*, and I(', ') = R where the asterisk is used to indicate the lower layer properties and, after some algebra involving matrix multiplication and S transformations, we can write 1 t (, ) ---- tt; (116) 1 * -- rdr  r,(-t, 'b) = r + tr -- rdr  t, (117) and 1 r,(-, 'b) = ra + t rta (118) 1 * -- rdr  1 * t, (119) t(, n) = tl_ rr which represents the adding formula for a two-layer atmo- sphere. The Adding Scheme for an Inhomogeneous Atmosphere With Sources Knowing the propagator, one can construct a particular solution for an arbitrary source function and write down the full solution as (33). For an inhomogeneous atmosphere the total propagator is given by product integral (20). In the same manner, it is natural to ask if the knowledge of F for individual, homogeneous layers is sufficient to reconstruct F for the entire inhomogeneous atmosphere. The answer is yes, and an equation similar to (20) can be derived as follows. Consider as above two layers (t, r,) and then from (34) we write I(') = P(% ',)I(',) + F(% 't) (120) and I() = P(, r)I(r) + r(r,, r) (121) Combining (120) and (121), we see that F(, ') = P(', ')F0', ) + F(',, ') (122) This is the adding formula for source terms in propagator formulation. Thus (84) and (122) may be used to obtain the adding of source terms in the extinction formulation. Doubling method. Doubling can be viewed as a partic- ular application of adding. In doubling we consider a ho- mogeneous layer for which A is constant and the upward and downward components of the extinction matrix are the same. For the sake of discussion, we choose to use the prop- erty of matrix exponentials, such as defined by (152), to derive the doubling algorithm. Consider a layer of optical thickness ' for which the extinction matrix R, defined by (65), is given. From (152) and (114) we get = (x2a) Applying the S transformation to R and using (70), we ob- tain s. yt(2)] = s. (124) Finally, from (68), it follows that [ t-(2 - ) --t- (2-)r(2-) ] r(2t)t -x(2t) t(2t) - r(2t)t -(2t)r(2t) _ [ t-(t) -t-(t)r(t) ]2 -- r(-)t-x( -) t(-) -- r()t-x()r() from which t(2t) = t(t)(1 - r(t)2)-xt(t) and r(2t) = r(t) + t(t)r(-)(1 - r(t)z)-xt(t) (125) (126) can be derived. This derivation defines the doubling scheme. It is usually employed to derive the extinction matrix for homogeneous layers. Clearly, an intimate relationship exists between Ward's [1977] algorithm for the calculation of a matrix exponent and the doubling formula. Initialization methods. So far in the discussion of the ex- tinction operator, it was assumed that the extinction matrix (127)

:LATAU AND STEPHENS: SOLUTION OF THE IADmTIVE TRANSFER EQUATION 11,047 R(') is known. We now discuss methods by which is evaluated or 'qnitiallzed\", so that the composition rules described earlier can be employed. The method of initial- ization provides such an estimate of this matrix under the constraint that the medium is homogeneous. To this end, we could use, for example, the Taylor approximation (45), say, to the second order in ', to approximate A A 2 [1 0] e I+A-+  +.-- 0 1 al a2 T2 1 + + -a -a  -a -ax (128) Again using the S transformation (68), we obtain the trans- reitrance function as T2 / --1 T2 d the reflectance fction  T2 r = a2 + y(axa2 + a2ax) (130) Equations (129) d (130) were derived by Wiscombe [1976] d are referred to  eded ond tia- eressio up to td-order tern  . Gnt and Hunt [1968a,b, 1969a,b] e a sar seme, but oy to a st or- der   d referred to ts method  the tes gen- erator tization (IGI). As o ysis shows, EDI d IGI e both Taylor series decompositions of eonti. As we ssed previoudy, Was [1977] ysis cates that P approts may play some role  tiaHzation semes. It is cle that the contio for sm optical t  d o no of , ny, II  /  1, d the error anysis of tiation semes shoed - dude such estate. e tiization semes considered up to now apply to a homogenous atmosphere without soces. Coider  extdon to ts w ows us to itize for soces. To ts end, we need to evuate tegr such  oc for the pti solution F giv by (32), i.e., F = eA('-')Sds (131) where both the 2n x 2n matr A d the 2n x I soce vector S are assreed to be cott. Vector S  a 2n-eiona] eqlt of (37) or (38) eluated for a pti sublayer. It c be se om (84) that om F we c recover the socrelated fction . ter van Loan [1978], coider the (2n + 1) x (2n + 1) blo trig matr 0 0 and it can be proved that (132) eHr [P F I (133) = 0 f where P and F follow our standard notation, with f a scalar function of '. By equating submatrices in the identity (150) specified for matrix H and with boundary condition exp-Ir] I=o= 1, we get ___d p = Ap and dF=AF+S and f = 1. The solution to these equations gives (133). Forming H and finding its matrix exponent with, for ex- ample, Ward's [1977] algorithm allows us to find both the fundamental solution and the source function integral F. Splitting Methods and the Commutator We begin our discussion by considering the fundamental solution for two combined layers extending from ('o, 'x) and ('x, '2). After (14), this solution requires the determination of the propagator P('-, ') -- e A('-n) e A'(\"-'ø) (134) where Ai and A are matrices appropriately determined from the optical properties of each layer. It is tempting now to generalize (134) (and (21)) and replace them with (/.') exp A.(') d' (135) which is a possibility motivated by intuition from the Baker- Hausdorf-Campbell(BHC) exponential-type identity (59). Examining (59), we notice that additional terms to (135) appear, thus = + ... (136) where these unspecified additional terms involve integrals of omm. to [(), ()] .a igh omm. to e. Although the commutator may seem only distantly rele- vant to radiative transfer, it convenlently describes sev- eral properties which we will now consider. We see that exp(f A(t)dr) describes only part of the interaction of ra- diation with the medium. The presence of a nonvanishing commutator is a statement about asymmetry in the distri- bution of optical properties. In other words, two clouds (say, one composed of large droplets in the upper layer and small droplets in a bottom layer, and the second cloud with the reversed structure)may have the same f. A(t)dr, but they still will interact differently with radiation. This difference is characterized by the the noncommutative property of matri- ces A, which tells us that two layers interact in a nonlinear way. Only in the small thickness limit is the interaction linear and do the matrices commute. We will now explore the two-stream version of (136) to illustrate this feature. Before doing so, it is convenient to consider several variants of (136) . Here we follow Magnus [1954] and Wilcox [1967] and assme that P(', ', -- 0) -- e v(') (137) where V(') is a series V(') = Z V,,(') (138) Using exponential identities [e.g. Stainberg, 1988] one can derive the following formulas [ Wilcox, 1967] for V

! ! ,048 ILATAU AND STEPHENS: SOLUTION OF THE t)DIATIVE TRANSFER EQUATION = - (140) = i v,.], + v,.], (141) ß ß where Ax -- A('x), etcß and the higher-order terms can be calculated from recursire procedures. The representation provided by (137) is of local nature as it applies to some local level '. Wei and Norman [1964] discuss a more global approach based on product of expo- nentials P(r, r, -- 0) = H eV\"(') (142) The advantage of the global representation (142) is that only finite amount of terms is required (see also Sauter [1979] and Mulholland [1971] ). Let us now consider an example of the commutator for the two-stream model. It follows from the definition of A for the two-stream model that (9) can be written in the form A=C+O (14a) where and --1/4 --7/4] (144) C -- 7/4 1/4  --(4 -- 3g) --(4 q- 3g) (145) Consider two overlying homogeneous cloudy layers, one of optical thickness 'x and characterized by a single scatter albedo &z and the second layer characterized by properties '2 and &2, and assume that the asymmetry parameter g is the same for both cloud layers. Such a situation might occur in reality, when, for instance, the upper portions of a cloud layer contain characteristically larger droplets (and hence smaller values of &o) than the lower layers of the cloud. A relevant question to ask of such a case is whether or not the radiation reflected from the cloud is altered significantly when the positions of these layers are reversed. Some assess- ment of this problem can be made using (17) and (13S). Suppose, for convenience, that rz ---- r --r. Then it follows that and Vx = (Xx q- X.)- (140) where, according to (143), the commutator becomes X,] -- (O, - It, O] [ 0 a(g- 1)] = 0 (147) (148) Applying now the S transformation (69), one can recover extinction matrix and, in turn, the reflection and transmis- sion properties. 7. SUMMARY AND CONCLUSIONS This paper outlines the general solution to the one- dimensional, azimuthally averaged radiative transfer equa- tion in terms of a matrix exponential. This solution is devek oped for a vertically nonhomogeneous atmosphere The link between the fundamental form of solution developed here to those more commonly used in radiative transfer (such as that of doubling/adding and Chandrasekhar's discrete ordinate approach) is established, and it is suggested that the matrix exponential formulations present the framework for the unification of the many different existing methods of solution. This paper demonstrates how the methods of invariant imbedding and discrete ordinates are in fact spe- cial forms of simi]arity transformation. It is apparent from discussions presented in this paper that the doubling and adding schemes provide a stable algorithm for the solution of the two-point boundary value problems typically encoun- tered in radiative transfer. It was also further demonstrated how the approach might be made more numerically efficient by using Pad approximant. The paper furthermore pro- vides examples of the two-stream solutions formulated as a matrix exponential and provides some new insight into the mathematical structure of the two-stream model. Some con- cepts, new to radiative transfer, which deal with the use of commutator and vertical structure equations, are described at the end of the paper. APPENDIX: SOME PROPERTIES OF PRODUCT INTEGRAL AND MATRIX EXPONENT Various properties of the product integral are used throughout this paper. We summarize some of these prop- erties hereß Refer to Dollard and Friedman [1979] for further discussion. Property 1. For a matrix A(') defined by (8) for optical thickness (-, 't) we get derivative of product integral II d- exp [A(s) ds] -- X(-) exp [X(s) ds] (149) For constant A, i.e., vertically homogeneous atmosphere, we get d X, xeX, (150) d' Property 2. For two levels defined by optical thicknesses 'b > ' > 't it follows from the deftrAtion of product integral (21) that 'b 'b H exp [A(s) ds] = H exp [A(s) ds] exp [A(s) ds] This was already discussed in some length when deriving (20). The special case of constant A and same optical thick- nesses gives (compare equation (159)) exp [A(s)r] = {exp [X(s)r/2]}' and a simple extension to m layers gives (152)

FLATAU AND STEPHENS' SOLUTION OF THE t:ADIATIVE TRANSFER EQUATION 11,049 exp IA() = exp IA()/m  (153) Property 3. Determinant of matrix exponent is given by [i ] = exp trA(s) ds (154) Where tr indicates trace (sum of diagonal elements) of a matrix. Equation (154) shows that the propagator (22) is nonsingular, i.e. der P  0. Property 4. Upward and downward propagators Pa and P, are related through matrix inversion by II exp [A(s) ds] - exp [A(s) ds] (155) In particular, for the homogenous layer, where A is constant, we have Pa = e A P,, = e -A (156) which satisfies (156). Also, A -A -A A e e = e e = I (157) Thus in terminology of radiative transfer, downward prop- agation followed by upward propagation leaves intensities unchanged. Property 5. If for all optical thicknesses % r f matrices A commute, i.e., A(r)A(r') = A(r')A(r), then the product integral of A can be evaluated as II exp [A(s) ds] = exp A(s) ds (158) In particular, tiffs formula is true for a homogeneous atmo- sphere, which gives Pt(r, r,) = exp [A(r - r,)] (159) Property 6. Certain bounds on the propagator can be established II 1-I I1_< exp II A(,)II ,a, (16o) Property 7 Similarity rule. Suppose matrix T(r) is de- fined for the layer and is nonsingular and has continuous derivative. Then T-X(r) H exp [A(s) ds] T(r) =H exp T-l (s)A(s)T(s) - T-l (r) ds which is a consequence of transformation of variables in (8). Acknowledgments. One of the authors (G.L.S.) would espe- cially like to acknowledge the contributions nd the influence of Rudolph Preisendorfer not only to the subject of radiative trans- fer in general, but to his own research in particular. The authors were aware that Rudy was actively involved in research into topics of radiative transfer that were closely related to work described in this paper before his untimely death. This work was partially supported by the NSF grant ATM-8519160 and by the Air Force grant AFOSR-88-0143. Most of the analytical calculations were done using symbolic manipulation programs: REDUCE, provided by A. C. Hearn, and the SMP. FORTRAN code of the matrix ex- ponentiation program was received from R. C. Ward. PcEFERENCES Aronson, R., ransfer-matrix solution of some one- and two- medium transport problem in slab geometry, Y. Math. Phys. N. ¾., 11,931-940, 1970. Aronson, R., Relation between the transfer matrix method and Case's method, Transp. Theory Star. Phys., 1, 209-224, 1971. Aronson, R., and D. L. Yapmush, asfer-matrix method for gamma-ray and neutron penetration, Y. Math. Phys. N.Y., 7, 221-237, 1966. Bellman, R., and R. Kalaba, Ymctional equations, wave propa- gation and invariant imbedding, Y. Math. Mech., 8, 683-704, 1959. Campbell, S. L., Singular Systems of Differential Equations, Pit- man, San I:ancisco, 1980. Case, K. M., and P. F. Zweifel, Linear Transport Theory, Addison-Wesley, Reading, Mass., 1967. Chandrasckhar, S., Radiative Transfer, Dover, New York, 1960. Coddington, E. A., and N. Levinson, Theory of Ordinary Differ- ential Equations, Mc-Graw-Hill, New York, 1955. Cole, R. H., Theory of Ordinary Differential Equations, Appleton-Century-Crofts, New York, 1968. de Hook, F., and R. Mattheij, The role of conditioning in shooting techniques, in Numerical Solution of B VODEs, Progress in Scientific Computing, Numerical Boundary Value ODEs, Proceedings of an International Workshop, edited by U. M. Ascher and R. D. Russell, Birkhauser, Boston, 1985, Vancouver, Canada. Denman, E. D., An overview of invariant imbedding algo- rithms and two-point boundary-value problems, in Codes for Boundary-Value Problems in Ordinary Differential Equa- tions, Proceedings of a Working Conference, edited by B. Childs, M. Scott, J. Daniel, E. Denman, and P. Nelson, Springer-Verlag, New York, 1979. Denman, E. D., Invariant imbedding, scattering processes and projectors in systems, Comp. math Appl., 1œA, 703-717, Dollard, J. D., and C. N. Friedman, Product integration with ap- plications to differential equations, Addison-Wesley, Read- ing, Mass., 1979. Duderstadt, J. J., and W. R. Martin, Transport Theory, John Wiley, New York, 1979. Evans, R. B., The decoupling of seismic waves, Wave Motion, 8, 321-328, 1986. Fryer, J. F., and L. N. Fazer, Seismic waves in stratified anisotropic media, Geophys. J. R. Astron. Soc., 78, 691-710, 1984. Gantmacher, F. R., Applications of the Theory of Matrices, Wiley-Interscience, New York, 1959. Gilbert, F., and G. E. Backus, Propagator matrices in elastic wave and vibration problems, Geophysics, 31, 326-332, 1966. Graham, M., Matriz Theory and Applications for Engineers and Mathematicians, Ellis Hopwood, Chichester, 1979. Grant, I. P., and G. E. Hunt, Solution of radiative transfer prob- lems in planetary atmospheres, Icarus, 9, 526-534, 1968a. Grant, L P., and G. E. Hunt, Solution of radiative transfer prob- lems using the invariant Sn method, Mon. Not. R. Astron. Soc., 141, 27-41, 1968b. Grant, I. P., and G. E. Hunt, Discrete space theory of radiative transfer, I. Fundamentals, Proc. R. Soc. London Set. A, 313, 183-197, 1969a. Grant, I. P., and G. E. Hunt, Discrete space theory of radiative transfer, II. Stability and non-negativity, Proc. R. Soc. Lon- don Set. A, 313, 199-216, 1969b.

11,050 FLATAU AND STEPHENS: SOLUTION OF THE lqADIATIVE TBANSFEI{ EQUATION Karp, A. H., J. Grecnstadt, and J. A. Filmore, Radiative trans= fer through an arbitrarily thick, scattering atmosphere, J. Quant. Spectrosc. Radiat. Transfer, 24, 391-406, 1980. Keller, H. B., Approximate solutions of transport problems, II, Convergence and applications of the discrete ordinate method, J. Soc. Ind. Appl. Math., 8, 43-73, 1960. Kermctt, B. L. N., Reflections, rays and reverberations, Bull. Seismol. Soc. Am., 64, 1685-1696, 1974. Kcrmett, B. L. N., and N.J. Kerry, Seismic waves in a strathqed half space, Geophys. J. R. astron. Soc., 57, 557-583, 1979. Lcnobl½, J. (ed.), Radiative Transfer in Scattering and Absorbing Atmospheres: Standard Computational Procedures, (Stud- ies in Geophysical Optics and Remote Sensing), A. Deepak, Hampton, Va., 1985. Liou, K., Application of the discrete-ordinate method for radia- tive transfer to inhomogeneous aerosol atmospheres, J. Geo- phys. Res., 80, 3434-3440, 1975. Magnus, W., On the exponential solution of differential equations for a linear operator, Comm. Put. Appl. Math., 7, 423-673, 1954. Matth½ij, R. M. M., Decoupling and stability of algorithms for boundary value problems, SIAM Rev., 27, 1-44, 1985. Matthe[i, R. M. M., and G. W. M. Staarink, An efficient algo- rithm for solving general linear two-point BVP, SIAM J. Sci. Star. Cornput., 5, 745-763, 1984. Meador, W. E., and W. R. Weaver, Two-stream approximation to radiative transfer in planetary atmospheres: A unified de- scription of existing methods and a new improvement, J. Attooz. Sci., 37, 630-643, 1980. Molcr, C. B., and C. van Loan, Nineteen dubious ways to compute the exponential of a matrix, SIAM Rev., 4,801-836, 1978. Mulholland, R. J., Exponential representation for linear systems, IEEE Tans. Aurora. Control, I6, 97-98, 1971. Nakajima, T., and M. Tanaka, Matrix formulations for the trans- fer of solar radiation in a plane-parallel scattering atmo- sphere, J. Quant. Spectrosc. Radiat. Transfer, 35, 13-21, 1986. Peebles, G. H., and M. Plesset, ransmission of gamma rays through large thicknesses of heavy metal, Phys. Rev., 8I, 430-439, 1951. Pennisi, L. L., Elements of Ordinary Differntial Equations, Holt, Rinehart and Winston, New York, 1972. Pollack, J. B., and C. P. McKay, The impact of polar stratospheric clouds on the heating rates of the winter polar stratosphere, J. Atmos. Sci., 42, 245-262, 1985. Preisendorfer, R. W., Invariant imbedding relation for the princi- ples of invariance, Proc. Natl. Acad. Sci., USA, 44, 320-323, 1958. Preisendorfer, R. W., Classic canal theory, Technical Report, Hawaii Inst. Geophys., Univ. of Hawaii, Honolulu, 1973. Preiscndorf'cr, R. W., Eigenmatr epresentation of radiance dis- tribution in layered natural waters with wind-uffled sur- faces, Technical Report, Pacific Marine Environmental Lab- oratory, NOAA, Seattle, WA., 1988, Contribution No. 841. Redheifer, R. M., On the relation of transmittion-line theory to scattering and transfer, J. Math. Phys. N.Y., 4I, 1-41, 1962. Russell, R., A unified view of some recent developments in numer- ical solution of bvodes, in Progress in Scientific Computing, Numerical Boundary Value ODEs, Proceedings of an Inter national Workshop, edited by U. M. Ascher and R. Russell, Birkhauser, Boston, 1985, pp 1-19, Vancouver, Canada. Saastamoinen, P. R., On propagators and scatterers in wave prob- lems of layered, elastic media- a spectral approach, Bull. Seismol. Soc. Am., 4, 1125-1135, 1980. Sautcr, E.G., Improvement on the WKB method using an al- gebraic formalism, J. Phys. A Math. Gen., I2, 1155-1164, 1979. Scott, M. R., On the Conversion of Boundary-Layer Problems Into Stable Initial-Value Problems via Several Invariant Imbedding Algorithms, Academic, San Diego, Calif., 1975. Shettle, E., and J. A. Weinman, The transfer of solar irradi- ance through inhomogeneous turbid atmospheres evaluated by Eddington's approximation, J. Atmos. Sci., 27, 1048- 1055, 1970. Smith, D. R., Decoupling and order reduction via Riccati trans- formation, SIAM Rev., 29, 91-113, 1987. Steinberg, S., Lie series, Lie transformations and their applica- tions, 1987, in Lie Methods in Optics, Lect. Notes in Phys., Springer-Verlag, New York, in press. Starnnes, K., The theory of multiple scattering of radiation in plane parallel atmospheres, Rev. Geophys., 24, 299-310, 1986. Stamnes, K., and R. A. Swanson, A new look at the discrete ordi- nate method for radiative transfer calculations in anisotrop- ically scattering atmospheres, J. Atmos. Sci., 38, 387-399, 1981. Stephens, G. L., Radiative transfer in spatJelly heterogeneous two-dimensional anisotropically scattering media, J. Quant. Spectrosc. Radiat. Transfer, 36, 51-67, 1986. Stephens, G. L., Radiative transfer through arbitrarily shaped op- tical media, I, A general method of solution, 1988, J. Atmos. Sci., in press, 1988. van Loan, C. F., Computing integrals involving the matrix expo- nential, IEEE Trans. Aurora. Control, 23, 395-404, 1978. van Loon, P., Riccati transformation: When and how to use ?, in Numerical solution of B VODEs, Progress in Scientific Computing, Numerical Boundary Value ODEs, Proceedings of an International Workshop, edited by U. M. Aschcr and R. D. Russell, Birkhauser, Boston, 1985, Vancouver, Canada. Ward, R. C., Numerical computation of the matrix exponential with accuracy estimate, SIAM J. Numer Anal., I4, 600- 609, 1977. Waterman, P. C., Matrix-exponential description of radiative transfer, J. Opt. Soc. Am., 71,410-422, 1981. Wei, J., and E. Norman, On global representation of the solutions of linear differential equations, Proc. Am. Math. Soc., 15, 327-334, 1964. Wilcox, R. M., Exponential operators and parameter differentia- tion in quantum physics, J. Math. Phys. N.Y., 8, 962-982, 1967. Wiscombc, W. J., On initialization, error and flux conservation in the doubling method, J. Quant. Spectrosc. Radiat. Transfer, 16, 637-658, 1976. P. J. Flatau and G. Stephens, Department of Atmospheric Sci= ence, Colorado State University, Fort Collins, CO 80523. (Received August 10, 1987; revised March 28, 1988; accepted April 6, 1988.)