Light transport in homogeneous tissue with m-dependent anisotropic scattering I: Case’s singular eigenfunctions solution and Chandrasekhar polynomials

This paper is the first of two deriving the analytical solutions for light transport in infinite homogeneous tissue with an azimuth-dependent (m-dependent) anisotropic scattering kernel by two approaches, Case’s singular eigenfuncions expansion and Fourier transform, as well as proving the consistence of the two solutions. In this paper, Case’s method was applied and extended to the general m-dependent anisotropic scattering case. The explicit Green’s function of radiance distributions, which was regarded as the comparative standard for the equivalent solution via Fourier transform and inversion in our second accompanying paper, was expanded into a complete set of the discrete and continuous eigenfunctions. Considering that the two kinds of m-dependent Chandrasekhar orthogonal polynomials that play vital roles in these analytical solutions are very sensitive to the typical optical parameters of biological tissue as well as the degrees or orders, four numerical evaluation methods were benchmarked to find the stable, reliable and feasible numerical evaluation methods in high degrees and high orders.


Introduction
Though for decades many approximation methods, such as the spherical harmonics ( N P ) [1,2], the simplified spherical harmonics ( N SP ) [3][4][5][6], the diffusion equation ( DE ) [7][8][9][10], N F method [11,12], etc., are applied to deal with the light transport equations and make great success, the traditional but effective Case's singular eigenfunctions (CSEs) expansion solutions [13] are still regarded as the exact analytical solutions and as the criteria for testing the approximate methods. Since K. M. Case reported the singular eigenfuntions expansion method of building the elementary solutions for the radiative transport equation (RTE) and proved the orthogonality and completeness of the singular eigenfunctions in isotropic scattering [13], the Case's method has become the hot hits and soon extended to linear anisotropic scattering [14,15], general anisotropic scattering with [16,17] and without [18] azimuthal symmetry. Due to its complexity, the numerical approximation was the main efforts in quite a long time. In recent years, the Case's method was revisited by some researchers [19][20][21] to seek the more "exact" solutions in biological optics, neutron transport theory, heat transport theory, and so on.
As for the analytical solutions, the Case's method and Fourier transform method are both indispensable. For the regular infinite space, the solutions obtained from the two methods should be consistent with each other. In fact, K. M. Case [22] applied the Fourier transform method to linear transport theory early. B. D. Ganapol [23][24][25][26] explicitly proved the consistence between the CSEs solution and the Fourier transform solution in the case 0 m = . The purpose of our two papers is to extend the consistent theory suitable for the neutron transport to the azimuth-dependent or m-dependent light transport in biological tissue. So in our first paper, we firstly derived the Green's function based on the CSEs expansion for the light transport in infinite homogeneous tissue with m-dependent anisotropic scattering kernel and made a comparative standard for the Fourier transform solution which is to be discussed in our second accompanying paper.
In Case's method, the Chandrasekhar orthogonal polynomials occupy the most important and fundamental position undisputedly. Without the orthogonal polynomials, it is impossible to build the eigenfunctions and elementary solutions, or to find the discrete eigenvalues. Inönü [27] and K. M. Case [28,29] discussed the polynomials thoroughly and laid a substantial mathematical foundation for successive researchers. C. E. Siewert, etc [30-32], proved many useful identities about Chandrasekhar polynomials and gave some feasible methods to compute the polynomials with high orders and high degrees for the sake of finding the discrete eigenvalues. In fact, in our reachable scope, there are not many literatures on the Chandrasekhar polynomials, especially in the evaluations of the polynomials. Recently, Ganapol [33] focused on the Chandrasekhar polynomials and benchmarked the various computational methods for zero order, 0 m = . However, our numerical evaluation experiments showed that in biological optics, unlike the neutron transport domain or the extremely critical testing parameters, the Chandrasekhar orthogonal polynomials are very sensitive to the optical parameters of biological tissue, such as single particle albedo ϖ and anisotropy factor g . So in this paper, we extended Ganapol's methods in the case 0 m = to benchmark the m-dependent Chandrasekhar polynomials with the typical optical parameters of biological tissue to find the stable, reliable, feasible numerical evaluation methods in high degrees and high orders.

