Eigen decomposition solution to the one-dimensional time-dependent photon transport equation

The time-dependent one-dimensional photon transport (radiative transfer) equation is widely used to model light propagation through turbid media with a slab geometry, in a vast number of disciplines. Several numerical and semi-analytical techniques are available to accurately solve this equation. In this work we propose a novel efficient solution technique based on eigen decomposition of the vectorized version of the photon transport equation. Using clever transformations, the four variable integrodifferential equation is reduced to a set of first order ordinary differential equations using a combination of a spectral method and the discrete ordinates method. An eigen decomposition approach is then utilized to obtain the closed-form solution of this reduced set of ordinary differential equations. © 2011 Optical Society of America OCIS codes: (010.5620) Radiative transfer, (290.7050) Turbid media, (080.2720) Mathematical methods. References and links 1. A. D. Klose, V. Ntziachristos, and A. H. Hielscher, “The inverse source problem based on the radiative transfer equation in optical molecular imaging,” J. Comput. Phys. 202, 323–345 (2005). 2. N. Y. Gnedin, “Multi-dimensional cosmological radiative transfer with a variable Eddington tensor formalism,” New Astron.6, 437–455 (2001). 3. D. M. O’Brien, “Accelerated quasi Monte Carlo integration of the radiative transfer equation,” J. Quant. Spectrosc. Radiat. Transf. 48, 41–59 (1992). 4. K. Stamnes, S. C. Tsay, W. Wiscombe, and K. Jayaweera, “Numerically stable algorithm for discrete-ordinatemethod radiative transfer in multiple scattering and emitting layered media,” Appl. Opt. 27, 2502–2509 (1988). 5. J. L. Deuz, M. Herman, and R. Santer, “Fourier series expansion of the transfer equation in the atmosphere-ocean system,” J. Quant. Spectrosc. Radiat. Transf. 41, 483–494 (1989). 6. C. E. Siewert and J. R. Thomas, “Radiative transfer calculations in spheres and cylinders,” J. Quant. Spectrosc. Radiat. Transf. 34, 59–64 (1985). 7. C. E. Siewert, “A radiative-transfer inverse-source problem for a sphere,” J. Quant. Spectrosc. Radiat. Transf. 52, 157–160 (1994). 8. E. W. Larsen, “The inverse source problem in radiative transfer,” J. Quant. Spectrosc. Radiat. Transf. 15, 1–5 (1975). 9. G. C. Pomraning and G. M. Foglesong, “Transport-diffusion interfaces in radiative transfer,” J. Comput. Phys. 32, 420–436 (1979). 10. Z. M. Tan and P. F. Hsu, “An integral formulation of transient radiative transfer,” ASME J. Heat Transfer 123, 466–475 (2001). #137125 $15.00 USD Received 25 Oct 2010; revised 27 Jan 2011; accepted 28 Jan 2011; published 1 Feb 2011 (C) 2011 OSA 14 February 2011 / Vol. 19, No. 4 / OPTICS EXPRESS 2922 11. J. A. Fleck, “The calculation of nonlinear radiation transport by a Monte Carlo method,” Methods Comput. Phys. 1, 43–65 (1963). 12. A. D. Kim and A. Ishimaru, “A Chebyshev spectral method for radiative transfer equations applied to electromagnetic wave propagation and scattering in discrete random media,” J. Comput. Phys. 152, 264–280 (1999). 13. A. D. Kim and M. Moscoso, “Chebyshev spectral methods for radiative transfer,” SIAM J. Sci. Comput. (USA) 23, 2074–2094 (2002). 14. C. C. Handapangoda, M. Premaratne, L. Yeo, and J. Friend, “Laguerre Runge-Kutta-Fehlberg method for simulating laser pulse propagation in biological tissue,” IEEE J. Sel. Top. Quantum Electron. 14, 105–112 (2008). 15. A. D. Kim and M. Moscoso, “Chebyshev spectral methods fro radiative transfer,” SIAM J. Sci. Comput. (USA) 23, 2074–2094 (2002). 16. A. B. Carlson,Communication Systems: An Introduction to Signals and Noise in Electrical Communication (McGraw-Hill, 1986). 17. S. Chandrasekhar, Radiative Transfer (Dover Publications, 1960). 18. K. Stamnes and R. A. Swanson, “A new look at the discrete ordinate method for radiative transfer calculations in anisotropically scattering atmospheres,” J. Atmos. Sci. 38, 387–399 (1981). 19. W. H. Press, S. A. Teukolsk, W. T. Vetterling, and B. P. Flannery, Numerical Recipes in C++, 2nd ed. (Cambridge University Press, 2003). 20. C. H. Edwards and D. E. Penney, Differential Equations: Computing and Modeling, 3rd ed. (Prentice Hall, 2004). 21. C. C. Handapangoda and M. Premaratne, “Implicitly causality enforced solution of multidimensional transient photon transport equation,” Opt. Express 17, 23423–23442 (2009).


