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

Chintha C. Handapangoda; Pubudu N. Pathirana; Malin Premaratne

doi:10.1364/OE.19.002922

1. 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.

2. Eigen decomposition method

The one dimensional transient photon transport equation for a medium without any internal sources is given by [14, 15]

\begin{array}{r} \frac{1}{v} \frac{\partial}{\partial t} I (z, u, ϕ, t) + u \frac{\partial}{\partial z} I (z, u, ϕ, t) - \frac{σ_{s}}{4 π} \int_{0}^{2 π} \int_{- 1}^{1} P (u^{'}, ϕ^{'}; u, ϕ) I (z, u^{'}, ϕ^{'}, t) d u^{'} d ϕ^{'} \\ + σ_{t} I (z, u, ϕ, t) = 0, \end{array}

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 I (z = 0, u₀, ϕ₀, t) = f (t)δ(u − u₀)δ(ϕ − ϕ₀), where δ(x) is the Dirac delta function [16] and f (t) is the temporal profile of the incident pulse. Using the transformation

τ = t - \frac{z}{v u}

which maps Eq. (1) to a moving reference frame with the pulse, we get

u \frac{\partial}{\partial z} I (z, u, ϕ, τ) + σ_{t} I (z, u, ϕ, τ) - \frac{σ_{s}}{4 π} \int_{0}^{2 π} \int_{- 1}^{1} P (u^{'}, ϕ^{'}; u, ϕ) I (z, u^{'}, ϕ^{'}, τ) d u^{'} d ϕ^{'} = 0,

