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

Lin Wang; Meng Rong; Kaiyang Li

doi:10.1364/BOE.9.004009

1. Introduction

Though for decades many approximation methods, such as the spherical harmonics ( $P_{N}$ ) [1, 2], the simplified spherical harmonics ( $S P_{N}$ ) [3–6], the diffusion equation ( $D E$ ) [7–10], $F_{N}$ 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–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–26] explicitly proved the consistence between the CSEs solution and the Fourier transform solution in the case $m = 0$ . 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, $m = 0$ . 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 $m = 0$ 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.

2. Case’s singular eigenfunctions solution

2.1 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]:

\begin{array}{l} μ \frac{\partial ψ (x, μ, φ; μ_{0}, φ_{0})}{\partial x} + μ_{t} ψ (x, μ, φ; μ_{0}, φ_{0}) = μ_{s} \int_{0}^{2 π} \int_{- 1}^{+ 1} p (\cos γ) ψ (x, μ^{'}, φ^{'}; μ_{0}, φ_{0}) d μ^{'} d φ^{'} \\ \begin{array}{l}  \end{array} + δ (x) δ (μ - μ_{0}) δ (φ - φ_{0}) \end{array}

where

ψ (x, μ, φ; μ_{0}, φ_{0})

is the radiance of the scattering light which only depends on the position coordinate

x \in ℝ

, the direction cosine

μ

measured with the angle between the propagating radiation and the positive

x

axis, and the azimuthal angle

φ

.

μ_{0}

,

φ_{0}

are the direction cosine and azimuthal angle of the source positioned at

x = 0

, respectively.

μ_{s}

is the scattering coefficient,

μ_{a}

is the absorption coefficient, and

μ_{t} = μ_{a} + μ_{s}

is the total attenuation coefficient [35].

p (\cos γ)

is the scattering phase function where

γ

is the scattering angle. We assume that the scatters are spherical symmetric and express the phase function as spherical harmonics expansion [21, 34]:

p (\cos γ) = \frac{1}{4 π} \sum_{l = 0}^{L} (2 l + 1) f_{l} \sum_{m = - l}^{l} \frac{(l - m)!}{(l + m)!} P_{l}^{m} (μ) P_{l}^{m} (μ^{'}) e^{i m (φ - φ^{'})}

where

L

is the truncated position of the scattering kernel and

L > 1

.

f_{0} = 1

,

0 < f_{l} < 1

for

1 \leq l \leq L

, and

f_{l} = 0

for

l > L

. In biological tissue optics, the Henyey-Greenstein anisotropic scattering kernel is often applied and then

f_{l} = g^{l}

[35]. The anisotropy factor

g

is the average cosine of angle of scattering [35].

P_{l}^{m} (μ)

is the associated Legendre functions of the first kind [36].

To develop the singular eigenfunctions solution of Eq. (1), we expand the $ψ (x, μ, φ; μ_{0}, φ_{0})$ as Fourier series [18]

ψ (x, μ, φ; μ_{0}, φ_{0}) = \frac{1}{2 π} \sum_{m = - \infty}^{\infty} ψ^{m} (x, μ; μ_{0}) e^{i m φ}

Substituting Eq. (3) and $δ (φ - φ_{0}) = \frac{1}{2 π} \sum_{m = - \infty}^{\infty} e^{i m (φ - φ_{0})}$ into Eq. (1), we get