Eigenvalues and singular eigenfunctions
The light transport in the infinite homogeneous tissue with an m-dependent anisotropic scattering kernel is governed by integro-differential equation [11,12,34,35]: where ( ) where L is the truncated position of the scattering kernel and where v is the eigenvalues and ( ) where s t ϖ μ μ = is the single particle albedo, and it is obvious that 0 1 ϖ < < in absorbing tissue, and 1 ϖ = when 0 a μ = in conservative tissue.
Considering the Hobson definition [36] of the associated Legendre functions of the first kind which is used by the mathematical computing softwares, such as Matlab, Mathematica, etc., we define g v is the Chandrasekhar polynomials that we will discuss in Section 3.

Then, when
[ ] and satisfy the normalization condition Here P denotes the Cauchy principal value. The CSEs for continuous spectrum From Eq. (12), we get

Orthogonality and norms
The orthogonality and full/partial-range completeness properties of the discrete and continuous eigenfunctions for isotropic scattering [13], for linear anisotropic scattering [14,15], and for general anisotropic scattering [16][17][18] are discussed in many early classical literatures. In our case the orthogonality and completeness are easy to be extended from the general anisotropic scattering case. However, the norms of the discrete and continuous eigenfunctions are rarely expressed explicitly. The discrete and continuous eigenfunctions for specific m satisfy the following orthogonalities: and the expansion of the directive of ( ) can also be expressed as [39,40] ( )

Singular eigenfunctions solution
Considering the boundary condition and the subcritical tissue, hence 1 ϖ < , the radiance must vanish at infinity: So the solution must be separated into two half spaces: So in the infinite homogeneous tissue with m-dependent anisotropic scattering, the component of the Green's function for specific m can be rewritten as (3) we can get the expected final result Green' function for Eq. (1). This singular eigenfunctions solution is the comparative standard for the Fourier transform solution which will be discussed in our second accompanying paper.

Chandrasekhar polynomials
The Chandrasekhar polynomials ( ) m l g z play a vital role in discrete eigenvalues finding and the numerical evaluation of the solution, so many researchers had studied them for decades. In this section we revisited the Chandrasekhar polynomials for 0 m ≥ and analyzed the stability of various computational method.

Definitions and recurrences of polynomials
To avoid confusing, following the naming conventions for Legendre polynomials, ( ) From Eq. (6), the famous three-term recurrence relation is emerged The matrix generated from Eq. (34) is not a real symmetric matrix form which is not convenient to compute the eigenvalues, so the normalized associated Chandrasekhar polynomials of the first kind are introduced and then Eq. (34) changes into with the initial term satisfy the following three-term recurrence relations

Determinant and eigenvalues of the characteristic matrix
The ]. For lower degrees or lower orders approximation evaluating the polynomials by determinant is still a good choice. However, they are inconvenient to compute and analyze due to their asymmetries. The associated Chandrasekhar polynomials of the first and second kinds which have got a real symmetric determinantal form can be expressed as where the characteristic matrix  As shown when g incre about 1 excep pointed out th From another biological tiss not be found.
Next, we the typical op without loss matrix to 500 varied from 2 sorted them in fill the extra p along the orde tendency of e the tendency polynomials.
the eigenvalue nd the discrete of eigenvalues ticle albedo ϖ from 0 to 30 ransport theory g are about 0 matrix were co 0.9 respectively 1. The 2-norms of subfigure denote th n in Fig. 1

