Infinite space Green’s function of the time-dependent radiative transfer equation

André Liemert; Alwin Kienle

doi:10.1364/BOE.3.000543

1. Introduction

The radiative transfer equation (RTE), an integro-partial differential equation, is widely used in many areas of physics to model the propagation of waves in random scattering media. Examples are the light propagation in biological tissue [1], the neutron transport theory [2, 3], the electromagnetic waves propagation in plasmas or in the atmosphere. Although this equation is used since many decades complete analytical solutions for the steady-state and time domains are only available for the simplest cases [2–4]. Due to this fact the RTE is usually solved by numerical approaches such as the Monte Carlo method [5], the discrete ordinate method [6], the finite-element method [7] or is approximated by the diffusion equation [8]. In general, numerical approaches are in contrast to analytical methods computationally very demanding.

Recently, elegant analytical approaches for solving the RTE in the steady-state domain were discussed in [9–11]. In principle, it is well-known that solutions obtained in the steady-state domain can be transformed into the frequency-domain resulting in a complex-valued absorption coefficient. Then, by performing numerically an inverse Fourier transform of the solution in the frequency domain the solution of the time-dependent RTE is obtained. However, this procedure involves in many cases difficulties due to the fact that the Green’s function of the steady-state RTE is already not a square integrable function [10]. Therefore the resulting time domain signal becomes a strongly oscillating curve which is inappropriate for an exact analysis or applications within inverse problems. Furthermore, the performance of the fast Fourier transform, which indeed would lead to decreasing calculation times, is circuitous if specified time values are required. In the publication of [12] the derivation of a cumulant solution for the time-dependent RTE was presented.

In this article an infinite space Green’s function of the time-dependent RTE is derived where the final solutions have an analytical dependence on the time variable and are well convergent due to the exponential decay of the time-dependent modes. Within the derivation no approximations are made so that also the causality condition is satisfied; see the Results section. The solutions are successfully verified by comparisons with Monte Carlo simulations [1]. The derived Green’s function cannot only be used for applications in infinitely extended scattering media but for all applications where the well-known approximated solution obtained from the diffusion equation is used [13], as was demonstrated for the case of isotropic scattering [14].

2. Theory

The time-dependent radiative transfer equation for the radiance I = I(r,μ,t) caused by the internal source density Q = Q(r,t) in a spherical symmetric medium has the form

\frac{1}{c} \frac{\partial I}{\partial t} + μ \frac{\partial I}{\partial r} + \frac{1 - μ^{2}}{r} \frac{\partial I}{\partial μ} + σ_{t} I = σ_{s} A I + \frac{Q}{4 π},

where A is the scattering operator defined by