\begin{array}{l} μ \frac{\partial ψ^{m} (x, μ; μ_{0})}{\partial x} + μ_{t} ψ^{m} (x, μ; μ_{0}) = \frac{μ_{s}}{2} \sum_{l = | m |}^{L} β_{l}^{m} P_{l}^{m} (μ) \int_{- 1}^{+ 1} P_{l}^{m} (μ^{'}) ψ^{m} (x, μ^{'}; μ_{0}) d μ^{'} \\ \begin{array}{l}  \end{array} + δ (x) δ (μ - μ_{0}) e^{- i m φ_{0}} \end{array}

where

β_{l}^{m} = \frac{(l - m)!}{(l + m)!} (2 l + 1) f_{l}

.Thus, according to the Case’s method [13] and separation of variables in homogeneous equation of Eq. (4), the solution in full-space has the form:

ψ^{m} (x, μ) \equiv \exp (- μ_{t} x / v) ϕ^{m} (v, μ)

where

v

is the eigenvalues and

ϕ^{m} (v, μ)

is the corresponding eigenfunctions. Substituting Eq. (5) into the homogeneous equation of Eq. (4), we get the equation for

ϕ^{m} (v, μ)

:

(1 - \frac{μ}{v}) ϕ^{m} (v, μ) = \frac{ϖ}{2} \sum_{l = | m |}^{L} β_{l}^{m} P_{l}^{m} (μ) \int_{- 1}^{+ 1} P_{l}^{m} (μ^{'}) ϕ^{m} (v, μ^{'}) d μ^{'}

where

ϖ = μ_{s} / μ_{t}

is the single particle albedo, and it is obvious that

0 < ϖ < 1

in absorbing tissue, and

ϖ = 1

when

μ_{a} = 0

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_{l}^{m} (v) \equiv {(- 1)}^{m} \int_{- 1}^{+ 1} P_{l}^{m} (μ) ϕ^{m} (v, μ) d μ

_{$g_{l}^{m} (v)$} is the Chandrasekhar polynomials that we will discuss in Section 3.

Then, when $v \notin [- 1, 1]$ , the discrete eigenfunctions are

ϕ^{m} (v, μ) = \frac{ϖ v}{2} \frac{g^{m} (v, μ)}{v - μ}

and satisfy the normalization condition

\int_{- 1}^{1} ϕ^{m} (v, μ) {(1 - μ^{2})}^{| m | / 2} d μ = 1

where

g^{m} (v, μ) = \sum_{l = | m |}^{L} β_{l}^{m} g_{l}^{m} (v) P_{l}^{m} (μ)

So we define the function

Λ^{m} (z) \equiv 1 - \frac{ϖ z}{2} \int_{- 1}^{1} \frac{g^{m} (z, μ)}{z - μ} {(1 - μ^{2})}^{| m | / 2} d μ

and get the discrete paired real eigenvalues

\pm v_{j}^{m} (1 \leq j \leq M)

which are determined by the roots of

Λ^{m} (\pm v_{j}^{m}) = 0

,

M

is the number of positive eigenvalues [29].

Λ^{m} (z)

is a holomorphic function in the complex z-plane excluding a cut line −1 to 1 on the real axis. When

z

approaches the cut line from above

(+)

and below

(-)

, by Plemelj’s formula [37], we get

{[Λ^{m} (v)]}^{\pm} = 1 - \frac{ϖ v}{2} P \int_{- 1}^{1} \frac{g^{m} (v, μ)}{v - μ} {(1 - μ^{2})}^{| m | / 2} d μ \pm \frac{i π ϖ v}{2} g^{m} (v, v) {(1 - v^{2})}^{| m | / 2}

Here

P

denotes the Cauchy principal value.

The CSEs for continuous spectrum $v \in (- 1, 1)$ are

ϕ^{m} (v, μ) = \frac{ϖ v}{2} P \frac{g^{m} (v, μ)}{v - μ} + \frac{λ (v)}{{(1 - v^{2})}^{| m | / 2}} δ (v - μ)

where

λ (v) = 1 - \frac{ϖ v}{2} P \int_{- 1}^{1} \frac{g^{m} (v, μ)}{v - μ} {(1 - μ^{2})}^{| m | / 2} d μ

From Eq. (12), we get

{[Λ^{m} (v)]}^{+} + {[Λ^{m} (v)]}^{-} = 2 λ (v)

{[Λ^{m} (v)]}^{+} {[Λ^{m} (v)]}^{-} = λ^{2} (v) + {[\frac{π ϖ v}{2} g^{m} (v, v) {(1 - v^{2})}^{| m | / 2}]}^{2}

2.2 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–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:

\int_{- 1}^{1} ϕ^{m} (v_{i}^{m}, μ) ϕ^{m} (v_{j}^{m}, μ) μ d μ = N_{0}^{m} (v_{j}^{m}) δ_{i j}

\int_{- 1}^{1} ϕ^{m} (v, μ) ϕ^{m} (v^{'}, μ) μ d μ = N^{m} (v) δ (v - v^{'})

\int_{- 1}^{1} ϕ^{m} (v_{j}, μ) ϕ^{m} (v, μ) μ d μ = 0

where

v_{i}^{m}

,

v_{j}^{m}

are any discrete eigenvalues,

v

,

v^{'}

are in continuous spectrum,

v, v^{'} \in (- 1, 1)

.

δ_{i j}

is the Kronecker delta,

δ_{i j} = 1

for

i = j

, and otherwise

δ_{i j} = 0

.

δ (v - v^{'})

is the Dirac delta.

N_{0}^{m} (v_{j}^{m})

is the norms of the discrete eigenfunctions

\begin{array}{l} N_{0}^{m} (\pm v_{j}^{m}) = \frac{\pm ϖ v_{j}^{m}}{2} [ϖ {(v_{j}^{m})}^{2} \sum_{l = | m |}^{L} β_{l}^{m} g_{l}^{m} (v_{j}^{m}) \sum_{k = | m |}^{L} β_{k}^{m} g_{k}^{m} (v_{j}^{m}) \frac{γ_{l m} δ_{m 0}}{{(v_{j}^{m})}^{2} - 1} \\ \begin{matrix}  \end{matrix} - \sum_{l = | m |}^{L} β_{l}^{m} {[g_{l}^{m} (v_{j}^{m})]}^{2} - (1 + δ_{m 0}) v_{j}^{m} \sum_{l = | m |}^{L} β_{l}^{m} g_{l}^{m} (v_{j}^{m}) \sum_{k = 0}^{[\frac{l - 1}{2}]} (2 l - 4 k - 1) g_{l - 2 k - 1}^{m} (v_{j}^{m})] \end{array}

where

δ_{m 0}

is the Kronecker delta,

δ_{m 0} = 1

for

m = 0

, and otherwise

δ_{m 0} = 0

.

γ_{l m} = 1

for

l + k

is even, and

γ_{l m} = v_{j}

for

l + k

is odd.

[(l - 1) / 2]

denotes the nearest integer less than or equal to

(l - 1) / 2

. The proof of

N_{0}^{m}

needs the integral relation for

z

that is not on the cut line −1 to 1 on the real axis [38]

\int_{- 1}^{1} \frac{P_{l}^{m} (μ) P_{k}^{m} (μ)}{z - μ} d μ = {\begin{matrix} 2 P_{l}^{m} (z) Q_{k}^{m} (z) & l \leq k \\ 2 Q_{l}^{m} (z) P_{k}^{m} (z) & l > k \end{matrix}

and the expansion of the directive of

P_{l}^{m} (x)

\frac{d P_{l}^{m} (x)}{d x} = \frac{1 + δ_{m 0}}{2} \sum_{k = 0}^{[\frac{l - 1}{2}]} (2 l - 4 k - 1) P_{l - 2 k - 1}^{m} (x)

_{$Q_{l}^{m} (z)$}is the associated Legendre functions of the second kind. When

m = 0

, it is easy to validate that Eq. (20) is really identical to the results by Mika [16] and by Van Den [17] for m-independent general anisotropic scattering.

N_{0}^{m} (\pm v_{j}^{m})

can also be expressed as [39, 40]

N_{0}^{m} (\pm v_{j}^{m}) = \pm \frac{1}{2} ϖ {(v_{j}^{m})}^{2} g^{m} (v_{j}^{m}, v_{j}^{m}) {[{(v_{j}^{m})}^{2} - 1]}^{- | m | / 2} {\frac{d Λ^{m} (z)}{d z} |}_{z = v_{j}^{m}}

_{$N^{m} (v)$} in Eq. (18) is the norm of eigenfuncions for continuous spectrum [39]

N^{m} (v) = \frac{v {[Λ^{m} (v)]}^{+} {[Λ^{m} (v)]}^{-}}{{(1 - v^{2})}^{| m |}}

2.3 Singular eigenfunctions solution

Now, by Case’s elementary solution theory [13], we can write the singular eigenfunctions elementary solution for homogeneous equation of Eq. (4):

\begin{array}{l} ψ^{m} (x, μ) = \sum_{j = 1}^{M} a_{+}^{m} ϕ^{m} (v_{j}^{m}, μ) \exp (- μ_{t} x / v_{j}^{m}) + \sum_{j = 1}^{M} a_{-}^{m} ϕ^{m} (- v_{j}^{m}, μ) \exp (μ_{t} x / v_{j}^{m}) \\ \begin{array}{l}  \end{array} + \int_{- 1}^{1} A^{m} (v) ϕ^{m} (v, μ) \exp (- μ_{t} x / v) d v \end{array}

Considering the boundary condition

ψ^{m} (0 +, μ) - ψ^{m} (0 -, μ) = \frac{1}{μ} δ (x) δ (μ - μ_{0}) e^{- i m φ_{0}}

and the subcritical tissue, hence

ϖ < 1

, the radiance must vanish at infinity:

\lim_{| x | \to \infty} ψ^{m} (x, μ) = 0

So the solution must be separated into two half spaces:

ψ^{m} (x, μ; μ_{0}) = {\begin{cases} \sum_{j = 1}^{M} a_{+}^{m} ϕ^{m} (v_{j}^{m}, μ) \exp (- μ_{t} x / v_{j}^{m}) \\ \begin{matrix}  \end{matrix} + \int_{0}^{1} A^{m} (v) ϕ^{m} (v, μ) \exp (- μ_{t} x / v) d v \begin{array}{l}  \end{array} (x > 0) \\ - \sum_{j = 1}^{M} a_{-}^{m} ϕ^{m} (- v_{j}^{m}, μ) \exp (μ_{t} x / v_{j}^{m}) \\ \begin{matrix}  \end{matrix} - \int_{- 1}^{0} A^{m} (v) ϕ^{m} (v, μ) \exp (- μ_{t} x / v) d v \begin{array}{l}  \end{array} (x < 0) \end{cases}

Using the orthogonalities Eq. (17) ~Eq. (19), the coefficients $a_{\pm}^{m}$ and $A^{m} (v)$ can be easily deduced:

a_{j \pm}^{m} = \frac{1}{N_{0}^{m} (\pm v_{j}^{m})} ϕ^{m} (\pm v_{j}^{m}, μ_{0}) e^{- i m φ_{0}}

A (v) = \frac{1}{N^{m} (v)} ϕ^{m} (v, μ_{0}) e^{- i m φ_{0}}

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

ψ^{m} (x, μ; μ_{0}) = e^{- i m φ_{0}} {\begin{cases} \sum_{j = 1}^{M} \frac{ϕ^{m} (v_{j}^{m}, μ_{0}) ϕ^{m} (v_{j}^{m}, μ)}{N_{0}^{m} (v_{j}^{m})} \exp (- μ_{t} x / v_{j}^{m}) \\ \begin{matrix}  \end{matrix} + \int_{0}^{1} \frac{ϕ^{m} (v, μ_{0}) ϕ^{m} (v, μ)}{N^{m} (v)} \exp (- μ_{t} x / v) d v \begin{array}{l}  \end{array} (x > 0) \\ - \sum_{j = 1}^{M} \frac{ϕ^{m} (- v_{j}^{m}, μ_{0}) ϕ^{m} (- v_{j}^{m}, μ)}{N_{0}^{m} (- v_{j}^{m})} \exp (μ_{t} x / v_{j}^{m}) \\ \begin{matrix}  \end{matrix} - \int_{- 1}^{0} \frac{ϕ^{m} (v, μ_{0}) ϕ^{m} (v, μ)}{N^{m} (v)} \exp (- μ_{t} x / v) d v \begin{array}{l}  \end{array} (x < 0) \end{cases}

From Eq. (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.

3. Chandrasekhar polynomials

The Chandrasekhar polynomials $g_{l}^{m} (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 $m \geq 0$ and analyzed the stability of various computational method.

3.1 Definitions and recurrences of polynomials

To avoid confusing, following the naming conventions for Legendre polynomials, $g_{l} (z)$ and $ρ_{l} (z)$ are called the Chandrasekhar polynomials of the first and second kinds respectively when $m = 0$ , $g_{l}^{m} (z)$ and $ρ_{l}^{m} (z)$ the associated Chandrasekhar polynomials of the first and second kinds respectively when $m > 0$ .

Let’s back to the definition the associated Chandrasekhar polynomials of the first kind in Eq. (7). Obviously, when $l = m$ , we get the first term for specific $m$ , $g_{m}^{m} (z) = (2 m - 1)!!$ 。Due to the attributes of $P_{l}^{m} (x)$ , $g_{l}^{m} (z)$ satisfies the following relations:

g_{l}^{- m} (z) = {(- 1)}^{m} \frac{(l - m)!}{(l + m)!} g_{l}^{m} (z)

g_{l}^{m} (- z) = {(- 1)}^{l - m} g_{l}^{m} (z)

From Eq. (6), the famous three-term recurrence relation is emerged

z h_{l} g_{l}^{m} (z) - (l - m + 1) g_{l + 1}^{m} (z) - (l + m) g_{l - 1}^{m} (z) = 0

where

h_{l} = (2 l + 1) (1 - ϖ f_{l}), m \geq 0, l \geq m

. When

m = 0

, Eq. (34) is the three-term recurrence relation for the associated Chandrasekhar polynomials of the first kind.

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

{\hat{g}}_{l}^{m} (z) = \sqrt{\frac{(l - m)!}{(l + m)!}} g_{l}^{m} (z)

and then Eq. (34) changes into

z h_{l} {\hat{g}}_{l}^{m} (z) - \sqrt{{(l + 1)}^{2} - m^{2}} {\hat{g}}_{l + 1}^{m} (z) - \sqrt{l^{2} - m^{2}} {\hat{g}}_{l - 1}^{m} (z) = 0

with the initial term

{\hat{g}}_{m}^{m} (z) = \frac{(2 m - 1)!!}{\sqrt{(2 m)!}} = \frac{\sqrt{(2 m)!}}{2^{m} m!}

Similarly we define the associated Chandrasekhar polynomials of the second kind $ρ_{l}^{m} (z)$ by three-term recurrence relations

z h_{l} ρ_{l}^{m} (z) - (l - m + 1) ρ_{l + 1}^{m} (z) - (l + m) ρ_{l - 1}^{m} (z) = - ρ_{m + 1}^{m} (z) δ_{l m}^{}

with initial two terms

ρ_{m}^{m} (z) = 0

and

ρ_{m + 1}^{m} (z) = z / (2 m - 1)!!

. So the normalized associated Chandrasekhar polynomials of the second kind

{\hat{ρ}}_{l}^{m} (z) = \sqrt{\frac{(l - m)!}{(l + m)!}} ρ_{l}^{m} (z)

satisfy the following three-term recurrence relations

z h_{l} {\hat{ρ}}_{l}^{m} (z) - \sqrt{{(l + 1)}^{2} - m^{2}} {\hat{ρ}}_{l + 1}^{m} (z) - \sqrt{l^{2} - m^{2}} {\hat{ρ}}_{l - 1}^{m} (z) = - \sqrt{2 m + 1} {\hat{ρ}}_{m + 1}^{m} (z) δ_{l m}

The Christoffel–Darboux (C-D) identity or Liouville-Ostrogradski (L-O) formula [27, 33, 34] between $g_{l}^{m} (z)$ and $ρ_{l}^{m} (z)$ when $m \geq 0$ are

g_{l - 1}^{m} (z) ρ_{l}^{m} (z) - g_{l}^{m} (z) ρ_{l - 1}^{m} (z) = \frac{(l + m - 1)!}{(l - m)!} \frac{ρ_{m + 1}^{m} (z) g_{m}^{m} (z)}{(2 m)!}

So for the normalized ${\hat{g}}_{l}^{m} (z)$ and ${\hat{ρ}}_{l}^{m} (z)$ we have

{\hat{g}}_{l - 1}^{m} (z) {\hat{ρ}}_{l}^{m} (z) - {\hat{g}}_{l}^{m} (z) {\hat{ρ}}_{l - 1}^{m} (z) = \frac{\sqrt{(2 m + 1)} {\hat{ρ}}_{m + 1}^{m} (z) {\hat{g}}_{m}^{m} (z)}{\sqrt{l^{2} - m^{2}}}

Equation (41) and Eq. (42) play an indispensable part in the following computation methods and in the derivation of the Fourier transform solution. It is noted that the initial terms of two kinds of polynomials, actually may be any nonzero constants. As shown in the following sections, we will let ${\hat{g}}_{m}^{m} (z) = 1$ without loss of generality. The first term of ${\hat{ρ}}_{l}^{m} (z)$ , ${\hat{ρ}}_{m}^{m} (z)$ , is usually set to zero, and the initial nonzero term, ${\hat{ρ}}_{m + 1}^{m} (z)$ , satisfies ${\hat{ρ}}_{m + 1}^{m} (z) {\hat{g}}_{m}^{m} (z) = z$ for convenience.

3.2 Determinant and eigenvalues of the characteristic matrix

The Chandrasekhar polynomials of the first and second kinds were expressed as the form of determinant by Inönü [27], Ganapol [23]. The explicit determinantal representation of the associated Chandrasekhar polynomials of first kind $g_{l}^{m} (z)$ was given by McCormick [32, 39]. 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

{\hat{g}}_{l}^{m} (z) = \frac{\sqrt{(2 m)!} {\hat{g}}_{m}^{m} (z)}{\sqrt{(l - m)! (l + m)!}} | A_{l - 1}^{m} (z) |

{\hat{ρ}}_{l}^{m} (z) = \frac{\sqrt{(2 m + 1)!} {\hat{ρ}}_{m + 1}^{m} (z)}{\sqrt{(l - m)! (l + m)!}} | A_{l - 1}^{m + 1} (z) |

respectively, where matrix

A_{l}^{m} (z) = (\begin{matrix} z h_{m} & \sqrt{{(m + 1)}^{2} - m^{2}} \\ \sqrt{{(m + 1)}^{2} - m^{2}} & z h_{m + 1} \\ ⋱ \\ z h_{l - 1} & \sqrt{l^{2} - m^{2}} \\ \sqrt{l^{2} - m^{2}} & z h_{l} \end{matrix})

The subscript $l$ denotes the maximum degree of $h$ in the lower right corner of the matrix, and the superscript $m$ denotes the minimum degree of $h$ in the upper left corner of the matrix. $| \cdot |$ denotes the determinant of the matrix. It is obvious that ${\hat{g}}_{l}^{m} (z)$ and ${\hat{ρ}}_{l}^{m} (z)$ are both polynomials of $l - m$ degrees of $z$ .

_{$A_{l}^{m} (z)$} is square matrix of order $l - m + 1$ , so

| A_{l}^{m} (z) | = (\prod_{j = m}^{l} h_{j}) | B_{l}^{m} + z I |

Where the characteristic matrix

B_{l}^{m} = (\begin{matrix} 0 & b_{m + 1}^{m} \\ b_{m + 1}^{m} & 0 & b_{m + 2}^{m} \\ b_{m + 2}^{m} & ⋱ \\ 0 & b_{l}^{m} \\ b_{l}^{m} & 0 \end{matrix})

where

b_{j}^{m} = \sqrt{(j^{2} - m^{2}) / (h_{j - 1} h_{j})}, (j = m + 1, m + 2, \dots, l)

and

I

is

l - m + 1

identity matrix. Because of the property of the real symmetry, the eigenvalues of

B_{l}^{m}

are real, simple,

\pm

paired, and

B_{l}^{m}

can be easily converted to diagonal matrix.

Because the eigenvalues of characteristic matrix $B_{l}^{m}$ will be the key in our following evaluations and the discrete eigenvalues finding in Case’s method, here we focused on the distributions of eigenvalues. Firstly, we were interested in how the parameters of $B_{l}^{m}$ , such as the single particle albedo $ϖ$ and anisotropy factor $g$ , affect the greatest eigenvalues.

We set $m$ from 0 to 30 with step 5 and built a $1001 \times 1001$ square matrix. Different from the neutron transport theory or the critical parameters given by Ganapol [33], in biological tissue, $ϖ$ and $g$ are about 0.7 ~0.9(such as epidermis and dermis, $g = 0.78$ [35]). So the 2-norms of the matrix were computed with anisotropy factor $g$ from 0.7 to 1 with step 0.02 and $ϖ = 0.7, 0.8, 0.9$ respectively.

As shown in Fig. 1, the greatest eigenvalues are very sensitive to the anisotropy factor $g$ when $g$ increases up to about 0.85. Furthermore, when $g \leq 0.85$ , the greatest eigenvalues are about 1 except in $m = 0$ case. In fact, this may be a bad news. Inönü [27] and Machida [12] pointed out that the discrete eigenvalues in Case’s method are the zeroes of $g_{l}^{m} (z)$ as $l \to \infty$ . From another point of view, the discrete eigenvalues must be greater than 1. So in some biological tissue there is the potential danger that these expected discrete eigenvalues would not be found.

Fig. 1 The 2-norms of the characteristic matrix vary with anisotropy factor $g$ , the curves in each subfigure denote the 2-norms when $m = 0, 5, 10, 15, 20, 25, 30$ from top to bottom

Download Full Size | PDF

Next, we benchmarked the distributions and variation tendency of the eigenvalues with the typical optical parameters of biological tissue. We let $m = 0$ , $ϖ = 0.8$ and $g = 0.7, 0.8, 0.9$ without loss of generality, and then enlarged the characteristic matrix from $2 \times 2$ square matrix to $500 \times 500$ square matrix successively. That is to say, the degree of the polynomials varied from 2 to 500. For each square matrix, we computed its positive eigenvalues and sorted them in descending order. If the number of positive eigenvalues is less than 500, then fill the extra position with zeroes. Then we observed the variation tendency of the eigenvalues along the order of the square matrix increasing. Figure 2 shows the distributions and variation tendency of eigenvalues at some sampled positions, and each curve in each subfigure denotes the tendency of the eigenvalues with the same position but the different degrees of polynomials.

Fig. 2 The distributions and variation tendency of the eigenvalues at some sampled positions, each curve in each subfigure corresponds to a certain sampled position

Download Full Size | PDF

Though for different $l$ , the number of eigenvalues of the matrix is, of course, different. But as shown in Fig. 2, the norms or the greatest eigenvalues (the top curve) will approach to a constant as $l$ increases. So do the other eigenvalues. This may be a good news. For a very large $l$ , let $g_{l + 1}^{m} (z) \equiv 0$ , as shown in the following linear system of equations method, and then the eigenvalues of $B_{l + 1}^{m}$ are just the zeroes of $g_{l + 1}^{m} (z)$ which are the discrete eigenvalues and can easily be found in the lower degrees of polynomials. This indicates that we can find Case’s discrete eigenvalues with a not much bigger $l$ .

3.3 Computation of the polynomials

Similar to the case $m = 0$ [33], there are various methods to compute the Chandrasekhar polynomials for high degrees or high orders. Here, we benchmarked the following four methods and compared the resulting sequences to evaluate the round off errors.

3.3.1 Computational method

(a) Determinant eigenvalues method
Theoretically, Eq. (43), Eq. (44) can be used to evaluate ${\hat{g}}_{l}^{m} (z)$ and ${\hat{ρ}}_{l}^{m} (z)$ , but the factorials will soon lead to overflow or underflow for high orders or high degrees. We rewrite them as
${\hat{g}}_{l}^{m} (z) = {\hat{g}}_{m}^{m} (z) \prod_{j = 0}^{l - m - 1} \frac{h_{m + j} (z + λ_{g, j}^{m})}{\sqrt{(j + 1) (j + 1 + 2 m)}}$ ${\hat{ρ}}_{l}^{m} (z) = {\hat{ρ}}_{m + 1}^{m} (z) \prod_{j = 1}^{l - m - 1} \frac{h_{m + j} (z + λ_{ρ, j}^{m})}{\sqrt{(j + 1) (j + 1 + 2 m)}}$
where $λ_{g, j}^{m}$ is the $j^{t h}$ eigenvalue of $B_{l - 1}^{m}$ and $λ_{ρ, j}^{m}$ the $j^{t h}$ eigenvalue of $B_{l - 1}^{m + 1}$ . Equation (48) and Eq. (49) are optimized to fit the numerical computation because each term in the product is not too big or too small.
(b) Forward recurrence method
The forward recurrence method is the most intuitive method of evaluating the polynomials because of the explicit form of the three-term recurrence relations. From Eq. (36) and Eq. (40), the forward recurrence formulas for ${\hat{g}}_{l}^{m} (z)$ and ${\hat{ρ}}_{l}^{m} (z)$ are
${\hat{d}}_{j}^{m} (z) = \frac{z h_{j - 1} {\hat{d}}_{j - 1}^{m} (z) - \sqrt{{(j - 1)}^{2} - m^{2}} {\hat{d}}_{j - 2}^{m} (z)}{\sqrt{j^{2} - m^{2}}}$
where $j = m + 2, m + 3, \dots, N, m \geq 0$ and ${\hat{d}}_{j}^{m} (z)$ are ${\hat{g}}_{j}^{m} (z)$ or ${\hat{ρ}}_{j}^{m} (z)$ . With the appropriate initial terms the resulting sequences will be got without any difficulties.
(c) Backward recurrence method
In the backward recurrence method, we set the degree $l$ from $m$ to a certain large integer $N$ and force, a ${\hat{g}}_{m + N + 1}^{m} (z) \equiv 0$ nd then evaluated all the ${\hat{g}}_{l}^{m} (z)$ in the reverse direction. From Eq. (36), we get
${\hat{g}}_{j}^{m} (z) = \frac{z h_{j + 1} {\hat{g}}_{j + 1}^{m} (z) - \sqrt{{(j + 2)}^{2} - m^{2}} {\hat{g}}_{j + 2}^{m} (z)}{\sqrt{{(j + 1)}^{2} - m^{2}}}$
where $j = m, m + 1, \dots, m + N - 1$ . Since ${\hat{g}}_{N}^{m} (z)$ is unknown, we define
${\hat{g}}_{j}^{m} (z) = γ_{j}^{m} {\hat{g}}_{m + N}^{m} (z)$
So the coefficients $γ_{j}^{m}$ satisfy
$γ_{j}^{m} (z) = \frac{z h_{j + 1} γ_{j + 1}^{m} (z) - \sqrt{{(j + 2)}^{2} - m^{2}} γ_{j + 2}^{m} (z)}{\sqrt{{(j + 1)}^{2} - m^{2}}}$
and $γ_{m + N}^{m} (z) = 1$ , $γ_{m + N + 1}^{m} (z) = 0$ . When the backward recurrence iterated to $j = m$ without overflows or underflows, we get
${\hat{g}}_{m}^{m} (z) = γ_{m}^{m} (z) {\hat{g}}_{N}^{m} (z)$
So all the terms will be got
${\hat{g}}_{j}^{m} (z) = \frac{γ_{j}^{m}}{γ_{m}^{m} (z)} {\hat{g}}_{m}^{m} (z)$
Since we get the resulting sequence of the associated Chandrasekhar polynomials of the first kind ${\hat{g}}_{l}^{m} (z)$ , the associated Chandrasekhar polynomials of the second kind ${\hat{ρ}}_{l}^{m} (z)$ may be got by C-D identity Eq. (42) along forward direction. If we are mandatory to apply the backward recurrence to ${\hat{ρ}}_{l}^{m} (z)$ ignoring the reasonability and validity, the iteration should be end at $j = m + 1$ because ${\hat{ρ}}_{m}^{m} (z) = 0$ .
(d) Linear system of equations method
In linear system of equations method, we still set the degree $l$ from $m$ to a certain large integer $m + N$ and assembled the $N + 1$ equations into a linear system
$[B_{m + N}^{m} (z) - z I] G_{m + N}^{m} (z) = {\bar{G}}_{m + N + 1}^{m} (z)$
where vector $G_{m + N}^{m} (z) = {(a_{m}, a_{m + 1}, \dots, a_{m + N})}^{T}$ , and vector ${\bar{G}}_{m + N + 1}^{m} (z) = {(0, 0, \dots, - c a_{m + N + 1})}^{T}$ . Here
$a_{j} = \sqrt{h_{j}} {\hat{g}}_{j}^{m} (z), j = m, m + 1, \dots, m + N + 1$
and $c$ is a certain coefficient.
Now, let ${\hat{g}}_{m + N + 1}^{m} (z) \equiv 0$ , so Eq. (56) changes to homogeneous linear equations
$[B_{m + N}^{m} - z I] G_{m + N}^{m} (z) = 0$
So, according to the theories of linear algebra, Eq. (58) has a nontrivial solution only when the determinant matrix of the coefficient equals to zero
$| B_{m + N}^{m} - z I | = 0$
Then $z_{i}, (i = 0, 1, \dots, N)$ that satisfy Eq. (59) are just the paired eigenvalues of $B_{m + N}^{m}$ and corresponding $G_{m + N}^{m} (z_{i})$ are eigenvectors. For various $z_{i}$ , the ratio of $a_{j}$ to $\sqrt{h_{j}} {\hat{g}}_{j}^{m} (z_{i})$ keeps a constant which can be determined by the known initial term ${\hat{g}}_{m}^{m} (z_{i})$ , so we get the approximation
${\hat{g}}_{j}^{m} (z_{i}) = \sqrt{\frac{h_{m}}{h_{j}}} \frac{a_{j}}{a_{m}} {\hat{g}}_{m}^{m} (z_{i})$
Different from the other three methods, here we only get the approximation values of ${\hat{g}}_{l}^{m} (z)$ at the specific eigenvalue $z_{i}$ .

3.3.2 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 $P_{l}^{m} (z)$ or $Q_{l}^{m} (z)$ , each computational method represents different stability for different $z$ because of the inescapable round off errors.

To measure the errors when ${\hat{g}}_{l}^{m} (z)$ and ${\hat{ρ}}_{l}^{m} (z)$ are evaluated independently, such as method (a) and (b), we define the first residual vector

R_{1, j}^{m} (z) \equiv {\hat{g}}_{j - 1}^{m} (z) {\hat{ρ}}_{j}^{m} (z) - {\hat{g}}_{j}^{m} (z) {\hat{ρ}}_{j - 1}^{m} (z) - \frac{\sqrt{2 m + 1} {\hat{ρ}}_{m + 1}^{m} (z) {\hat{g}}_{m}^{m} (z)}{\sqrt{j^{2} - m^{2}}}

where

j = m + 1, m + 2, \dots, l

. If the computational method is stable,

{\hat{g}}_{l}^{m} (z)

and

{\hat{ρ}}_{l}^{m} (z)

will be evaluated precisely, and then they should confirm well to Eq. (42) and

R_{1, j}^{m} (z)

should be equal or very close to zero.

When ${\hat{g}}_{l}^{m} (z)$ is evaluated firstly and then ${\hat{ρ}}_{l}^{m} (z)$ by Eq. (42), we can define the second residual vector by the three-term recurrence of ${\hat{ρ}}_{l}^{m} (z)$

R_{2, j}^{m} (z) \equiv z h_{j} {\hat{ρ}}_{j}^{m} (z) - \sqrt{{(j + 1)}^{2} - m^{2}} {\hat{ρ}}_{j + 1}^{m} (z) - \sqrt{j^{2} - m^{2}} {\hat{ρ}}_{j - 1}^{m} (z)

where

j = m + 1, m + 2, \dots, l - 1

. Similarly, when

{\hat{ρ}}_{l}^{m} (z)

is evaluated firstly and then

{\hat{g}}_{l}^{m} (z)

by Eq. (42), we define the third residual vector by the three-term recurrence of

{\hat{g}}_{l}^{m} (z)

R_{3, j}^{m} (z) \equiv z h_{j} {\hat{g}}_{j}^{m} (z) - \sqrt{{(j + 1)}^{2} - m^{2}} {\hat{g}}_{j + 1}^{m} (z) - \sqrt{j^{2} - m^{2}} {\hat{g}}_{j - 1}^{m} (z)

For the same kind of resulting sequences of ${\hat{g}}_{l}^{m} (z)$ or ${\hat{ρ}}_{l}^{m} (z)$ from different computational method, we define the absolute error vector

R_{3, j}^{m} (z) \equiv {\hat{d}}_{j, 1}^{m} (z) - {\hat{d}}_{j, 2}^{m} (z)

where

{\hat{d}}_{j, 1}^{m} (z)

,

{\hat{d}}_{j, 2}^{m} (z)

denotes the values of

{\hat{g}}_{j}^{m} (z)

or

{\hat{ρ}}_{j}^{m} (z)

from different computational method, and the relative error vector

R_{4, j}^{m} (z) \equiv \frac{{\hat{d}}_{j, 1}^{m} (z) - {\hat{d}}_{j, 2}^{m} (z)}{{\hat{d}}_{j, 1}^{m} (z)}

After the residual or error vectors are got, we can measure the average of errors by

D_{i}^{m} (z) \equiv - \log | R_{i, j}^{m} (z) |

or

{\bar{D}}_{i}^{m} (z) \equiv - \log {R_{i, j}^{m} (z) {[R_{i, j}^{m} (z)]}^{T}}

where

i = 1, 2, 3, 4

.

3.3.3 Numerical evaluation results

We have done a plenty of numerical experiments for the associated Chandrasekhar polynomials of the first and second kinds based on the previous computational methods and measured the errors by the previous measurements. In our numerical experiments, we set ${\hat{g}}_{m}^{m} (z) \equiv 1$ and ${\hat{ρ}}_{m}^{m} (z) \equiv z$ .

1) Method (a) and method (b) benchmark
Based on the experience in the case $m = 0$ , firstly, we set a certain $| z | < 1$ . We let $N = 1000$ , $m = 15$ then $l = 15, 16, \dots, 1015$ and let anisotropy factor $g = 0.8$ , single particle albedo $ϖ = 0.9$ , and $z_{0} = 0.48$ , and then the results are shown in Fig. 3.

As shown in Fig. 3 (a), the sequence of $g_{l}^{m} (z_{0})$ , and (b), the sequence of $ρ_{l}^{m} (z_{0})$ , even for a larger $N$ , the two methods can complete successfully and the envelopes of two kinds of polynomials decrease as the degree $l$ increases. Figure 3(c) shows the logarithmic absolute errors $D_{3}^{m} (z_{0})$ of the two sequences of $g_{l}^{m} (z_{0})$ from method (a) and method (b) respectively. Figure 3(d) shows the $D_{3}^{m} (z_{0})$ for $ρ_{l}^{m} (z_{0})$ .The extremely small errors show that the resulting sequences of the associated Chandrasekhar polynomials of the first/second kind computed by method (a) and (b) are almost identical. Figure 3(e) and (f) are logarithmic residual errors $D_{1}^{m} (z_{0})$ of $g_{l}^{m} (z_{0})$ and $ρ_{l}^{m} (z_{0})$ respectively, where each polynomials is evaluated by method (a) and (b). It is apparent that they agree with C-D formula perfectly.

Fig. 3 The result sequences and errors of the first and second kinds of Chandrasekhar polynomials as the degree (x-axis) increases when $z = 0.48$ by method (a) or (b).

Download Full Size | PDF

Next, since when $| z | < 1$ , the method (a) and method (b) seem to be stable, we iterated $z$ from 0.01 to 0.99 with step 0.02 and $m$ from 0 to 30 with step 5. To save the computational time, we set $N = 599$ and focused on the average errors measurement. For each pair of $z$ and $m$ , we computed $g_{l}^{m} (z)$ and $ρ_{l}^{m} (z)$ with the degree $l$ from $m$ to $m + N$ by method (a) and (b), and then got four resulting sequences. For each method, we got the residual vector $R_{1, j}^{m} (z), j = m, m + 1, \dots, m + N$ , and then the average residual errors . Figure 4(a) and (b) describe the changes of the average residual errors with $z$ and $m$ . When $m = 0$ , the errors is still extremely small even when $z \to 1$ . So in Ref [33], Ganapol can safely conclude that the forward recurrence is stable for $g_{l} (z)$ when $| z | < 1$ . However, as shown in Fig. 4 (a) and (b), for a specific $m$ , the errors will be increasing as $z$ is increasing, and for a not very big order $m$ , the logarithmic average residual errors will deteriorate rapidly to under zero. That is to say, the resulting sequences do not confirm to the C-D formula at all. Similarly, we computed the absolute errors $R_{3, j}^{m} (z)$ by Eq. (64), and then the average absolute errors ${\bar{D}}_{3}^{m} (z)$ shown in Fig. 4 (c) and (d). The changes of the average absolute errors are similar to the average residual errors.

Fig. 4 The average errors of the first (left) and second (right) kinds of Chandrasekhar polynomials as $z$ from 0 to 1 with step 0.2 and m from 0 to 30 with step 5 when $N = 599$

Download Full Size | PDF

As shown in Fig. 4, when $| z | < 1$ and $z$ is not approaching $\pm 1$ , for a not very big $m$ , we can conclude that method (a) and method (b) are stable for two kinds of polynomials. But we must be also clearly aware that in certain aforementioned cases the stabilities of these two methods are broken remarkably and we can’t determine that which method is unstable or both methods are unstable. Then the backward recurrence must come on stage. And when $| z | > 1$ , the two methods will soon lead to overflow, and whether they are stable is to be determined.

Since in a non-strict sense the method (a) and (b) are identical and stable when $| z | < 1$ , we can easily evaluate the Chandrasekhar polynomials with variable $z$ and observe the complex curves to verify the abstract attributes of zeroes that summarized by Inönü [27] which are common to all orthogonal polynomials in general [36]. For example, as shown in Fig. 5, the zeroes of the five curves are interlacing, i.e., between two consecutive zeroes of ${\hat{g}}_{l}^{m} (z)$ , there must be one and only one zero of ${\hat{g}}_{l - 1}^{m} (z)$ , and vice versa.

Fig. 5 The curves of the first (left) and second (right) kinds of Chandrasekhar polynomials as $z$ from −1 to 1. $g (l, m)$ denotes ${\hat{g}}_{l}^{m} (z)$ , and $ρ (l, m)$ denotes ${\hat{ρ}}_{l}^{m} (z)$

Download Full Size | PDF

2) Method (b) and method (c) benchmark
When $| z | < 1$ , since in most cases the forward recurrence methods for two kinds of polynomials are stable and reliable, we introduced the backward recurrence and applied to the two kinds of polynomials to determine the stability of method (c) when $| z | < 1$ . We iterated $z$ from 0.01 to 0.99 with step 0.02 and $m$ from 0 to 30 with step 5 as previous parameters and computed, for each kind of polynomials, the average logarithmic absolute errors of the two sequences from method (b) and (c) by Eq. (64) and Eq. (66). As shown in Fig. 6, for two kinds of polynomials, the average absolute errors are too large to be accepted even in the case $m = 0$ . So the backward recurrence method is not suitable for the two kinds of polynomials when $| z | < 1$ .

Thus, when $| z | < 1$ , for each kind of polynomials, we got two series of resulting sequences, one is from forward recurrence, and the other is from backward recurrence. So the average residual errors can be evaluated by Eq. (61) and Eq. (66) to determine how they confirm to C-D formula. The average residual errors for ${\hat{g}}_{l}^{m} (z)$ and ${\hat{ρ}}_{l}^{m} (z)$ which are both from method (a) are discussed before. The other combinations are shown in Fig. 7. As shown in Fig. 7, if the backward recurrence is applied to the polynomials of the second kind, the errors are still very large. It can also be seen that the backward recurrence method is not suitable for the second kind of polynomials when $| z | < 1$ . Though Fig. 7 (c) shows that the errors are very small, considering the results of Fig. 6 and comparing the resulting sequences, we still believe that the backward recurrence method is not suitable for the first kind of polynomials when $| z | < 1$ .

Fig. 6 The average logarithmic absolute errors of the two sequences of the first (left) and second (right) kind of Chandrasekhar polynomials as z from 0 to 1 with step 0.2 and m from 0 to 30 with step 5 when $N = 599$ .

Download Full Size | PDF

Fig. 7 The average logarithmic residual errors of the first and second kind of Chandrasekhar polynomials as z from 0 to 1 with step 0.2 and m from 0 to 30 with step 5 when $N = 599$ .

Download Full Size | PDF

Thus, it is obvious that the forward recurrence method is suitable for two kinds of polynomials when $| z | < 1$ , but the backward recurrence method is not. We went on making the benchmarks when $| z | > 1$ . When $| z | > 1$ , the benchmarks will tend to suffer from the parameters severely. When the backward recurrence is applied to the first kind of polynomials the evaluation will soon exceed the machine precision. So we must adjust the appropriate parameters, such as $m, ϖ, z$ , and decrease $N$ . We let $N = 251$ , $z$ from 1.25 to 9.00 with step 0.25 and $m$ from 0 to 30 with step 5 as previous parameters.

As shown in Fig. 8 (a), the backward recurrence is applied to ${\hat{g}}_{l}^{m} (z)$ and ${\hat{ρ}}_{l}^{m} (z)$ , the very large errors show that the backward recurrence is not simultaneously suitable for two kinds of polynomials. So as shown in Fig. 8 (b), the combination, that the backward recurrence is applied to ${\hat{g}}_{l}^{m} (z)$ and the forward recurrence is applied to ${\hat{ρ}}_{l}^{m} (z)$ , is the best choice because of the extremely small errors.

Fig. 8 The average logarithmic residual errors of the first and second kind of Chandrasekhar polynomials as z from 1.25 to 9.00 with step 0.25 and m from 0 to 30 with step 5 when $N = 251$ .

Download Full Size | PDF

Thus, we can safely conclude (not show) that (1) in the case $m = 0$ , the forward recurrence method is stable for the polynomials of the first kind only when $| z | < 1$ and the backward recurrence is suitable for them only when $| z | > 1$ , (2) in the case $m \neq 0$ , the forward recurrence method is only valid for the polynomials of the first kind when $| z | < 1$ and not close to $\pm 1$ , and it should be aware that the bigger $m$ , the worse the stability as $z$ is approaching $\pm 1$ , so when $| z | > 1$ or close to $\pm 1$ , the backward recurrence should be applied to them, (3) 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.

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 $| z | < 1$ , this method is equivalent to the method (a) and method (b), and when $| z | > 1$ it is equivalent to the method (c). But due to the rapid overflow or underflow, we don’t recommend this method.

3) Method (d) benchmark
From Eq. (59), the main task in method (d) is to find the eigenvalues of $B_{m + N}^{m}$ , and the approximate values of the polynomials of the first kind at each eigenvalue $z_{i} (i = 0, 1, \dots, N)$ can be evaluate by Eq. (60). The distributions and variation tendency of the eigenvalues are discussed in section 3.2. Here, the question is whether these approximate values are accurate for $| z_{i} | > 1$ or $| z_{i} | < 1$ .
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.
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 $m = 0$ . In a typical case, we let $ϖ = 0.9, g = 0.8, m = 0,$ and $N = 251$ . There are four positive eigenvalues greater than 1. For every eigenvalue, to verify if the value vector of the polynomials of first kind from method (d) satisfies the C-D formula, we checked the logarithmic residual errors by Eq. (61) and Eq. (67).

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.

Fig. 9 The logarithmic residual errors and absolute errors in method (d).

Download Full Size | PDF

So method (d) is also a feasible method to evaluate the polynomials of the first kind, especially when the other methods are not available.

4. 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 $m = 0$ , the forward recurrence method is stable for the polynomials of the first kind only when $| z | < 1$ , and the backward recurrence method is suitable for them only when $| z | > 1$ .
3) In the case $m \neq 0$ , the forward recurrence method is only valid for the polynomials of the first kind when $| z | < 1$ and not close to $\pm 1$ , and it should be noticed that the bigger $m$ , the worse the stability as $z$ is approaching $\pm 1$ , when $| z | > 1$ or close to $\pm 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.

Funding

National Major Special Program of Scientific Instrument & Equipment Development of China (No. 2012YQ160203); Talent Plan 3551 of Wuhan City [No. WuXinGuan-(2010) 210], the Technology Innovation Fund of the National Science and Technology Ministry of China (No. 08C26224211088).

Acknowledgments

The main results of our two papers were obtained by extending derivation and benchmarks for the azimuth-independent Chandrasekhar polynomials by B D Ganapol to the m-dependent case and adapting these theories to biological optics.

Disclosures

The authors declare that there are no conflicts of interest related to this article.

References and links

1. S. Wright, M. Schweiger, and S. R. Arridge, “Reconstruction in optical tomography using the P_N approximations,” Meas. Sci. Technol. 18(1), 247–260 (2007). [CrossRef]

2. 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(19), 7364–7383 (2011). [CrossRef]

3. A. D. Klose and E. W. Larsen, “Simplified spherical harmonics methods for modeling light transport in biological tissue,” in Biomedical Optics, Technical Digest (CD) (Optical Society of America, 2006), paper MH3. [CrossRef]

4. A. Liemert and A. Kienle, “Analytical solutions of the simplified spherical harmonics equations,” Opt. Lett. 35(20), 3507–3509 (2010). [CrossRef] [PubMed]

5. H. Zheng and W. Han, “On simplified spherical harmonics equations for the radiative transfer equation,” J. Math. Chem. 49(8), 1785–1797 (2011). [CrossRef]

6. A. Liemert and A. Kienle, “Comparison between radiative transfer theory and the simplified spherical harmonics approximation for a semi-infinite geometry,” Opt. Lett. 36(20), 4041–4043 (2011). [CrossRef] [PubMed]

7. M. Machida, G. Y. Panasyuk, J. C. Schotland, and V. A. Markel, “Diffusion approximation revisited,” J. Opt. Soc. Am. A 26(5), 1291–1300 (2009). [CrossRef] [PubMed]

8. A. Liemert and A. Kienle, “Light diffusion in a turbid cylinder. I. Homogeneous case,” Opt. Express 18(9), 9456–9473 (2010). [CrossRef] [PubMed]

9. A. Liemert and A. Kienle, “Light diffusion in a turbid cylinder. II. Layered case,” Opt. Express 18(9), 9266–9279 (2010). [CrossRef] [PubMed]

10. A. Kienle and A. Liemert, “Light diffusion in turbid media of different geometries in the steady-state, frequency, and time domains,” in Biomedical Optics and 3-D Imaging, OSA Technical Digest (CD) (Optical Society of America, 2010), paper BSuD41. [CrossRef]

11. C. Devaux and C. E. Siewert, “The F_N method for radiative transfer problems without azimuthal symmetry,” Z. Angew. Math. Phys. 31(5), 592–604 (1980). [CrossRef]

12. M. Machida, “An F_N method for the radiative transport equation in three dimensions,” J. Phys. A Math. Theor. 48(32), 325001 (2015). [CrossRef]

13. K. M. Case, “Elementary solutions of the transport equation and their applications,” Ann. Phys. 9(1), 1–23 (1960). [CrossRef]

14. N. J. Mccormick and N. Joseph, “One-speed neutron transport problems in plane geometry,” Univ. of Michigan Thesis, (1965).

15. N. J. McCormick and I. Kuščer, “Half-space neutron transport with linearly anisotropic scattering,” J. Math. Phys. 6(12), 1939–1945 (1965). [CrossRef]

16. J. R. Mika, “Neutron transport with anisotropic scattering,” Nucl. Sci. Eng. 11(4), 415–427 (1961). [CrossRef]

17. G. Van den Eynde, “Neutron transport with anisotropic scattering theory and applications,” Thesis, (2005).

18. N. V. Sultanov, “Elementary solution of the neutron transport equation with anisotropic scattering,” Sov. Atom Energ. 32(6), 539–546 (1972). [CrossRef]

19. A. D. Kim, “Transport theory for light propagation in tissues,” in Biomedical Topical Meeting, OSA Technical Digest (Optical Society of America, 2004), paper SD4.

20. A. Liemert and A. Kienle, “Light transport in three-dimensional semi-infinite scattering media,” J. Opt. Soc. Am. A 29(7), 1475–1481 (2012). [CrossRef] [PubMed]

21. M. Machida, “Singular eigenfunctions for the three-dimensional radiative transport equation,” J. Opt. Soc. Am. A 31(1), 67–74 (2014). [CrossRef] [PubMed]

22. K. M. Case and R. D. Hazeltine, “Fourier transform methods in linear transport theory,” J. Math. Phys. 12(9), 1970–1980 (1971). [CrossRef]

23. B. D. Ganapol, “A consistent theory of neutral particle transport in an infinite medium,” Transp. Theory Stat. Phys. 29(1–2), 43–68 (2000). [CrossRef]

24. B. D. Ganapol, “Fourier transform transport solutions in spherical geometry,” Transp. Theory Stat. Phys. 32(5–7), 587–605 (2003). [CrossRef]

25. B. D. Ganapol, “The Fourier transform solution for the Green’s function of monoenergetic neutron transport theory,” in 2014 Hawaii University International Conferences on Science, Technology, Engineering, Math & Education (2014).

26. B. D. Ganapol, “The infinite medium Green’s function of monoenergetic neutron transport theory via Fourier transform,” Nucl. Sci. Eng. 180(2), 224–246 (2015). [CrossRef]

27. Inönü and Erdal, “Orthogonality of a set of polynomials encountered in neutron transport and radiative transfer theories,” J. Math. Phys. 11(2), 568–577 (1970). [CrossRef]

28. K. M. Case, “Orthogonal polynomials from the viewpoint of scattering theory,” J. Math. Phys. 15(12), 2166–2174 (1974). [CrossRef]

29. K. M. Case, “Scattering theory, orthogonal polynomials, and the transport equation,” J. Math. Phys. 15(7), 974–983 (1974). [CrossRef]

30. R. D. M. Garcia and C. E. Siewert, “On discrete spectrum calculations in radiative transfer,” J. Quant. Spectrosc. Radiat. Transf. 42(5), 385–394 (1989). [CrossRef]

31. R. D. M. Garcia and C. E. Siewert, “On computing the Chandrasekhar polynomials in high order and high degree,” J. Quant. Spectrosc. Radiat. Transf. 43(3), 201–205 (1990). [CrossRef]

32. C. E. Siewert and N. J. Mccormick, “Some identities for Chandrasekhar polynomials,” J. Quant. Spectrosc. Radiat. Transf. 57(3), 399–404 (1997). [CrossRef]

33. B. D. Ganapol, “Chandrasekhar polynomials and the solution to the transport equation in an infinite medium,” J. Comput. Theor. Trans. 43(1–7), 433–473 (2014). [CrossRef]

34. M. Machida, “The Green’s function for the three-dimensional linear Boltzmann equation via Fourier transform,” J. Phys. A Math. Theor. 49(17), 175001 (2016). [CrossRef]

35. A. J. Welch and M. J. C. V. Gemert, Optical-Thermal Response of Laser-Irradiated Tissue (Springer Netherlands, 2011).

36. Z. X. Wang and D. R. Guo, Special Functions (World Scientific, 1989).

37. N. I. Muskhelishvili, Singular Integral Equations: Boundary Problems of Function Theory and Their Application to Mathematical Physics (Dover Publications, 2011).

38. D. Zwillinger and V. Moll, Table of Integrals, Series, and Products (Eighth Edition) (Academic Press, 2015).

39. I. Kuščer and N. J. McCormick, “Some analytical results for radiative transfer in thick atmospheres,” J. Comput. Theor. Trans. 20(5–6), 351–381 (1974).

40. N. J. McCormick and I. Kuščer, “Bi-orthogonality relations for solving half-space transport problems,” J. Math. Phys. 7(11), 2036–2045 (1966). [CrossRef]

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

Abstract

1. Introduction

2. Case’s singular eigenfunctions solution

2.1 Eigenvalues and singular eigenfunctions

2.2 Orthogonality and norms

2.3 Singular eigenfunctions solution

3. Chandrasekhar polynomials

3.1 Definitions and recurrences of polynomials

3.2 Determinant and eigenvalues of the characteristic matrix

3.3 Computation of the polynomials

3.3.1 Computational method

3.3.2 Error measurements

3.3.3 Numerical evaluation results

4. Conclusion

Funding

Acknowledgments

Disclosures

References and links

Cited By

Figures (9)

Equations (67)

Biomedical Optics Express