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

Optics Express, Vol. 19, Issue 4, pp. 2922-2927 (2011)

http://dx.doi.org/10.1364/OE.19.002922

Acrobat PDF (649 KB)

### Abstract

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 integro-differential 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

## 1. Introduction

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]

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]

*F*method [6

_{N}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]

*et al.*[8

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

*et al.*[9

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]

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]

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

*et al.*[14

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]

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]

## 2. Eigen decomposition method

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]

*I*(

*z*,

*u*,

*ϕ*,

*t*) is the light intensity (radiance), (

*z*,

*θ*,

*ϕ*) are the standard spherical coordinates,

*u*= cos

*θ*,

*t*denotes time,

*σ*and

_{t}*σ*are attenuation and scattering coefficients, respectively. The speed of light in the medium is denoted by

_{s}*v*and

*P*(

*u*′,

*ϕ*′;

*u*,

*ϕ*) is the phase function.

*I*(

*z*= 0,

*u*

_{0},

*ϕ*

_{0},

*t*) =

*f*(

*t*)

*δ*(

*u*−

*u*

_{0})

*δ*(

*ϕ*−

*ϕ*

_{0}), 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 subject to

*I*(

*z*= 0,

*u*,

*ϕ*,

*τ*) =

*f*(

*τ*)

*δ*(

*u*−

*u*

_{0})

*δ*(

*ϕ*−

*ϕ*

_{0}). Use of Eq. (2) in Eq. (1) eliminates the time derivative term in Eq. (1).

*ϕ*, can be discretized using the discrete ordinates method [4

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]

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]

*ϕ*.

*W*(

*x*) is the weight function and

*f*(

*x*) is any arbitrary function, by a summation such that

*r*= 1, . . . ,

*L*,

*ϕ*is the

_{j}*j*quadrature point of

^{th}*ϕ*and

*u*, of Eq. (4) using the discrete ordinates method, which results in where

*i*= 1, . . . ,

*K*,

*r*= 1, . . . ,

*L*,

*u*is the

_{k}*k*quadrature point of

^{th}*u*and

*τ*, we can expand

*I*(

*z*,

*u*,

_{i}*ϕ*,

_{r}*τ*) and

*f*(

*τ*) using Laguerre polynomials with respect to

*τ*, such that

*B*and

_{k}*F*are the coefficients corresponding to the Laguerre polynomial

_{k}*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*and Eq. (6) corresponds to the

_{n}*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*repeating

_{K}*L*times. The matrix

**W**is also a (

*K*×

*L*) by (

*K*×

*L*) diagonal matrix with diagonal elements

**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

**Y**can be written as a product of three matrices composed of its eigen values and eigen vectors, as follows. where

**V**= [

**v**

_{1}

**v**

_{2}. . .

**v**

*].*

_{m}**v**

*is the*

_{i}*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,

*λ*, of

_{i}**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 The column matrix

**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**+

**b**

*i*and

**v**

*=*

_{j}**a**−

**b**

*i*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*

^{λ1z}

**v**

_{1}

*e*

^{λ2z}

**v**

_{2}. . .

**w**

_{i}**w**

*. . .*

_{j}*e*

^{λmz}

**v**

*] and*

_{m}**C**is a column matrix with its

*i*element equal to

^{th}*c*. The constants in