A I (r, μ, t) = \int f (\hat{s} \cdot {\hat{s}}^{'}) I (r, μ^{'}, t) d^{2} {\hat{s}}^{'} .

Here σ_t = σ_a + σ_s is the total interaction coefficient, σ_a is the absorption coefficient, σ_s is the scattering coefficient and μ = ŝ·r̂ is the cosine of the angle between the direction of propagation ŝ and the unit vector of the position vector r̂. The function f (ŝ · ŝ′) describes a rotationally invariant phase function which can be expanded in Legendre polynomials P_l(x), where the expansion coefficients f_l are defined by

f_{l} = 2 π \int_{- 1}^{1} f (μ) P_{l} (μ) d μ .

For finding the impulse response in space and time the source term becomes to

Q (r, t) = \frac{δ (r)}{4 π r^{2}} δ (t) .

In the following a plane wave decomposition of a rotationally invariant function in spherical polar coordinates (r,θ,φ) is derived. The radiance can be written as a three-dimensional Fourier integral

I (r, cos θ_{r}, t) = \frac{1}{{(2 π)}^{3}} \int I (k, cos θ_{k}, t) exp (j k \cdot r) d^{3} k .

Note that the indices r and k of the angular variables θ_k and φ_k indicate the spatial and transformed space, respectively. By using the relation [15]

e^{jkr cos (θ_{r} - θ_{k})} = 4 π \sum_{l = 0}^{\infty} j^{l} j_{l} (k r) \sum_{m = - l}^{l} Y_{l m}^{*} (θ_{k}, φ_{k}) Y_{l m} (θ_{r}, φ_{r}),

where Y_lm(θ,φ) are the spherical harmonics, and defining expansion coefficients

I_{l} (k, t) = 2 π \int_{0}^{π} I (k, cos θ_{k}, t) P_{l} (cos θ_{k}) sin θ_{k} d θ_{k}

yields the Legendre polynomial series

I (r, μ, t) = \sum_{l = 0}^{\infty} \frac{2 l + 1}{4 π} I_{l} (r, t) P_{l} (μ),

where the radial dependent kernel has the form of an inverse spherical Hankel transform

I_{l} (r, t) = \frac{j^{l}}{2 π^{2}} \int_{0}^{\infty} I_{l} (k, t) j_{l} (k r) k^{2} d k .

The belonging forward transform is found via orthogonality as

I_{l} (k, t) = 4 π {(- j)}^{l} \int_{0}^{\infty} I_{l} (r, t) j_{l} (k r) r^{2} d r .

The function j_l(x) denotes the spherical Bessel function of the first kind. Inserting the series representation with the belonging radial dependent coefficients leads to a linear and time invariant system (LTI system)

\frac{1}{c} \frac{d I_{l} (k, t)}{dt} + \frac{l}{2 l + 1} j k I_{l - 1} (k, t) + σ_{l} I_{l} (k, t) + \frac{l + 1}{2 l + 1} j k I_{l + 1} (k, t) = δ (t) q_{l},

where σ_l = σ_a + (1 – f_l)σ_s and q_l = δ_l₀. For the Henyey-Greenstein model [8] the expansion coefficients become f_l = g^l ∀ l ∈ 𝕅₀. In P_N approximation (N odd) equations must be considered only for l = 0,...,N with the initial conditions I₋₁(k,t) = I_N₊₁(k,t) = 0. The time-dependent impulse response of a LTI system can be obtained via the one-sided Laplace transform with zero initial conditions

I_{l} (k, s) = 𝒧 {I_{l} (k, t)} = \int_{0}^{\infty} I_{l} (k, t) e^{- s t} d t .

Due to the above integral transforms the LTI system (11) is converted to a set of linear equations. After a few simple matrix conversions the system of linear equations becomes

Q {[I_{0} (k, s), I_{1} (k, s), \dots, I_{N} (k, s)]}^{T} = q,

where q = [1,0,...,0]^T and

Q = (A + s / c I) \cdot diag (\sqrt{1}, \sqrt{3}, \dots, \sqrt{2 N + 1})

. Here the matrix I denotes the identity matrix. The complex symmetric tridiagonal matrix A = A^T has the form

A = (\begin{matrix} σ_{0} & \frac{j k}{\sqrt{1 \cdot 3}} & 0 & 0 & \dots & 0 \\ \frac{j k}{\sqrt{1 \cdot 3}} & σ_{1} & \frac{j 2 k}{\sqrt{3 \cdot 5}} & 0 & \dots & ⋮ \\ 0 & \frac{j 2 k}{\sqrt{3 \cdot 5}} & ⋱ & ⋱ & \dots & 0 \\ 0 & 0 & ⋱ & ⋱ & ⋱ & 0 \\ ⋮ & \dots & \dots & ⋱ & ⋱ & \frac{j N k}{\sqrt{4 N^{2} - 1}} \\ 0 & \dots & 0 & 0 & \frac{j N k}{\sqrt{4 N^{2} - 1}} & σ_{N} \end{matrix}) .

The solution of the matrix equation Eq. (13) is formally obtained as [16]

{[I_{0} (k, s), I_{1} (k, s), \dots, I_{N} (k, s)]}^{T} = D {[A + \frac{s}{c} I]}^{- 1} q .

Here

D = diag (\frac{1}{\sqrt{1}}, \frac{1}{\sqrt{3}}, \dots, \frac{1}{\sqrt{2 N + 1}})

. The inverse Laplace transform can be performed using a well-known transform pair within the linear system theory for the so-called transition matrix

exp (- A t) = 𝒧^{- 1} {{(s I + A)}^{- 1}},

Thus the time-dependent vector components become for t > 0

{[I_{0} (k, t), I_{1} (k, t), \dots, I_{N} (k, t)]}^{T} = c D exp (- c A t) q .

A simple but inefficient method for evaluation of the matrix exponential is given by the power series

exp (- c A t) = \sum_{n = 0}^{\infty} \frac{{(- c t)}^{n}}{n!} A^{n} = I - \frac{c t}{1!} A + \frac{{(c t)}^{2}}{2!} A^{2} - \frac{{(c t)}^{3}}{3!} A^{3} \pm \dots,

which converge for an arbitrary square complex-valued matrix. Despite the fact that mathematical software generally provide the matrix exponential it is maybe advantageous with respect to the computation time to avoid formally this function and performing the eigenvalue decomposition A=UΛU⁻¹, where U is the transformation matrix, which contains the eigenvectors |u_n〉 (n = 0,1,...,N) of A. The diagonal matrix Λ =diag(λ₀,λ₁,...,λ_N) (λ_n ∈ 𝔺) contains the eigenvalues of the equation det(A −λI) = 0. Due to the above eigenvalue decomposition equation Eq. (17) simplifies to

{[I_{0} (k, t), I_{1} (k, t), \dots, I_{N} (k, t)]}^{T} = c D U exp (- c Λ t) U^{- 1} q,

where the matrix exponential can now be evaluated as

exp (- c Λ t) = (\begin{matrix} exp (- c λ_{0} t) & 0 & \dots & 0 \\ 0 & exp (- c λ_{1} t) & \dots & 0 \\ ⋮ & ⋮ & ⋱ & ⋮ \\ 0 & 0 & \dots & exp (- c λ_{N} t) \end{matrix}) .

Note that all expansion coefficients of the Legendre series Eq. (8) have a time decaying behaviour which arises from the sign of the real part of the complex eigenvalues. After some algebraic conversions the time-dependent modes needed for the radiance Eq. (8) are obtained as

I_{l} (k, t) = \frac{c}{\sqrt{2 l + 1}} \sum_{n = 0}^{N} 〈 n | {\tilde{u}}_{0} 〉 〈 u_{n} | l 〉 exp (- c λ_{n} t),

where the vector |ũ₀〉 denotes the first column vector of the matrix U⁻¹. Note that for obtaining an expression for the radiance caused by an isotropic emitting spherical shell source of the form

Q (r, t) = \frac{δ (r - r^{'})}{4 π r^{2}} δ (t),

where r′ denotes the radius of the spherical shell, the source coefficients q_l = δ_l₀ in system (11) must be replaced by q_l = (− j)^l j_l(kr′). For a source distribution with an arbitrary spatial dependence the resulting radiance can be obtained via convolution.

Furthermore, the result for the fluence caused by an unidirectional source which illuminates in direction ŝ₀ is the same as for the obtained radiance caused by a isotropic source multiplied by the factor 4π. Then the above defined cosine of the angle between the direction of propagation and the vector of the position has now the meaning μ = ŝ₀ · r̂.

2.1. Numerical implementation

For a fast and accurate evaluation of the obtained solution we replace the continuous Hankel transform over an infinite range, equation Eq. (10), by the corresponding finite version

I_{k l} (t) = 4 π {(- j)}^{l} \int_{0}^{R} I_{l} (r, t) j_{l} (ξ_{l k} r) r^{2} d r,

where R equals the radius of a large sphere. Now the continuous wave number k in relation Eq. (10) becomes to discrete values ξ_kl which are dependent on the transform order l and are given by the positive roots of the equation j_l(ξ_klR) = 0. Inverting the finite transform via the completeness relation of the spherical Bessel functions yields an alternative inverse transform for the integral representation Eq. (9)

I_{l} (r, t) = \frac{j^{l}}{2 π R^{3}} \sum_{ξ_{l k} > 0} I_{k l} (t) \frac{j_{l} (ξ_{l k} r)}{j_{l + 1}^{2} (ξ_{l k} R)} .

Due to the exponential decay of the time-dependent modes I_kl(t) the solution, in general, converges already with a small number of discrete wave numbers. Additionally, we performed comparisons between the finite Hankel transform and the continuous version over an infinite range (R → ∞) which was evaluated by a Gaussian quadrature resulting in an excellent agreement. Note that especially for the evaluation of the fluence the required roots are easily given by the positive zeros of the sine function ξ_k₀ = kπ/R ∀ k ∈ 𝕅. For every root the system of linear equations must be solved via eigenvalue decomposition of the matrix A. The whole procedure can be easily implemented within a small MATLAB script, see e.g. on the internet [17]. In contrast to the steady-state theory here the eigenvalues and eigenvectors become additionally complex-valued. However, if an eigenvalue λ_n becomes complex then

λ_{n}^{*}

is also an eigenvalue due to the symmetry to the real axis in the complex plane so that in all cases the functions I_l(r,t) are automatically real-valued. The fluence [5] which is given by the lowest order coefficient I₀(r,t) becomes the simple form

Φ (r, t) = \int I (r, μ^{'}, t) d^{2} {\hat{s}}^{'} = \frac{1}{2 r R^{2}} \sum_{k = 1}^{\infty} k I_{k 0} (t) sin (ξ_{k 0} r) .

Regarding the computation time the fluence in P₂₅ approximation using 250 discrete wave numbers can be evaluated for 500 time values in ≈ 1 s using a state of the art personal computer and a MATLAB scrip. Note that such an approximation order is much more than needed in most relevant cases. It is mentionable that the fluence in P₃ approximation, which is a real alternative to the often used diffusion Green’s function, can be evaluated completely analytically without a numerical eigenvalue decomposition.

3. Results

In this section the obtained solutions are validated against Monte Carlo simulations and partly the diffusion theory. Monte Carlo simulations are in the limit of an infinitely large number of simulated photons an exact solution of the RTE. The corresponding equations for the infinite space fluence and radiance within the diffusion approximation can be obtained from [8]. In the following comparisons a isotropic emitting point source distribution and the Henyey-Greenstein phase function are considered. The refractive index of the infinite medium is assumed as n = 1.4. First, we compared the obtained expression for the fluence, equation Eq. (25), with the Monte Carlo method for an anisotropically scattering medium in Fig. 1. It can be seen that the derived solution and the Monte Carlo simulation converge to the same fluence. Due to the fact that the two curves are practically the same, the relative differences between both methods are shown in the inset of Fig. 1. By increasing the number of simulated photons used in the Monte Carlo simulation these differences can be made arbitrarily small. For the next comparison the fluence for an higher absorbing medium and two different anisotropy factors is shown in Fig. 2.

Fig. 1 Time-resolved fluence at r = 3.05mm in an anisotropically scattering medium with properties σ_a = 0.01mm⁻¹, σ_s = 1mm⁻¹ and g = 0.9.

Download Full Size | PDF

Fig. 2 Time-resolved fluence at r = 3.05mm for two different anisotropy factors in a scattering medium with properties σ_a = 0.1mm⁻¹ and σ_s = 1mm⁻¹.

Download Full Size | PDF

The derived fluence (solid curves) agrees again with the Monte Carlo method (symbols) for both anisotropy factors. Figure 3 shows the fluence obtained from the analytical solution and the Monte Carlo method for a relatively small distance to the isotropic emitting point source. Additionally, the fluence within the diffusion approximation is included for two different diffusion coefficients D (D = 1/(3σ₁) and D = 1/(3σ′_s)) reported in literature [8], where σ′_s = (1 − g)μ_s and σ₁ = μ_a + (1 − g)μ_s.

Fig. 3 Time-resolved fluence at r = 1.05mm in an anisotropically scattering medium with properties σ_a = 0.1mm⁻¹, σ_s = 1mm⁻¹ and g = 0.9.

Download Full Size | PDF

The derived analytical Green’s function of the fluence (solid curve) converges also for the relatively small source detector distances to the same result as obtained from the Monte Carlo simulation (symbols) whereas the diffusion approximation (dashed-curves) shows for both diffusion coefficients large differences compared to the radiative transfer calculations.

Finally we give a comparison of the time-resolved radiance for three different directions of propagation. The results are shown in Fig. 4 for the angles θ = 0.5°, θ = 90.5° and θ = 179.5° (indicated in the legend by θ = 0, π/2 and π).

Fig. 4 Time-resolved radiance at r = 2.05mm in an anisotropically scattering medium with properties σ_a = 0.1mm⁻¹, σ_s = 1mm⁻¹ and g = 0.8.

Download Full Size | PDF

Similar as for the fluence both methods result within the stochastic nature of the Monte Carlo simulations in the same radiance.

4. Conclusion

In this article time-dependent Green’s functions of the RTE for the radiance and the fluence were derived and validated with the Monte Carlo method resulting in an excellent agreement. The obtained solution is completely analytical up to the P₃ approximation and enables the consideration of an arbitrary phase function which must be expanded in the basis of Legendre polynomials.

The relative differences between the independent solutions of the time-dependent RTE were only caused by the statistics of the Monte Carlo simulations. The main computational cost of the derived Green’s function is due to the numerical determination of several partial fraction coefficients which results in an eigenvalue decomposition of a complex symmetric tridiagonal matrix.

The applications of the derived equations are manifold. They can be directly used for evaluating measurements inside the investigated scattering media, like biological tissue or tissue phantoms [18, 19]. Further, the derived equations can be applied for obtaining solutions for boundary-value problems, e.g. the semi-infinite geometry, considering approximative boundary conditions [14]. In addition, the obtained solutions are important for validation of numerical methods for solving the RTE, e.g. Monte Carlo simulations. Finally, the derived Green’s functions can be applied for the boundary element method to solve the transport equation for an arbitrary geometry [20].

References and links

1. A. Kienle and R. Hibst, “Light Guiding in Biological Tissue due to Scattering,” Phys. Rev. Lett. 97, 018104 (2006). [CrossRef] [PubMed]

2. K. M. Case and P. F. Zweifel, Linear Transport Theory (Addison-Wesley, 1967).

3. J. J. Duderstadt and W. R. Martin, Transport Theory (John Wiley & Sons, 1979).

4. J. C. J. Paasschens, “Solution of the time-dependent Boltzmann equation,” Phys. Rev. E 56, 1135–1141 (1997). [CrossRef]

5. F. Martelli, M. Bassani, L. Alianelli, L. Zangheri, and G. Zaccanti, “Accuracy of the diffusion equation to describe photon migration through an infinite medium: numerical and experimental investigation,” Phys. Med. Biol. 45, 1359–1373 (2000). [CrossRef] [PubMed]

6. A. D. Klose, U. Netz, and A. H. Hielscher, “Iterative reconstruction scheme for optical tomography based on the equation of radiative transfer,” Med. Phys. 26, 1698–1707 (1999). [CrossRef] [PubMed]

7. P. S. Mohan, T. Tarvainen, M. Schweiger, A. Pulkkinen, and S. R. Arridge, “Variable order spherical harmonic expansion scheme for the radiative transport equation using finite elements,” J. Comput. Phys. 230, 7364–7383 (2011). [CrossRef]

8. F. Martelli, S. Del Bianco, A. Ismaelli, and G. Zaccanti, Light Propagation Through Biological Tissue (SPIE Press, 2010). [CrossRef]

9. V. A. Markel, “Modified spherical harmonics method for solving the radiative transport equation,” Waves Random Complex Media 14, L13–L19 (2004).

10. G. Panasyuk, J. C. Schotland, and V. A. Markel, “Radiative transport equation in rotated reference frames,” J. Phys. A 39, 115–137 (2006). [CrossRef]

11. M. Machida, G. Y. Panasyuk, J. C. Schotland, and V. A. Markel, “The Green’s function for the radiative transport equation in the slab geometry,” J. Phys. A 43, 065402 (2010). [CrossRef]

12. W. Cai, M. Lax, and R. R. Alfano, “Cumulant solution of the elastic Boltzmann transport equation in an infinite uniform medium,” Phys. Rev. E 61, 3871–3876 (2000). [CrossRef]

13. M. S. Patterson, B. Chance, and C. Wilson, “Time resolved reflectance and transmittance for the non-invasive measurement of tissue optical properties,” Appl. Opt. 12, 2331–2336 (1989). [CrossRef]

14. F. Martelli, A. Sassaroli, A. Pifferi, A. Torricelli, L. Spinelli, and G. Zaccanti, “Heuristic Greens function of the time dependent radiative transfer equation for a semi-infinite medium, experimental validation,” Opt. Express 15, 18168–18175 (2007). [CrossRef] [PubMed]

15. N. Baddour, “Operational and convolution properties of three-dimensional Fourier transforms in spherical polar coordinates,” J. Opt. Soc. Am. A 27, 2144–2155 (2010). [CrossRef]

16. F. R. Gantmacher, The Theory of Matrices (AMS Chelsea Publishing, 1959).

17. http://www.uni-ulm.de/ilm/index.php?id=10020200.

18. L. C. L. Chin, A. E. Worthington, W. M. Whelan, and I. A. Vitkin, “Determination of the optical properties of turbid media using relative interstitial radiance measurements: Monte Carlo study, experimental validation,” J. Biomed. Opt. 12, 064027 (2007). [CrossRef]

19. L. C. L. Chin, B. Lloyd, W. M. Whelan, and I. A. Vitkin, “Interstitial point radiance spectroscopy of turbid media,” J. Appl. Phys. 105, 102025 (2009). [CrossRef]

20. S. Srinivasan, B. W. Pogue, C. Carpenter, P. K. Yalavarthy, and K. Paulsen, “A boundary element approach for image-guided near-infrared absorption and scatter estimation,” Med. Phys. 34, 4545–4557 (2007). [CrossRef] [PubMed]

Infinite space Green’s function of the time-dependent radiative transfer equation

Abstract

1. Introduction

2. Theory

2.1. Numerical implementation

3. Results

4. Conclusion

References and links

Cited By

Figures (4)

Equations (25)

Biomedical Optics Express