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

This study contains the derivation of an infinite space Green's function of the time-dependent radiative transfer equation in an anisotropically scattering medium based on analytical approaches. The final solutions are analytical regarding the time variable and given by a superposition of real and complex exponential functions. The obtained expressions were successfully validated with Monte Carlo simulations.

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 where A is the scattering operator defined by Here σ = σ + σ is the total interaction coefficient, σ is the absorption coefficient, σ 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(x), where the expansion coefficients f are defined by For finding the impulse response in space and time the source term becomes to 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 Note that the indices r and k of the angular variables θ and φ indicate the spatial and transformed space, respectively. By using the relation [15] where Y(θ,φ) are the spherical harmonics, and defining expansion coefficients yields the Legendre polynomial series where the radial dependent kernel has the form of an inverse spherical Hankel transform The belonging forward transform is found via orthogonality as The function j(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) where σ = σ + (1 – f)σ and q = δ0. For the Henyey-Greenstein model [8] the expansion coefficients become f = g ∀ l ∈ ℕ0. In P approximation (N odd) equations must be considered only for l = 0,...,N with the initial conditions I−1(k,t) = I+1(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 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 where q = [1,0,...,0]T and . Here the matrix I denotes the identity matrix. The complex symmetric tridiagonal matrix A = AT has the form The solution of the matrix equation Eq. (13) is formally obtained as [16] Here . The inverse Laplace transform can be performed using a well-known transform pair within the linear system theory for the so-called transition matrix Thus the time-dependent vector components become for t > 0 A simple but inefficient method for evaluation of the matrix exponential is given by the power series 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−1, where U is the transformation matrix, which contains the eigenvectors |u〉 (n = 0,1,...,N) of A. The diagonal matrix Λ =diag(λ0,λ1,...,λ) (λ ∈ ℂ) contains the eigenvalues of the equation det(A −λI) = 0. Due to the above eigenvalue decomposition equation Eq. (17) simplifies to where the matrix exponential can now be evaluated as 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 where the vector |ũ0〉 denotes the first column vector of the matrix U−1. Note that for obtaining an expression for the radiance caused by an isotropic emitting spherical shell source of the form where r′ denotes the radius of the spherical shell, the source coefficients q = δ0 in system (11) must be replaced by q = (− j) j(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 ŝ0 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 μ = ŝ0 · 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 where R equals the radius of a large sphere. Now the continuous wave number k in relation Eq. (10) becomes to discrete values ξ which are dependent on the transform order l and are given by the positive roots of the equation j(ξ) = 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) Due to the exponential decay of the time-dependent modes I(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 ξ0 = 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 λ becomes complex then is also an eigenvalue due to the symmetry to the real axis in the complex plane so that in all cases the functions I(r,t) are automatically real-valued. The fluence [5] which is given by the lowest order coefficient I0(r,t) becomes the simple form Regarding the computation time the fluence in P25 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 P3 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 σ = 0.01mm−1, σ = 1mm−1 and g = 0.9.

Fig. 2

Time-resolved fluence at r = 3.05mm for two different anisotropy factors in a scattering medium with properties σ = 0.1mm−1 and σ = 1mm−1.

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σ1) and D = 1/(3σ′)) reported in literature [8], where σ′ = (1 − g)μ and σ1 = μ + (1 − g)μ.

Fig. 3

Time-resolved fluence at r = 1.05mm in an anisotropically scattering medium with properties σ = 0.1mm−1, σ = 1mm−1 and g = 0.9.

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 σ = 0.1mm−1, σ = 1mm−1 and g = 0.8.

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