Introduction
The photon transport equation (PTE), also known as the radiative transfer equation (RTE), is used in various different disciplines to model light propagation through turbid media and to characterize the scattering properties of media.For example, biomedical applications such as near-infrared fluorescence tomography [1], astrophysical and cosmological applications such as problems related to the reionization of the intergalactic medium and absorption line signatures of high redshift structures [2], meteorological applications such as visibility studies (e.g.modeling haze) and modeling shortwave and longwave fluxes transmitted by fields of broken clouds [3], all require the solution to the radiative transfer equation.
A number of models for solving the steady state (time-independent) PTE has been developed over the past four decades [4,5], ignoring the time dependency of the intensity profile.Some of these models include the discrete ordinates method [4], integral transformation techniques and the F N method [6].Siewert [7] and Larsen et al. [8] have worked on the inverse-source problem where the source term is determined from the known angular distribution of radiation that exits the surface.Pomraning et al. [9] have coupled the diffusion approximation and the radiative transfer equation for boundary treatment for the steady state case.
With the advent of short laser pulses that are applied in various different disciplines, researchers have started developing techniques to solve the transient (time-dependent) PTE.Tan and Hsu [10] developed an integral equation formulation to treat the general transient PTE, and, Fleck used the Monte Carlo simulations [11].Kim, Ishimaru and Moscoso have applied Chebyshev collocation techniques for solving the PTE [12,13].Kim [13] has used Chebyshev expansion for the spatial discretization and Crank-Nicholson method for time marching.Handapangoda et al. [14] proposed a technique to transform the four-variable transient PTE to a set of coupled first order ordinary differential equations.However, in that work they solved this set of equations using the Runge-Kutta-Fehlberg method, which is a numerical technique.In this paper we propose an improvement to this technique, which gives a closed-form solution to the transformed PTE based on eigen decomposition.This approach improves the computational efficiency significantly as shown later.
This paper is organized as follows.In section 2 we present the eigen decomposition method to solve the reduced vectorized version of the PTE.Section 3 presents some numerical results.A comparison of the output values and the execution time of the algorithm is carried out against a those of the standard solution, which is obtained using the Laguerre Runge-Kutta-Fehlberg (LRKF) method [14].Section 4 carries some concluding remarks.