_{i}**C**can be found using Eq. (15) (ie.

*λ*is repeated 3 times and the corresponding generalized eigenvectors,

_{i}**v**

*,*

_{i}**v**

_{i}_{+1}and

**v**

_{i}_{+2}are found such that (

**Y**−

*λ*

_{i}**I**)

^{2}

**v**

_{i}_{+2}=

**0**, (

**Y**−

*λ*

_{i}**I**)

**v**

_{i}_{+2}=

**v**

_{i+1}, (

**Y**−

*λ*

_{i}**I**)

**v**

_{i}_{+1}=

**v**

*and (*

_{i}**Y**−

*λ*

_{i}**I**)

**v**

*=*

_{i}**0**as illustrated in [20].

## 3. Results and discussion

**14**, 105–112 (2008). [CrossRef]

*I*

_{0}intensity (of arbitrary units). The Henyey-Greenstein phase function with an asymmetry factor of 0.7 together with

*σ*= 0.5 mm

_{a}^{−1},

*σ*= 0.8 mm

_{s}^{−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.

*θ*and

*ϕ*. However, as discussed in detail in [21

21. C. C. Handapangoda and M. Premaratne, “Implicitly causality enforced solution of multidimensional transient photon transport equation,” Opt. Express **17**, 23423–23442 (2009). [CrossRef]

**14**, 105–112 (2008). [CrossRef]

## 4. Conclusion

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

2. | N. Y. Gnedin, “Multi-dimensional cosmological radiative transfer with a variable Eddington tensor formalism,” New Astron. |

3. | D. M. O’Brien, “Accelerated quasi Monte Carlo integration of the radiative transfer equation,” J. Quant. Spectrosc. Radiat. Transf. |

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

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

6. | C. E. Siewert and J. R. Thomas, “Radiative transfer calculations in spheres and cylinders,” J. Quant. Spectrosc. Radiat. Transf. |

7. | C. E. Siewert, “A radiative-transfer inverse-source problem for a sphere,” J. Quant. Spectrosc. Radiat. Transf. |

8. | E. W. Larsen, “The inverse source problem in radiative transfer,” J. Quant. Spectrosc. Radiat. Transf. |

9. | G. C. Pomraning and G. M. Foglesong, “Transport-diffusion interfaces in radiative transfer,” J. Comput. Phys. |

10. | Z. M. Tan and P. F. Hsu, “An integral formulation of transient radiative transfer,” ASME J. Heat Transfer |

11. | J. A. Fleck, “The calculation of nonlinear radiation transport by a Monte Carlo method,” Methods Comput. Phys. |

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

13. | A. D. Kim and M. Moscoso, “Chebyshev spectral methods for radiative transfer,” SIAM J. Sci. Comput. (USA) |

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

15. | A. D. Kim and M. Moscoso, “Chebyshev spectral methods fro radiative transfer,” SIAM J. Sci. Comput. (USA) |

16. | A. B. Carlson, |

17. | S. Chandrasekhar, |

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

19. | W. H. Press, S. A. Teukolsk, W. T. Vetterling, and B. P. Flannery, |

20. | C. H. Edwards and D. E. Penney, |

21. | C. C. Handapangoda and M. Premaratne, “Implicitly causality enforced solution of multidimensional transient photon transport equation,” Opt. Express |

**OCIS Codes**

(080.2720) Geometric optics : Mathematical methods (general)

(290.7050) Scattering : Turbid media

(010.5620) Atmospheric and oceanic optics : Radiative transfer

**ToC Category:**

Scattering

**History**

Original Manuscript: October 25, 2010

Revised Manuscript: January 27, 2011

Manuscript Accepted: January 28, 2011

Published: February 1, 2011

**Virtual Issues**

Vol. 6, Iss. 3 *Virtual Journal for Biomedical Optics*

**Citation**

Chintha C. Handapangoda, Pubudu N. Pathirana, and Malin Premaratne, "Eigen decomposition solution to the one-dimensional time-dependent photon transport equation," Opt. Express **19**, 2922-2927 (2011)

http://www.opticsinfobase.org/oe/abstract.cfm?URI=oe-19-4-2922

Sort: Year | Journal | Reset

### References

- 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]
- N. Y. Gnedin, “Multi-dimensional cosmological radiative transfer with a variable Eddington tensor formalism,” N. Astron. 6, 437–455 (2001). [CrossRef]
- D. M. O’Brien, “Accelerated quasi Monte Carlo integration of the radiative transfer equation,” J. Quant. Spectrosc. Radiat. Transf. 48, 41–59 (1992). [CrossRef]
- 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). [CrossRef] [PubMed]
- 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]
- C. E. Siewert and J. R. Thomas, “Radiative transfer calculations in spheres and cylinders,” J. Quant. Spectrosc. Radiat. Transf. 34, 59–64 (1985). [CrossRef]
- C. E. Siewert, “A radiative-transfer inverse-source problem for a sphere,” J. Quant. Spectrosc. Radiat. Transf. 52, 157–160 (1994). [CrossRef]
- E. W. Larsen, “The inverse source problem in radiative transfer,” J. Quant. Spectrosc. Radiat. Transf. 15, 1–5 (1975). [CrossRef]
- G. C. Pomraning and G. M. Foglesong, “Transport-diffusion interfaces in radiative transfer,” J. Comput. Phys. 32, 420–436 (1979). [CrossRef]
- Z. M. Tan and P. F. Hsu, “An integral formulation of transient radiative transfer,” ASME J. Heat Transfer 123, 466–475 (2001). [CrossRef]
- J. A. Fleck, “The calculation of nonlinear radiation transport by a Monte Carlo method,” Methods Comput. Phys. 1, 43–65 (1963).
- 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]
- A. D. Kim and M. Moscoso, “Chebyshev spectral methods for radiative transfer,” SIAM J. Sci. Comput. (USA) 23, 2074–2094 (2002). [CrossRef]
- 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]
- A. D. Kim and M. Moscoso, “Chebyshev spectral methods fro radiative transfer,” SIAM J. Sci. Comput. (USA) 23, 2074–2094 (2002). [CrossRef]
- A. B. Carlson, Communication Systems: An Introduction to Signals and Noise in Electrical Communication (McGraw-Hill, 1986).
- S. Chandrasekhar, Radiative Transfer (Dover Publications, 1960).
- 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]
- W. H. Press, S. A. Teukolsk, W. T. Vetterling, and B. P. Flannery, Numerical Recipes in C + +, 2nd ed. (Cambridge University Press, 2003).
- C. H. Edwards and D. E. Penney, Differential Equations: Computing andModeling, 3rd ed. (Prentice Hall, 2004).
- C. C. Handapangoda and M. Premaratne, “Implicitly causality enforced solution of multidimensional transient photon transport equation,” Opt. Express 17, 23423–23442 (2009). [CrossRef]

## Cited By |
Alert me when this paper is cited |

OSA is able to provide readers links to articles that cite this paper by participating in CrossRef's Cited-By Linking service. CrossRef includes content from more than 3000 publishers and societies. In addition to listing OSA journal articles that cite this paper, citing articles from other participating publishers will also be listed.

« Previous Article | Next Article »

OSA is a member of CrossRef.