The analysis of problems of radiative and particle transport provides us with a vast literature with history going well back into the nineteenth century. Many well known names are associated with this development: Stokes, Rayleigh, Eddington, Milne, Hopf, Ambartsumian, and Chandrasekhar, to mention only a few. The subject is a difficult one; classical methods aim to manipulate the equations into the form of an integral equation whose solution can only be obtained analytically if the most restrictive assumptions of homogeneity of the region of interest are made. Even then, when the results have been expressed in terms of tabulated standard functions, practical computations are difficult and tedious.
Numerical techniques have been applied in transport problems, the usual approach being to approximate a differential or integral equation as a system of algebraic equations. For the most part, the theoretical understanding of the convergence, stability, and accuracy so necessary in making practical calculations has been lacking. This is because of the peculiar unwieldiness of the equations when scattering terms are present. Thus the only methods that have been available to check the validity of numerical computations have been heuristic in nature.
The methods to be outlined here had their origin in dissatisfaction with current numerical techniques for solving radiative transfer problems both from a theoretical and a practical point of view. The feeling that, since the equations of radiative transfer express the conservation of radiant energy, any accurate numerical technique must also express this conservation was another motive for our investigation. The correctness of this conjecture has been shown by subsequent work as we shall describe.
For simplicity, we shall confine our treatment to the case of steady monochromatic radiative transfer in plane-parallel media. In itself, the mathematics of this case is already quite complex. A thorough understanding of this problem is essential before going on to more complicated situations.
The situation we wish to study is depicted in Figure 1. We shall suppose, in accordance with observation, that the radiation emerging from the slab, represented by vectors u2 +, u1 -, is related linearly to the radiation incident on its bounding surfaces x = x1, x = x2, x2 > x1, represented by vectors u1 +, u1 - and to the radiation generated within the slab. We express this by writing
...(1)
The elements of the vectors ui +, ui - are the intensities in m discrete directions making angles with respect to the direction of x increasing whose cosines may be denoted by μj, -μj, j=1,2,...m, 0 < μj ≤ 1, thus
...(2)
The coefficients of u1 +, u2 - in (1) are therefore m × m matrices which may be interpreted as operators for diffuse transmission (t) and diffuse reflection (r). The effective sources within the slab are denoted by Σ½ + and Σ½ -, vectors which are to be interpreted in the same fashion as u+ and u- vectors.
Equation (1) contains the entire physics of the problem. It is an exact relation, valid whether the slab is homogeneous or inhomogeneous, scattering or non-scattering, optically thick or optically thin. However, to make use of it, we must be able to compute the r and t matrices and source vectors for slabs of arbitrary dimensions. To begin with, we wish equations (1) to pass over in the limit to the usual equation of transfer
...(3)
where 0 < μ ≤ 1, 0 ≤ ≤ 1, and B+(x), B-(x), p(x;μ,μ') are all prescribed functions of their arguments. The quantity , usually called the albedo for single scattering, gives the expected fraction of radiant energy scattered at interactions with the material, and p(x;μ,μ') describes how the direction is changed by scattering. The phase function p(x;μ,μ') is normalized so that
and it usually satisfies the symmetry condition
p(x; μ,μ') = p(x; -μ,-μ') ...(4)
In discrete ordinate methods we approximate to integrals over (say) [0,1] by writing them as Gaussian sums in the form
...(5)
where {μj}, and {cj} are some set of Gaussian division points and weights and εm(f) is an expression for the error. Setting μ equal to each of the m values of μj in turn, we may write equations (3) in the matrix form
...(6)
where M and c are diagonal m × m matrices
M=[μjδjk], c=[cjδjk]
and
P+(x)=[p(x;μj, μk)], P-(x)=[p(x;-μj, +μk)]
are symmetric m × m matrices, whilst u+ (x), u-(x) are defined as in (2). It is now easy to see that when x2 - x1 → 0, equations (1) converge to equations (6) provided
...(7)
whilst
...(8)
where x1≤ x3/2 ≤ x2. These conditions are necessary, but one prefers to have expressions which possess a higher order truncation error in practice, and these are available [1].
Equations (7) give us the response of a thin finite slab to an arbitrary external flux. We can build thicker slabs by superposition- adding layers is the traditional term - and so if we write (1) in the form
and likewise
the response of the composite system will be given by an equation of the form
All that is necessary is to eliminate u2 +, u2 - between the four equations. The matrix S(1, 3) may be defined in terms of the star product
S(1, 3) = S(1, 2)*S(2, 3). .......(9)
Writing the S-matrices in partitioned form,
we can express the star-product through the relations
t(3,1) = t(3,2) [I - r(1,2)r(3,2)]-1 t(2,1), t(l,3) = t(l,2) [I - r(3,2)r(l,2)]-1 t(2,3), r(3,1) = r(2,1) + t(l,2)r(3,2)[I - r(1,2)r(3,2)]-1 t(2,1), ......(10) r(1,3) = r(2,3) + t(3,2)r(l,2)[I - r(3,2)r(l,2)]-1 t(2, 3),
expressions which are well-known from transmission-line theory [2]. Also if
then
Σ+(3,1) = Σ+(3,2) + t(3,2) [I - r(1,2) r(3,2)]-1[Σ+(2,1) + r(1,2)Σ-(2,3)] ....(11) Σ1(1,3) = Σ-(1,2) + t(1,2) [I - r(3,2) r(1,2)]-1[Σ-(2,3) + r(3,2)Σ+(2,1)]
define
Σ(1,3) = Σ(1,2)*Σ(2,3)
The relations (11) appear to have been ignored up till now.
The matrices S(m,n) form a semigroup with respect to the star product: however we shall make little use of this fact in the sequel.
For homogeneous media, it is clear that S(m,n) depends only on the thickness of the layer xn-xm. Since, if we superpose two layers of thickness x and y, the order in which the layers are taken cannot matter, equation (9) will give
S(x + y) = S(x)*S(y) = S(y)*S(x),
which is Stokes' commutativity property [3]. If, in this relation, we let y → ∞, we get Ambartsumian's principles of invariance [4],
S(∞) = S(x)*S(∞) = S(∞)*S(x)
valid for any finite x.
In the general inhomogeneous case, we may insert equations (7) into (9) and (10) and, if we take the limit as x2-x1 → 0, we obtain a system of differential equations expressing Chandrasekhar's principles of in variance for finite slabs [5]. Most of the numerical work that has been done on the external response of scattering media has been directed towards the solution of Chandrasekhar's equations. More recently, a number of people, particularly van de Hulst and his collaborators, have realized that there is no difficulty in applying equations (10) once an initial approximation is available for thin slabs [6]. These equations can be used in two ways: first, by straightforward addition of successive layers; second, for homogeneous slabs, by using (10) to give the response of a double layer of thickness 2d from two layers of thickness d, followed by repetition to give successively matrices for thicknesses 22d, 23d,... ,2pd in p cycles. This latter technique is particularly useful in certain types of practical problem, for example the study of non-coherent scattering [10].
We have seen in the last section how to calculate matrices and sources for finite layers by building them up from layers of smaller thickness. In doing so, we deliberately eliminated the internal fluxes which we now wish to compute. Suppose that we have a slab divided into layers by planes at positions x1, x2,...,xN+1 = a. Each layer or cell between the planes at xn and xn+1 will be associated with an equation of the form
...(12)
We have shown, [1], that the complete solution of this system of equations can be obtained recursively as follows:
For n = 1,2,..., N, compute successively
r(1,1) = 0, r(1,n+1) = r(n,n+1).
+t(n+1,n)r(1,n)[I - r(n+l,n)r(l,n)]-1t(n,n+1) ...(13)
and
...(14)
where
Then, with u- N+1 given, compute successively
...(15)
This system of equations is particularly simple and rapid to compute in practice.
If, instead of having uN+1 - given, we have a reflecting boundary, as at the base of a planetary atmosphere, we have a simple modification. Suppose that the reflection may be described by a matrix rG so that
u-
N+1 = rGu+
N+1 ...(16)
Then (15) gives
u-
N+1=[I-r(1,N+1)rG]-1V+
N+½
and the remaining internal fluxes can be computed using (15) as before. The effect on external response can be described by the relation
SG(l,2) = S(1,2)*SG
where
and S(1,2) is the matrix for the slab without the reflecting termination.
If we use the limiting form (7) for the elements of S(n,n+1), and pass to the limit as xn+1 - xn → 0, we obtain the differential equation system first suggested by Rybicki and Usher [7], which in the present notation becomes
...(17)
and
u+(x) = r(x) u-(x) + V+(x).
Here
and
as in (6)-(8).
Standard techniques for numerical solution of initial-value problems can be used to solve (17).
Before proceeding to numerical calculations we should try to assure ourselves that the equations possess a unique solution. Also that the algorithm used to generate the solutions should be stable, and not allow errors to build up. The physical specific intensity is certainly non-negative, and we want to ensure that the mathematical solution shares this property.
For the purpose of this analysis, we define a radiometric norm; for example the norm of a vector u+ will be taken to be
...(18)
which is an approximation to the net positive flux,
when u+ has non-negative elements. The matrix norm of an operator Q consistent and subordinate with this definition is
...(19)
In applying this to equations (1) with u2 +, u1 - having non-negative elements (we write u2 + ≥ 0, u1 - ≥ 0), we find if
then
...(20)
Using (7), we see that, when
... (21)
and that, for x2-x1 sufficiently small, S(1,2) is a non-negative matrix, S(1, 2) ≥ 0. It then follows that matrices S(m,n) for finite layers, built up from elementary non-negative cell matrices with norm less than unity by means of star products, are also non-negative and have norm less then unity. In particular submatrices t(m,n), r(m,n) are non-negative, and have norms strictly less than unity. These conditions are sufficient so that the solutions of (12) (or (13)-(15)) exist, are non-negative and are uniquely determined [8].
Taken with the knowledge that the finite difference equations (12) are consistent in the limit with the conventional differential equations, these results give us confidence in our methods and the numerical solutions to which they give rise. It is, however, possible to say more. Suppose we now restrict ourselves to a homogeneous slab, divided into N identical layers. Then, if by A ≤ B we mean that each matrix element aij ≤ bij, we have
0 ≤ r(1,2) ≤ r(1,3) ≤ . .. ≤ r∞ ) = sup r(1, n).
The supremum exists since
for all n. Thus
Now for i < j
If the maximum eigenvalue of is denoted by λ ≥ 0, then it can be proved that, for k sufficiently large, there exists a positive constant v such that
Thus, for any i ≠ j,
.....(22)
We have proved that this is strictly less than unity, and so λ < 1 also. From this we get an estimate
.....(23)
for the norm of the reflection matrix, which is uniformly bounded with respect to n for λ < 1.
If the basic cell has thickness h, then for h → 0,
and so there exists a constant γ > 0 such that
Keeping a = Nh fixed, and letting N → ∞, we find
...(24)
for transmission and reflection of slabs of thickness a. Similar bounds can be computed for norms of fluxes at a point, but these are not very informative.
Although the equations of the method, (1), (9) and (12) are formally exact, the matrices S(n,n+1) and source vectors Σ(n,n+1) which must be used in practice will be correct only to some given order in layer thickness h. Indeed it is only necessary that (7) and (8) should hold, and then we can express the local truncation error in the form
where k ≥ 2, for some matrices θ, ρ and vectors σ±. We then find that perturbations of t(a) and r(a) are bounded by
and
which go to zero as hk- 1 for fixed a. It can also be shown that the relative error in internal fluxes can be written
for some constant C.
It is therefore important to choose matrices S(n,n+1) and source Σ(n,n+1) which are correct to second order in h if possible. Expressions for which k = 3 are given in [1], so that the resulting errors converge quadratically to zero as h → 0.
We have only dealt with the spatial truncation error. The error εm(f) arising from the Gaussian quadrature over μ, equation (5), could also be studied if we had a suitable expression for this error. Standard treatments express Gaussian quadrature errors either in terms of a high order derivative of the integrand, or in terms of its values at points on a contour in the complex μ-plane [9]. These do not seem particularly suitable to apply to functions defined only at the Gaussian quadrature points, and it is clear that further work is required on this aspect.
We have now done a number of computations, both of external response functions and of internal fluxes when internal emission sources are present. Most of these computations have had the aim of checking the correctness of the theory and of providing figures for comparison with other published calculations. The data available for this purpose mainly relate to reflection from isotropically scattering homogeneous atmospheres, computed mainly by differential equation methods based on the imbedding equations of Chandrasekhar [5]. Using seven Gaussian points for the angular resolution, we have found good agreement (5 or 6 figures) with computations of Bellman et al. [11]. Some of our results appear in [1].
Calculations relating to scattering by planetary atmospheres with reflecting lower boundaries appear in [12]. In addition we have extended the programs to handle polarization in Rayleigh scattering problems, although we have not as yet included all the Stokes parameters [13]. The results appear to be consistent with figures given by Chandrasekhar [5].
It is less easy to compare calculations relating to internal light fields. We have done calculations both with our own equations (13)-(15) and with the Rybicki-Usher equations (17) and obtained good agreement. In general, our method is more flexible, since the techniques of Section 3 allow us to circumvent restrictions on layer thickness by way of the star product. Without this, we would have to ensure that the optical thickness h would not be so large that the matrices defined by (7) failed to be nonnegative:
For the case of the Rybicki-Usher method stability of the numerical solution requires a step size
for most methods. Thus we can handle a layer of given thickness with fewer subdivisions. This is of particular importance in non-coherent line formation calculations [1O], where the atmosphere may be optically thin in the wings of the line, but extremely thick (about 103, say) in the core. The coupling between frequencies in different parts of the line then presents a problem which is not difficult to handle using our methods but is much harder to solve using the Rybicki-Usher equations.
A practical problem for which we have used equations (l3)-(15) is related to a joint proposal by atmospheric physicists at Oxford and Reading Universities to observe the global temperature structure of the terrestrial atmosphere at altitudes up to 50 km over an extended period of time by means of a selective chopper radiometer mounted on the Nimbus-E weather satellite [14]. As part of this proposal, it is necessary to make supporting observations of water vapour distribution and to try to infer something about the density of cirrus cloud. We have made some preliminary calculations on what signals might be expected for radiation in the 8-150 μ region due to emission from molecular bands of water in the presence of cirrus cloud. We have used a model in which the water vapour band is represented by an equivalent single line giving the same curve of growth. The data were provided by Dr. C. D. Rodgers of Oxford University. Scattering within the line is assumed to be coherent. Results have suggested that the experiment is feasible; changes in the signal corresponding to reductions in blackbody temperature from 1 or 2°C to about 20°C in different cloud models have been found. These calculations could not so easily have been made by other methods.
The methods described in this paper provide us with a powerful new weapon for attacking transport problems. Not only can we now circumvent many step size difficulties that plagued us in the past. We also have available, for the first time, a penetrating insight into the structure of the available numerical methods. The new discrete methods are simple to program, and compete well in terms of computer storage and execution time with rival techniques [1].
There is still much to be done in the way of development. We have presented here a theory for anisotropic scattering restricted to axi-symmetric phase functions. We are endeavouring at present to find an economic technique for handling the general anisotropic case. We also want to be able to handle the full problem of polarization including all the Stokes parameters. We expect our method to be of great assistance in non-coherent scattering problems, and we intend to verify this. Finally, there are many problems, particularly in neutron transport theory, where it is necessary to work in coordinate systems with less symmetry than the plane parallel slab-spherical coordinates, cylindrical and so on. There is no reason to think that our treatment will not generalize to such cases.
This work has been done with the active collaboration of Dr. G. E. Hunt of Reading University. Mrs. E. Gill has assisted with the programming.
1. Grant, I. P. and Hunt, G. E. 1968a. Mon. Not. R. Astr. Soc. 141, 27--48.
2. Redheffer, R.M. 1962. J. Math. Phys. 41, 1-41.
3. Stokes, G. G. 1862. Proc. Roy. Soc. 11, 545-556.
4. Ambartsumian, V. A. 1943. Dokl. Akad. Nauk SSSR 38, 229-232.
5. Chandrasekhar, S. 1950. 'Radiative Transfer,' Ch. VII. O.U.P.
6. van de Hulst, H. C. 1962. 'A New Look at Multiple Scattering' (NASA, Institute for Space Studies, New York); also H. C. van de Hulst and K. Grossman, 1968.. in 'Atmospheres of Mars and Venus' (edited by J. C. Brandt and M. B. McElroy) pp. 35-55. Gordon and Breach Science Publishers, New York.
7. Rybicki, G. B. and Usher, P. D. 1966. Ap. J. 146, 871-879.
8. Grant, I. P. 1968. Unpublished material.
9. Stroud, A. H. and Secrest, D. 1966. 'Gaussian Quadrature Formulas,' Ch. 4. Prentice Hall
10. Hummer, D. G. and Rybicki, G. B. 1967. In 'Methods of Computational Physics,' Vol 7 (ed. B. Alder et al.) Academic Press.
11. Bellman, R., Kalaba, R. and Prestrud, M. 1963. 'Invariant Imbedding and Radiative Transfer in Slabs of Finite Thickness.' Elsevier, New York.
12. Grant, I. P. and Hunt, G. E. 1968b. Icarus, submitted for publication.
13. Grant, I. P. and Hunt, G. E. 1968c. J. Quantit. Spectrosc. rad. Transfer 8,1817-1832.
14. Houghton, J. T., et al. 'Selective Chopper Radiometer for Water-Vapour, Cloud and Atmospheric Temperature Sounding for the Nimbus-E Satellite.' (Scientific section of joint proposal to NASA by Oxford and Reading Universities, 4th March 1968.)