Eigen decomposition method
The one dimensional transient photon transport equation for a medium without any internal sources is given by [14,15] where I (z, u, φ ,t) is the light intensity (radiance), (z, θ , φ ) are the standard spherical coordinates, u = cos θ , t denotes time, σ t and σ s are attenuation and scattering coefficients, respectively.The speed of light in the medium is denoted by v and P (u ′ , φ ′ ; u, φ ) is the phase function.Suppose Eq. ( 1) is subject to the boundary condition , where δ (x) is the Dirac delta function [16] and f (t) is the temporal profile of the incident pulse.Using the transformation which maps Eq. ( 1) to a moving reference frame with the pulse, we get 2) in Eq. ( 1) eliminates the time derivative term in Eq. ( 1).Then, the azimuthal angle, φ , can be discretized using the discrete ordinates method [4,17,18] by applying the Gaussian quadrature rule [19] for the integral that corresponds to φ .
Gaussian quadrature [19] is used to approximate an integral of the form b a W (x) f (x) dx, where W (x) is the weight function and f (x) is any arbitrary function, by a summation such that b a W (x) f (x) dx ≈ ∑ N−1 j=0 w j f (x j ).Application of the discrete ordinates method to the azimuthal angle of Eq. (3) will result in a set of uncoupled equations where r = 1, . . ., L, φ j is the j th quadrature point of φ and w φ j is the corresponding Gaussian weight.We then discretize the cosine of the zenith angle, u, of Eq. ( 4) using the discrete ordinates method, which results in where i = 1, . . ., K, r = 1, . . ., L, u k is the k th quadrature point of u and w u k is the corresponding Gaussian weight.
In order to remove the dependency on τ, we can expand I (z, u i , φ r , τ) and f (τ) using Laguerre polynomials with respect to τ, such that where B k F k are the coefficients corresponding to the Laguerre polynomial L k (τ).Expanding the transformed one-dimensional PTE, Eq. ( 5) and the boundary condition using a Laguerre basis and taking moments we get, subject to where n = 1, . . ., N, i = 1, . . ., K and r = 1, . . ., L. Thus, there are K × L coupled equations for each Laguerre coefficient, B n and Eq. ( 6) corresponds to the i th discrete ordinate of u and r th discrete ordinate of φ .We can write this set of equations in matrix form as; where B n = [B n (z, u i , φ r )] K×L,1 , P = [P (u k , φ j ; u i , φ r )] K×L,K×L , and, A is a (K × L) by (K × L) diagonal matrix with diagonal elements u 1 to u K repeating L times.The matrix W is also a Rearranging, Eq. ( 8) can be written as where Y = A −1 σ s 4π PW − σ t I and I denotes the identity matrix.Hence, the original transient PTE is reduced to a one-variable ordinary differential equation.The boundary condition given by Eq. ( 7) can be decomposed into a set of values on the boundary corresponding to each discrete ordinate, which can be represented in a column matrix as B 0 n = [B n (z = 0, u i , φ r )] K×L,1 .
An explicit solution to Eq. ( 9) can be found using eigen decomposition.The square matrix Y can be written as a product of three matrices composed of its eigen values and eigen vectors, as follows. where . v i is the i th eigen vector of Y and D is a diagonal matrix with its i th diagonal element equal to the i th eigen value, λ i , of Y. Using Eq. ( 10) in Eq. ( 9) we get where Since D in Eq. ( 11) is a diagonal matrix, each row of Eq. ( 11) can be solved independently resulting in x 0 2 e λ 2 z . . .
The column matrix X 0 = V −1 B 0 n is composed of the boundary values x 0 1 to x 0 m .Once X is evaluated using Eq. ( 13), the solution to Eq. ( 9) can be found using the transformation in Eq. (12).That is, B n = VX.Thus, the explicit solution to Eq. ( 9) is given by Without loss of generality, here we have assumed that the eigen values of matrix Y are real and distinct, which is the case for most of the widely used phase functions in practice.However, if for a particular application the eigen values and eigen vectors contain complex values (in the form of complex conjugate pairs), a real solution to Eq. ( 9) can be obtained as illustrated in [20].For example, suppose λ i = p + qi and λ j = p − qi are two complex eigen values and v i = a + bi and v j = a − bi are their corresponding complex eigen vectors.It can be shown that we can then obtain two real valued solutions to Eq. ( 9), w i = e pz (a cos qz − b sin qz) and w j = e pz (b cos qz + a sin qz) [20].The general solution to Eq. ( 9) then becomes Equation ( 14) can be written in matrix notation as where V = [e λ 1 z v 1 e λ 2 z v 2 . . .w i w j . . .e λ m z v m ] and C is a column matrix with its i th element equal to c i .The constants in C can be found using Eq. ( 15) (ie.C = inv(V)B 0 n ).For the case of repeated eigen values, the general solution to Eq. ( 9) will be of the form [20] where the eigen value λ i is repeated 3 times and the corresponding generalized eigenvectors, v i , v i+1 and v i+2 are found such that (Y − λ i I) [20].

Results and discussion
Figure 1 shows a comparison of irradiance profiles obtained from the work proposed in this paper (Eigen method) and the standard solution (obtained using the LRKF method [14]), at the exit plane of a layer of a human skin specimen of 1 mm thickness, illuminated by a Gaussian pulse of I 0 intensity (of arbitrary units).The Henyey-Greenstein phase function with an asymmetry factor of 0.7 together with σ a = 0.5 mm −1 , σ s = 0.8 mm −1 ,v = 0.215 mm/ps and T = 1.5 was used for this simulation.The results produced by the two methods for the same number of discrete ordinates were exactly the same as shown by Fig. 1.In the proposed algorithm numerical errors occur due to the truncation of the Laguerre series and due to the finite number of discrete ordinates chosen for θ and φ .However, as discussed  [21] and [14], it is possible to obtain a very accurate Laguerre fit, in the observation window, with only 63 terms for a Gaussian pulse with an arbitrary width by using a suitable scaling factor.Hence, we can consider that the numerical errors are introduced only due to the discretization of the zenith and azimuthal angles.Figure 2 shows a comparison of execution times of the Eigen method proposed in this paper with the standard solution.The rate of increase of the execution time with the number of discrete ordinates is significantly low for the Eigen method compared to the standard solution.

Conclusion
From the simulation results it can be concluded that the proposed algorithm is accurate to the same degree compared to the standard solution which is based on a numerical solution to the reduced vectorized version of the PTE.However, the proposed algorithm is significantly faster in execution time compared to the standard solution method.For example, increasing the number of directions from 25 to 196 increased the execution time of the proposed algorithm by only about 11 seconds (6.4%)where as that of the standard method was increased by about 166 seconds (97.0%).