Modified spherical harmonics method for solving the radiative transport equation

A new effective approach to solving the three-dimensional radiative transport equation with an arbitrary phase function is proposed. The solution depends on eigenvectors and eigenvalues of several symmetrical tridiagonal matrices of infinite size. The matrices must be truncated and diagonalized numerically. Then, given eigenvectors and eigenvalues of these matrices, the dependence of the solution on position and direction is found analytically. The approach is based on expanding the angular part of the specific intensity in q-dependent spherical functions for each spatial Fourier component characterized by the vector q. Apart from the truncation of the matrices, no other approximations are made. The radiative transport equation (RTE) is widely used to describe propagation of waves and particles in random scattering media. Applications include biomedical imaging of tissues with near-IR light [1, 2], propagation of waves in atmosphere and oceans and astrophysics [3–5], and nuclear reactor theory [6]. However, analytical solutions to the RTE cannot be obtained even in very simple cases. In fact, the only known analytical solution to the three-dimensional RTE is for the case of an infinite homogeneous medium and a constant phase function (isotropic scattering) [5, 7]. Even in this case the solution is expressed as a quadrature which must be evaluated numerically. Therefore, current research into methods of solving the RTE is mainly focused on numerical or approximate methods (see, for example, [8–10]). When analytical solutions are unavailable, a number of numerical methods can be applied, such as the discrete ordinate method or the method of spherical harmonics [6]. In particular, the method of spherical harmonics can be effectively used in the cases with special symmetry, such as an isotropic source in the infinite medium (spherical symmetry) or a plane wave incident on a slab or half-space (cylindrical symmetry). In a more general situation, when none of the above symmetries are present, the mathematical formulation of the standard spherical harmonics method is very complicated and hardly usable, except in conjunction with the Pl approximation [6]. In this letter, I derive a modification of the spherical harmonics method. The modification leads to significant mathematical simplifications and can be used to calculate the 0959-7174/04/010013+07$30.00 © 2004 IOP Publishing Ltd Printed in the UK L13 L14 Letter to the Editor Green’s function of the RTE for a point unidirectional source in infinite homogeneous medium. The obtained solutions depend on eigenvectors and eigenvalues of certain tridiagonal matrices which, in turn, depend only on the form of the phase function. After these quantities are computed numerically, the dependence of solutions on the position and direction of propagation is found analytically. The approach developed below can be applied to the Boltzmann equation describing transport of any waves or particles that move with a fixed absolute velocity c but can change direction as a result of collisions or scattering events. To be more specific, we will assume that the RTE describes propagation of light in a multiply scattering medium. In this case c is the average velocity of light in this medium (set to unity everywhere below) and the quantity of interest is the specific intensity I (r, ŝ) at the point r and in the direction specified by the unit vector ŝ. The RTE has the form ŝ · ∇ I + (μa + μs)I = μs  I + ε. (1) Here μa and μs are the absorption and scattering coefficients, respectively,  is an integral operator defined by  I (r, ŝ) = ∫ A(ŝ, ŝ′)I (r, ŝ′) d2ŝ′, (2) A(ŝ, ŝ′) is the phase function and,finally, ε = ε(r, ŝ) is the source. We limit consideration to the case when the phase function depends only on the angle between ŝ and ŝ′: A(ŝ, ŝ′) = f (ŝ · ŝ′). This fundamental assumption is often used and corresponds to scattering by spherically symmetrical particles. We also require that the phase function is normalized by the condition ∫ A(ŝ, ŝ′) d2ŝ′ = 1. We start with expressing all position-dependent functions as Fourier integrals, according to I (r, ŝ) = ∫ Ĩ (q, ŝ) exp(iq · r) dq, (3) ε(r, ŝ) = ∫ ε̃(q, ŝ) exp(iq · r) dq. (4) Substituting (3) and (4) into (1), we obtain (iq · ŝ + μt) Ĩ = μs  Ĩ + ε̃, (5) where we have introduced the notation μt = μa + μs. Next we expand all angular-dependentquantities in spherical harmonics which are defined in the reference frame whose z-axis coincides with the direction of the vector q. We denote such functions as Ylm(ŝ; q̂) and write Ĩ (q, ŝ) = ∑ lm Flm(q)Ylm(ŝ; q̂), (6)