subject to I (z = 0, u, ϕ,τ) = f (τ)δ(u − u₀)δ(ϕ − ϕ₀). Use of Eq. (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 $\int_{a}^{b} W (x) f (x) d x$ , where W (x) is the weight function and f (x) is any arbitrary function, by a summation such that $\int_{a}^{b} W (x) f (x) d x \approx \sum_{j = 0}^{N - 1} 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

u \frac{\partial}{\partial z} I (z, u, ϕ_{r}, τ) + σ_{t} I (z, u, ϕ_{r}, τ) - \frac{σ_{s}}{4 π} \sum_{j = 1}^{L} w_{j}^{ϕ} \int_{- 1}^{1} P (u^{'}, ϕ_{j}; u, ϕ_{r}) I (z, u^{'}, ϕ_{j}, τ) d u^{'} = 0,

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

u_{i} \frac{\partial}{\partial z} I (z, u_{i}, ϕ_{r}, τ) + σ_{t} I (z, u_{i}, ϕ_{r}, τ) - \frac{σ_{s}}{4 π} \sum_{j = 1}^{L} w_{j}^{ϕ} \sum_{k = 1}^{K} w_{k}^{u} P (u_{k}, ϕ_{j}; u_{i}, ϕ_{r}) I (z, u_{k}, ϕ_{j}, τ) = 0,

where i = 1, . . . , K, r = 1, . . . , L, u_k is the k^th quadrature point of u and

w_{k}^{u}

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 $I (z, u_{i}, ϕ_{r}, τ) \approx \sum_{k = 0}^{N} B_{k} (z, u_{i}, ϕ_{r}) L_{k} (τ)$ and $f (τ) \approx \sum_{k = 0}^{N} F_{k} L_{k} (τ)$ where B_k and 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,

u_{i} \frac{\partial}{\partial z} B_{n} (z, u_{i}, ϕ_{r}) + σ_{t} B_{n} (z, u_{i}, ϕ_{r}) - \frac{σ_{s}}{4 π} \sum_{j = 1}^{L} w_{j}^{ϕ} \sum_{k = 1}^{K} w_{k}^{u} P (u_{k}, ϕ_{j}; u_{i}, ϕ_{r}) B_{n} (z, u_{k}, ϕ_{j}) = 0,

subject to

B_{n} (z = 0, u_{i}, ϕ_{r}) = F_{n} δ (u_{i} - u_{0}) δ (ϕ_{r} - ϕ_{0}),

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;

\frac{d}{d z} {AB}_{n} + σ_{t} B_{n} - \frac{σ_{s}}{4 π} PW B_{n} = 0,

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₁ to u_K repeating L times. The matrix W is also a (K × L) by (K × L) diagonal matrix with diagonal elements

w_{r}^{ϕ} \times w_{k}^{u}

with the pattern

w_{1}^{ϕ} \times w_{1}^{u}

,

w_{1}^{ϕ} \times w_{2}^{u}

, . . . ,

w_{1}^{ϕ} \times w_{K}^{u}

,

w_{2}^{ϕ} \times w_{1}^{u}

, . . . ,

w_{L}^{ϕ} \times w_{K}^{u}

.

Rearranging, Eq. (8) can be written as

\frac{d}{d z} B_{n} = {YB}_{n},

where

Y = A^{- 1} [\frac{σ_{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_{n}^{0} = {[B_{n} (z = 0, u_{i}, ϕ_{r})]}_{K \times 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.

Y = {VDV}^{- 1},

where V = [v₁ v₂ . . . v_m]. 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

\begin{array}{l} \frac{d}{d z} B_{n} & = & Y B_{n} \\ = & {VDV}^{- 1} B_{n} \\ V^{- 1} \frac{d}{d z} B_{n} & = & D V^{- 1} B_{n} \\ \frac{d}{d z} X & = & DX \end{array}

where

X = V^{- 1} B_{n} .

Since D in Eq. (11) is a diagonal matrix, each row of Eq. (11) can be solved independently resulting in

X = (\begin{matrix} x_{1} \\ x_{2} \\ ⋮ \\ x_{m} \end{matrix}) = (\begin{matrix} x_{1}^{0} e^{λ_{1} z} \\ x_{2}^{0} e^{λ_{2} z} \\ ⋮ \\ x_{m}^{0} e^{λ_{m} z} \end{matrix}) .

The column matrix

X^{0} = V^{- 1} B_{n}^{0}

is composed of the boundary values

x_{1}^{0}

to

x_{m}^{0}

. 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

B_{n} = x_{1}^{0} e^{λ_{1} z} v_{1} + x_{2}^{0} e^{λ_{2} z} v_{2} + \dots + x_{m}^{0} e^{λ_{m} z} v_{m} .

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

B_{n} = c_{1} e^{λ_{1} z} v_{1} + c_{2} e^{λ_{2} z} v_{2} + \dots + c_{i} w_{i} + c_{j} w_{j} + \dots + c_{m} e^{λ_{m} z} v_{m} .

Equation (14) can be written in matrix notation as

B_{n} = VC,

where V = [e^λ₁zv₁ e^λ₂zv₂ . . . w_i w_j . . . e^λ_mzv_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_{n}^{0}

).

For the case of repeated eigen values, the general solution to Eq. (9) will be of the form [20]

\begin{array}{r} B_{n} = c_{1} e^{λ_{1} z} v_{1} + c_{2} e^{λ_{2} z} v_{2} + \dots + c_{i} e^{λ_{i} z} v_{i} + c_{i + 1} e^{λ_{i} z} (v_{i} z + v_{i + 1}) \\ + c_{i + 2} e^{λ_{i} z} (\frac{1}{2} v_{i} z^{2} + v_{i + 1} z + v_{i + 2}) + \dots + c_{m} e^{λ_{m} z} v_{m}, \end{array}

where the eigen value λ_i is repeated 3 times and the corresponding generalized eigenvectors, v_i, v_i₊₁ and v_i₊₂ are found such that (Y − λ_iI)² v_i₊₂ = 0, (Y − λ_iI) v_i₊₂ = v_i+1, (Y − λ_iI) v_i₊₁ = v_i and (Y − λ_iI) v_i = 0 as illustrated in [20].

3. 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₀ intensity (of arbitrary units). The Henyey-Greenstein phase function with an asymmetry factor of 0.7 together with σ_a = 0.5 mm⁻¹, σ_s = 0.8 mm⁻¹, 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.

Fig. 1 Irradiance profile transmitted from a human skin specimen

Download Full Size | PDF

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 in detail in [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.

Fig. 2 A comparision of execution times for different numbers of directions

Download Full Size | PDF

4. 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%).

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). [CrossRef]

2. N. Y. Gnedin, “Multi-dimensional cosmological radiative transfer with a variable Eddington tensor formalism,” New Astron. 6, 437–455 (2001). [CrossRef]

3. D. M. O’Brien, “Accelerated quasi Monte Carlo integration of the radiative transfer equation,” J. Quant. Spectrosc. Radiat. Transf. 48, 41–59 (1992). [CrossRef]

4. K. Stamnes, S. C. Tsay, W. Wiscombe, and K. Jayaweera, “Numerically stable algorithm for discrete-ordinate-method radiative transfer in multiple scattering and emitting layered media,” Appl. Opt. 27, 2502–2509 (1988). [CrossRef] [PubMed]

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). [CrossRef]

6. C. E. Siewert and J. R. Thomas, “Radiative transfer calculations in spheres and cylinders,” J. Quant. Spectrosc. Radiat. Transf. 34, 59–64 (1985). [CrossRef]

7. C. E. Siewert, “A radiative-transfer inverse-source problem for a sphere,” J. Quant. Spectrosc. Radiat. Transf. 52, 157–160 (1994). [CrossRef]

8. E. W. Larsen, “The inverse source problem in radiative transfer,” J. Quant. Spectrosc. Radiat. Transf. 15, 1–5 (1975). [CrossRef]

9. G. C. Pomraning and G. M. Foglesong, “Transport-diffusion interfaces in radiative transfer,” J. Comput. Phys. 32, 420–436 (1979). [CrossRef]

10. Z. M. Tan and P. F. Hsu, “An integral formulation of transient radiative transfer,” ASME J. Heat Transfer 123, 466–475 (2001). [CrossRef]

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). [CrossRef]

13. A. D. Kim and M. Moscoso, “Chebyshev spectral methods for radiative transfer,” SIAM J. Sci. Comput. (USA) 23, 2074–2094 (2002). [CrossRef]

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). [CrossRef]

15. A. D. Kim and M. Moscoso, “Chebyshev spectral methods fro radiative transfer,” SIAM J. Sci. Comput. (USA) 23, 2074–2094 (2002). [CrossRef]

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). [CrossRef]

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). [CrossRef]

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

Abstract

1. Introduction

2. Eigen decomposition method

3. Results and discussion

4. Conclusion

References and links

Cited By

Figures (2)

Equations (17)

Optics Express