Error measurements
From above four computational methods, the resulting sequences of the associated Chandrasekhar polynomials of the first and second kinds can be evaluated. However, as many other similar polynomials, such as    Since the forward recurrence is a perfect evaluation method for the polynomials of the second kind, we may compute them firstly, and then compute the polynomials of the first kind by C-D formula. After numerical experiments, we found that when 1 z < , this method is equivalent to the method (a) and method (b), and when 1 z > it is equivalent to the method (c). But due to the rapid overflow or underflow, we don't recommend this method. Firstly, we noticed that the parameters , , , N m g ϖ had quite an influence on the arithmetic or the eigenvalues and eigenvectors of characteristic matrix. In some cases there are no eigenvalues greater than 1, and in other cases there are eigenvalues greater than 1 or less than 1, as indicated in section 3.2. After the eigenvectors were computed, we found that in many cases, when the eigenvalues are approaching the maximum and minimum eigenvalue, the corresponding eigenvectors have some zero components (or very close to zero), and the number of the zero components increases as the eigenvalues increase. These zero components is caused by the fact, as discussed in section 3.2, that many of the eigenvalues with the same position will tend to a constant as the order of the characteristic matrix increases. So the corresponding components, or polynomials with different degrees, have common zeroes.

3) Method (d) benchmark
Next, we tested many characteristic matrix with various parameters like before and got their eigenvalues and eigenvectors. On one hand, by method (d), we evaluated the values of the first kind of polynomials for the corresponding degrees and orders. On the other hand, for every eigenvalue we applied forward recurrence to the polynomials of the first and second kinds with the same parameters and without considering whether the eigenvalues is less than 1 or not. Then we find that the method (d), in fact, is equivalent to the forward recurrence method, especially when the eigenvalue is close to 0 or  As shown in Fig. 9, the x-axis, eigenvalue index, is the position number of the positive eigenvalues in ascending order, and the y-axis, error vector index, is the index of the residual errors computed by Eq. (61) and Eq. (67). In fact, though when the eigenvalue is greater than 1 the errors is not much small, as shown in the Fig. 9 (a), the results of method (d) confirm well to the C-D formula. Then, we checked the logarithmic absolute errors between the matrix of the polynomials from eigenvectors matrix by method (d) and the matrix of the polynomials by method (a). As shown in Fig. 9 (b), it is obvious that when the eigenvalue is less than 1, the errors are very small. When the eigenvalue is greater than 1, the forward recurrence is unstable, so it is no meaning to focus on it.
So method (d) is also a feasible method to evaluate the polynomials of the first kind, especially when the other methods are not available.

Conclusion
The Case's method is used to find the radiance distributions of the light transport in the infinite homogeneous tissue with m-dependent anisotropic scattering kernel, and the final analytical solution is expanded by the contributions of the discrete eigenfunctions and the continuous singular eigenfunctions. In our second accompanying paper, we will apply the Fourier transform to obtain another analytical solution of the light transport equation and prove the consistence between the two solutions in a mathematical manner.
Considering the vital importance of the two kinds of m-dependent Chandrasekhar orthogonal polynomials, we made great efforts on them. After the benchmark, the following guiding rules may be helpful: 1) The determinant eigenvalues method is identical to the forward recurrence method for two kinds of Chandrasekhar polynomials.
2) In the case 0 m = , the forward recurrence method is stable for the polynomials of the first kind only when 1 z < , and the backward recurrence method is suitable for them only when 1 z > .
3) In the case 0 m ≠ , the forward recurrence method is only valid for the polynomials of the first kind when 1 z < and not close to 1 ± , and it should be noticed that the bigger m , the worse the stability as z is approaching 1 ± , when 1 z > or close to 1 ± , the backward recurrence should be applied to them. 4) The forward recurrence is stable for the polynomials of the second kind for any z and any m , and the backward recurrence is not suitable for the polynomials of the second kind.
5) The linear system of equations method is also a feasible method for the larger degrees of the polynomials of the first kind, it is not suitable for the polynomials of the second kind. Based on the conclusions, our future wok is to establish the measurable standards or formulas for the stabilities in m-dependent case. Since the Chandrasekhar polynomials can be evaluated effectively, the solutions of the light transport integro-differential equations in one/two/three-dimension full space, half space, and slab space etc. with general scattering kernel may be benchmarked to seek the more "exact" radiance or heat distributions in various tissue.