An F_N method for the radiative transport equation in three dimensions

We consider light propagating in a homogeneous random medium occupying the half-space ${{\mathbb{R}}}_{+}^{3}$ ($=\{{\bf{r}}\in {{\mathbb{R}}}^{3};{\bf{r}}=({\boldsymbol{\rho }},z),{\boldsymbol{\rho }}\in {{\mathbb{R}}}^{2},z\gt 0\}$) with the boundary at z = 0. The specific intensity $I({\bf{r}},\hat{{\bf{s}}})$ (${\bf{r}}\in {{\mathbb{R}}}_{+}^{3}$, $\hat{{\bf{s}}}\in {{\mathbb{S}}}^{2}$) of light obeys the following radiative transport equation.

$$\begin{eqnarray}\left\{\begin{array}{ll}\hat{{\bf{s}}}\cdot \nabla I({\bf{r}},\hat{{\bf{s}}})+I({\bf{r}},\hat{{\bf{s}}})=\varpi {\displaystyle \int }_{{{\mathbb{S}}}^{2}}p(\hat{{\bf{s}}},\hat{{\bf{s}}}\prime )I({\bf{r}},\hat{{\bf{s}}}\prime )\;{\rm{d}}\hat{{\bf{s}}}\prime +S({\bf{r}},\hat{{\bf{s}}}), & z\gt 0,\\ I({\bf{r}},\hat{{\bf{s}}})=f({\boldsymbol{\rho }},\hat{{\bf{s}}}), & z=0,\mu \in (0,1],\\ I({\bf{r}},\hat{{\bf{s}}})\to 0, & z\to \infty ,\end{array}\right.\end{eqnarray} \tag{ 1 }$$

where $f({\boldsymbol{\rho }},\hat{{\bf{s}}})$ is the incident beam and $S({\bf{r}},\hat{{\bf{s}}})$ is the internal source. Let μ and $\varphi$ be the cosine of the polar angle and the azimuthal angle of $\hat{{\bf{s}}}\in {{\mathbb{S}}}^{2}$. Here $\varpi \in (0,1)$ is the albedo for single scattering. Using the absorption and scattering parameters ${\mu }_{a}$ and ${\mu }_{s}$, we have $\varpi ={\mu }_{s}/{\mu }_{t}$, where ${\mu }_{t}={\mu }_{a}+{\mu }_{s}$ is the total attenuation. The above form (1) implies that ${\bf{r}}$ is normalized by ${\mu }_{t}$. Furthermore $p(\hat{{\bf{s}}},\hat{{\bf{s}}}\prime )$ is the scattering phase function which is normalized as

$$\begin{eqnarray*}&&{\displaystyle \int }_{{{\mathbb{S}}}^{2}}p(\hat{{\bf{s}}}\prime ,\hat{{\bf{s}}})\;{\rm{d}}\hat{{\bf{s}}}\prime =1,\quad \hat{{\bf{s}}}\in {{\mathbb{S}}}^{2}.\end{eqnarray*}$$

The radiative transport equation or the linear Boltzmann equation governs transport processes of noninteracting particles such as neutrons in a reactor as well as light propagation in random media such as fog, clouds, and biological tissue.

In this paper we will present a numerical method of solving (1) by extending the F_N method (F stands for facile) to three dimensions. The F_N method first developed by Siewert [51] is a method of obtaining the specific intensity in one dimension making use of orthogonality relations of singular eigenfunctions [4, 6, 10]. The use of rotated reference frames [43, 48, 50] makes it possible to extend the F_N method to three dimensions.

In 1960 Case considered the time-independent one-dimensional radiative transport equation with isotropic scattering and solved the equation with separation of variables by finding singular eigenfunctions [4]. The method was soon extended to the case of anisotropic scattering without [44, 47] and with [45] azimuthal dependence. Such singular-eigenfunction approach is sometimes called Caseology. In this method, solutions to the one-dimensional radiative transport equation are given by a superposition of singular eigenfunctions. The existence and uniqueness of such solutions were proved [25–28]. In the F_N method, there is no need of evaluating singular functions although the fact that the specific intensity consists of singular eigenfunctions is used. In one dimension, the radiative transport equation was solved by the F_N method in the slab geometry for isotropic scattering [12, 52] and anisotropic scattering without [9, 16, 51] and with [19, 20] azimuthal dependence. The method was also extended to multigroup [14]. After finding the specific intensity on the boundary, we can further calculate the specific intensity inside the medium [16]. The uniqueness of the solution to the key F_N equation was proved [29]. For isotropic scattering, the three dimensional radiative transport equation was solved with the F_N method [11, 53] using the pseudo-problem [55], which is based on plane-wave decomposition. See the review article by Garcia [13].

In 1964 Dede used rotated reference frames to solve the three-dimensional radiative transport equation with the P_N method [8]. Dede pointed out that equations in three dimensions reduce to one-dimensional equations if reference frames are rotated in the direction of the Fourier vector. Kobayashi developed Dedeʼs calculation and computed coefficients in the P_N expansion by solving a three-term recurrence relation recursively starting with the initial term [24]. In 2004 Markel obtained the coefficients in terms of eigenvalues and eigenvectors of the tridiagonal matrix originating from the three-term recurrence relation, and showed that the specific intensity can be efficiently computed [43]. With the use of eigenvalues, the relation to Caseʼs method became visible. This new formulation can be viewed as separation of variables in which the eigenvalues are separation constants [50]. Moreover it was found that any complex unit vector can be used to rotate reference frames [48]. This generalization makes it possible to solve boundary value problems in the form of plane-wave decomposition [41]. It was then found that the structure of separation of variables implies Caseʼs method in rotated reference frames [42]. Thus the singular-eigenfunction approach was extended to three dimensions. Indeed the method of rotated reference frames is a three-dimensional extension of the spherical-harmonic expansion [1, 49] in Caseology.

The usefulness of the method of rotated reference frames has been numerically justified for a two-dimensional rectangular domain [24], a three-dimensional infinite medium [43, 48], the slab geometry in three dimensions [41], in flatland [30, 31, 38], in the half-space geometry [33, 35–37, 39], and the time-dependent equation in an infinite medium [32, 34]. The method was also used to experimentally determine optical properties of turbid media [56, 57]. It is expected that more accurate numerical values are obtained if higher terms in the series are taken into account. Although the method of rotated reference frames is an efficient method, the obtained values become unstable when high-degree spherical harmonics are used. The three-dimensional F_N method developed in the present paper does not suffer from this instability.

By assuming that scatterers are spherically symmetric, we model $p(\hat{{\bf{s}}},\hat{{\bf{s}}}\prime )$ as

$$\begin{eqnarray}&&p(\hat{{\bf{s}}},\hat{{\bf{s}}}\prime )=\displaystyle \frac{1}{4\pi }\displaystyle \sum _{l=0}^{L}{\beta }_{l}{P}_{l}(\hat{{\bf{s}}}\cdot \hat{{\bf{s}}}\prime )=\displaystyle \sum _{l=0}^{L}\displaystyle \sum _{m=-l}^{l}\displaystyle \frac{{\beta }_{l}}{2l+1}{Y}_{{lm}}(\hat{{\bf{s}}}){Y}_{{lm}}^{*}(\hat{{\bf{s}}}\prime ),\end{eqnarray} \tag{ 2 }$$

where $L\geqslant 1$, and ${\beta }_{0}=1$, $0\lt {\beta }_{l}\lt 2l+1$ for $l\geqslant 1$. Moreover P_l are Legendre polynomials and Y_lm are spherical harmonics. We introduce the scattering asymmetry parameter ${\rm{g}}$ as ${\beta }_{l}=(2l+1){{\rm{g}}}^{l}$ ($0\lt {\rm{g}}\lt 1$). The Henyey–Greenstein model [22] is obtained in the limit $L\to \infty$.

Let us define

$$\begin{eqnarray*}&&\tilde{I}({\bf{q}},z,\hat{{\bf{s}}})={\displaystyle \int }_{{{\mathbb{R}}}^{2}}{{\rm{e}}}^{{\rm{i}}{\bf{q}}\cdot {\boldsymbol{\rho }}}I({\bf{r}},\hat{{\bf{s}}})\;d{\boldsymbol{\rho }},\quad {\bf{q}}\in {{\mathbb{R}}}^{2}.\end{eqnarray*}$$

We similarly define $\tilde{f}({\bf{q}},\hat{{\bf{s}}})$ and $\tilde{S}({\bf{q}},z,\hat{{\bf{s}}})$. Let us express the upper and lower hemispheres as ${{\mathbb{S}}}_{\pm }^{2}=\{\hat{{\bf{s}}}\in {{\mathbb{S}}}^{2};\pm \mu \gt 0\}$. We expand the Fourier transform of the reflected light $\tilde{I}({\bf{q}},0,-\hat{{\bf{s}}})$ ($\hat{{\bf{s}}}\in {{\mathbb{S}}}_{+}^{2}$) as

$$\begin{eqnarray}&&\tilde{I}({\bf{q}},0,-\hat{{\bf{s}}})\approx \displaystyle \sum _{m=-{l}_{\mathrm{max}}}^{{l}_{\mathrm{max}}}\displaystyle \sum _{\alpha =0}^{\left\lfloor ({l}_{\mathrm{max}}-| m| )/2\right\rfloor}{c}_{| m| +2\alpha ,m}({\bf{q}}){Y}_{| m| +2\alpha ,m}(\hat{{\bf{s}}}),\end{eqnarray} \tag{ 3 }$$

where ${l}_{\mathrm{max}}$ is the highest degree of the expansion (${l}_{\mathrm{max}}\geqslant L$). Only same-parity degrees are taken because the three-term recurrence relation of associated Legendre polynomials implies that Y_lm of opposite-parity l are not independent [20, 48]. This expansion in (3) can be compared to the P_N method [6], but the F_N method is more efficient because the spatial dependence of the specific intensity is analytically given and the orthogonality relation among three-dimensional singular eigenfunctions can be used (see section 2.3). On the other hand, $\tilde{I}({\bf{q}},0,-\hat{{\bf{s}}})$ is given as a linear combination of eigenmodes ${{\mathcal{R}}}_{\hat{{\bf{k}}}(\nu ,{\bf{q}})}{\Phi }_{\nu }^{{m}^{\prime }}(\hat{{\bf{s}}})$ [42], for which notations are introduced in section 2.3. They satisfy orthogonality relations. For simplicity let us assume $S({\bf{r}},\hat{{\bf{s}}})=0$. Making use of the fact that $\tilde{I}$ contains only decaying modes, we have (See (34) for the general case)

$$\begin{eqnarray}&&{\displaystyle \int }_{{{\mathbb{S}}}^{2}}\mu \left({{\mathcal{R}}}_{\hat{{\bf{k}}}(-\xi ,{\bf{q}})}{\Phi }_{-\xi }^{m\prime *}(\hat{{\bf{s}}})\right)\tilde{I}({\bf{q}},0,\hat{{\bf{s}}})\;{\rm{d}}\hat{{\bf{s}}}=0,\qquad \xi \gt 0.\end{eqnarray} \tag{ 4 }$$

The above equation results in a linear system for ${c}_{| m| +2\alpha ,m}({\bf{q}})$. The specific intensity of the reflected light is then calculated as

$$\begin{eqnarray*}&&I({\boldsymbol{\rho }},0,-\hat{{\bf{s}}})\approx \displaystyle \frac{1}{{(2\pi )}^{2}}{\displaystyle \int }_{{{\mathbb{R}}}^{2}}{{\rm{e}}}^{-{\rm{i}}{\bf{q}}\cdot {\boldsymbol{\rho }}}\displaystyle \sum _{m=-{l}_{\mathrm{max}}}^{{l}_{\mathrm{max}}}\displaystyle \sum _{l=| m| ,| m| +2,...}{c}_{{lm}}({\bf{q}}){Y}_{{lm}}(\hat{{\bf{s}}})\;{\rm{d}}{\bf{q}},\end{eqnarray*}$$

where $\mu \in (0,1]$.

Remark 1.1. Isotropic scattering ${\rm{g}}=0$ is possible. However we need to change the collocation scheme for obtaining c_lm. For the sake of simplicity, we assume ${\rm{g}}\gt 0$ in this paper.

Remark 1.2. The expansion in (3) can be compared to the method of rotated reference frames, which expands every eigenmode with spherical harmonics:

$$\begin{eqnarray}&&{{\mathcal{R}}}_{\hat{{\bf{k}}}(\nu ,{\bf{q}})}{\Phi }_{\nu }^{{m}^{\prime }}(\hat{{\bf{s}}})\approx \displaystyle \sum _{l=0}^{{l}_{\mathrm{max}}}\displaystyle \sum _{m=-l}^{l}{c}_{{lm}}^{{m}^{\prime }}(\nu ){{\mathcal{R}}}_{\hat{{\bf{k}}}(\nu ,{\bf{q}})}{Y}_{{lm}}(\hat{{\bf{s}}}),\end{eqnarray} \tag{ 5 }$$

with some coefficients ${c}_{{lm}}^{{m}^{\prime }}(\nu )$. This causes numerical instability regardless of $f({\boldsymbol{\rho }},\hat{{\bf{s}}})$ and $S({\bf{r}},\hat{{\bf{s}}})$ when ${l}_{\mathrm{max}}$ is increased to achieve higher precision. For example, let us consider a simple case of L = 0, $m\prime =0$, $\mathrm{cos}(\varphi -{\varphi }_{{\bf{q}}})=0$, and $\nu \ne \mu {\hat{k}}_{z}(\nu q)$. Noting that ${{\mathcal{R}}}_{\hat{{\bf{k}}}(\nu ,{\bf{q}})}{Y}_{{lm}}(\hat{{\bf{s}}})=\frac{2l+1}{4\pi }{P}_{l}^{m}\left(\mu {\hat{k}}_{z}(\nu q)\right){{\mathcal{R}}}_{\hat{{\bf{k}}}(\nu ,{\bf{q}})}{{\rm{e}}}^{{\rm{i}}m\varphi }$, we see that the right-hand side of (5) is a polynomial of $\mu {\hat{k}}_{z}(\nu q)$. On the left-hand side, we have ${{\mathcal{R}}}_{\hat{{\bf{k}}}(\nu ,{\bf{q}})}{\Phi }_{\nu }^{{m}^{\prime }}(\hat{{\bf{s}}})=\frac{\varpi \nu }{2}{\left[\nu -\mu {\hat{k}}_{z}(\nu q)\right]}^{-1}=\frac{\varpi }{2}\left[1+(1/\nu )\mu {\hat{k}}_{z}(\nu q)+{(1/\nu )}^{2}{\left(\mu {\hat{k}}_{z}(\nu q)\right)}^{2}+\cdots \right]$. This series is divergent if $\nu -\mu {\hat{k}}_{z}(\nu q)\lt 0$. In general, the instability takes place due to the same mechanism. Figure 1 shows the exiting current on the boundary (z = 0) as a function of ${l}_{\mathrm{max}}$. See section 4 for the details.

The remainder of the paper is organized as follows. In section 2 we introduce singular eigenfunctions and rotated reference frames. In section 3 we consider the F_N method in three dimensions. The key F_N equation is obtained in (35), from which the coefficients c_lm in (3) are computed. In section 4 the three-dimensional F_N method is numerically tested for structured illumination. Section 5 is devoted to concluding remarks. Finally structured illumination by the method of rotated reference frames is summarized in appendix A.

Zoom In Zoom Out Reset image size

To develop the F_N method in three dimensions in section 3, we give brief reviews and define our notations in this section. In section 2.1, we introduce polynomials ${g}_{l}^{m}$ and ${p}_{l}^{m}$. In section 2.2, Caseʼs singular-eigenfunction approach is explained. In section 2.3, we give a review on singular eigenfunctions in three dimensions. In section 2.4, it is sketched how the method of rotated reference frames is obtained using three-dimensional singular eigenfunctions.

2.1. Polynomials

Definition 2.1. We introduce h_l ($l=0,1,...$) as

$$\begin{eqnarray*}&&{h}_{l}=2l+1-\varpi {\beta }_{l}{\chi }_{[0,L]}(l),\end{eqnarray*}$$

with ${\chi }_{[0,L]}(l)$ the step function ($\chi =1$ for $0\leqslant l\leqslant L$ and $\chi =0$ otherwise).

([17, 18]).

Definition 2.2 The normalized Chandrasekhar polynomials ${g}_{l}^{m}(\xi )$ ($m\geqslant 0$, $l\geqslant m$, $\nu \in {\mathbb{R}}$) are given by the three-term recurrence relation

$$\begin{eqnarray}&&\nu {h}_{l}{g}_{l}^{m}(\nu )=\sqrt{{(l+1)}^{2}-{m}^{2}}{g}_{l+1}^{m}(\nu )+\sqrt{{l}^{2}-{m}^{2}}{g}_{l-1}^{m}(\nu ),\end{eqnarray} \tag{ 6 }$$

with the initial term

$$\begin{eqnarray}&&{g}_{m}^{m}(\nu )=\displaystyle \frac{(2m-1)!!}{\sqrt{(2m)!}}=\displaystyle \frac{\sqrt{(2m)!}}{{2}^{m}m!}.\end{eqnarray} \tag{ 7 }$$

We note that

$$\begin{eqnarray*}&&{g}_{l}^{-m}(\nu )={(-1)}^{m}{g}_{l}^{m}(\nu ),\qquad {g}_{l}^{m}(-\nu )={(-1)}^{l+m}{g}_{l}^{m}(\nu ).\end{eqnarray*}$$

The polynomials ${g}_{l}^{m}$ are obtained if we multiply Chandrasekhar polynomials [7] by $\sqrt{(l-m)!/(l+m)!}$ [54].

Definition 2.3. The polynomials ${p}_{l}^{m}(\mu )$ ($m\geqslant 0$, $l\geqslant m$) are introduced as

$$\begin{eqnarray}&&{p}_{l}^{m}(\mu )={(-1)}^{m}\sqrt{\displaystyle \frac{(l-m)!}{(l+m)!}}{P}_{l}^{m}(\mu ){(1-{\mu }^{2})}^{-m/2}=\sqrt{\displaystyle \frac{(l-m)!}{(l+m)!}}\displaystyle \frac{{d}^{m}}{d{\mu }^{m}}{P}_{l}(\mu ),\end{eqnarray} \tag{ 8 }$$

where ${P}_{l}(\mu )$ is the Legendre polynomial of degree l and ${P}_{l}^{m}(\mu )$ is the associated Legendre polynomial of degree l and order m.

We have

$$\begin{eqnarray*}&&{p}_{l}^{-m}(\mu )={(-1)}^{m}{p}_{l}^{m}(\mu ).\end{eqnarray*}$$

The polynomials satisfy the three-term recurrence relation

$$\begin{eqnarray}&&\sqrt{{l}^{2}-{m}^{2}}{p}_{l-1}^{m}(\mu )-(2l+1)\mu {p}_{l}^{m}(\mu )+\sqrt{{(l+1)}^{2}-{m}^{2}}{p}_{l+1}^{m}(\mu )=0,\end{eqnarray} \tag{ 9 }$$

with

$$\begin{eqnarray*}&&{p}_{| m| }^{| m| }(\mu )=\displaystyle \frac{(2| m| -1)!!}{\sqrt{(2| m| )!}}=\displaystyle \frac{\sqrt{(2| m| )!}}{{2}^{| m| }| m| !},\end{eqnarray*}$$

and the orthogonality relation

$$\begin{eqnarray*}&&{\displaystyle \int }_{-1}^{1}{p}_{l}^{m}(\mu ){p}_{{l}^{\prime }}^{m}(\mu ){\left(1-{\mu }^{2}\right)}^{| m| }\;{\rm{d}}\mu =\displaystyle \frac{2}{2l+1}{\delta }_{{{ll}}^{\prime }}.\end{eqnarray*}$$

2.2. Singular eigenfunctions for one dimension

We will first investigate the one-dimensional homogeneous radiative transport equation (10) and then consider the three dimensional equation (26). Let us begin with

$$\begin{eqnarray}&&\mu \displaystyle \frac{\partial }{\partial z}I(z,\hat{{\bf{s}}})+I(z,\hat{{\bf{s}}})=\varpi {\displaystyle \int }_{{{\mathbb{S}}}^{2}}p(\hat{{\bf{s}}},\hat{{\bf{s}}}\prime )I(z,\hat{{\bf{s}}}\prime )\;{\rm{d}}\hat{{\bf{s}}}\prime ,\end{eqnarray} \tag{ 10 }$$

where $z\in {\mathbb{R}}$, $\hat{{\bf{s}}}\in {{\mathbb{S}}}^{2}$ and $p(\hat{{\bf{s}}},\hat{{\bf{s}}}\prime )$ is given in (2). Separated solutions to (10) are given by [4, 45, 47]

$$\begin{eqnarray}&&I(z,\hat{{\bf{s}}})={\Phi }_{\nu }^{m}(\hat{{\bf{s}}}){{\rm{e}}}^{-z/\nu },\end{eqnarray} \tag{ 11 }$$

where $\nu \in {\mathbb{R}}$ is a separation constant, m ($| m| \leqslant L$) is an integer, and

$$\begin{eqnarray}&&{\Phi }_{\nu }^{m}(\hat{{\bf{s}}})={\phi }^{m}(\nu ,\mu ){\left(1-{\mu }^{2}\right)}^{| m| /2}{{\rm{e}}}^{{\rm{i}}m\varphi }.\end{eqnarray} \tag{ 12 }$$

Here ${\phi }^{m}(\nu ,\mu )$ satisfies

$$\begin{eqnarray*}&&{\displaystyle \int }_{-1}^{1}{\phi }^{m}(\nu ,\mu ){\left(1-{\mu }^{2}\right)}^{| m| }\;{\rm{d}}\mu =1.\end{eqnarray*}$$

By plugging (11) into (10) we obtain

$$\begin{eqnarray}&&\left(1-\displaystyle \frac{\mu }{\nu }\right){\Phi }_{\nu }^{m}(\hat{{\bf{s}}})=\varpi {\displaystyle \int }_{{{\mathbb{S}}}^{2}}p(\hat{{\bf{s}}},\hat{{\bf{s}}}\prime ){\Phi }_{\nu }^{m}(\hat{{\bf{s}}}\prime )\;{\rm{d}}\hat{{\bf{s}}}\prime .\end{eqnarray} \tag{ 13 }$$

We multiply (13) by ${Y}_{l\prime m\prime }^{*}(\hat{{\bf{s}}})$ and integrate both sides over ${{\mathbb{S}}}^{2}$. By noticing the expression of $p(\hat{{\bf{s}}},\hat{{\bf{s}}}\prime )$ in (2) and rearranging terms, we obtain

$$\begin{eqnarray}&&\nu \left(1-\displaystyle \frac{\varpi {\beta }_{{l}^{\prime }}}{2l\prime +1}{\chi }_{[0,L]}(l\prime )\right){\displaystyle \int }_{{{\mathbb{S}}}^{2}}{Y}_{l\prime m\prime }^{*}(\hat{{\bf{s}}}){\Phi }_{\nu }^{m}(\hat{{\bf{s}}})\;{\rm{d}}\hat{{\bf{s}}}={\displaystyle \int }_{{{\mathbb{S}}}^{2}}\mu {Y}_{l\prime m\prime }^{*}(\hat{{\bf{s}}}){\Phi }_{\nu }^{m}(\hat{{\bf{s}}})\;{\rm{d}}\hat{{\bf{s}}}.\end{eqnarray} \tag{ 14 }$$

Using the recurrence relation $(2l+1)\mu {P}_{l}^{m}(\mu )=(l+1-m){P}_{l+1}^{m}(\mu )+(l+m){P}_{l-1}^{m}(\mu )$, we see that (14) becomes the three-term recurrence relation (6) for $m\prime =m$. That is, we obtain

$$\begin{eqnarray}&&{g}_{l}^{m}(\nu )={(-1)}^{m}\sqrt{\displaystyle \frac{(l-m)!}{(l+m)!}}{\displaystyle \int }_{-1}^{1}{\phi }^{m}(\nu ,\mu ){(1-{\mu }^{2})}^{| m| /2}{P}_{l}^{m}(\mu )\;{\rm{d}}\mu .\end{eqnarray} \tag{ 15 }$$

Noting that ${P}_{m}^{m}(\mu )={(-1)}^{m}(2m-1)!!{(1-{\mu }^{2})}^{m/2}$ ($m\geqslant 0$), we see that (15) satisfies (7).

Let us rewrite (13) as

$$\begin{eqnarray*}&&\left(1-\displaystyle \frac{\mu }{\nu }\right){\phi }^{m}(\nu ,\mu )=\displaystyle \frac{\varpi }{2}\displaystyle \sum _{l\prime =| m| }^{L}{\beta }_{{l}^{\prime }}{p}_{{l}^{\prime }}^{m}(\mu ){g}_{{l}^{\prime }}^{m}(\nu ).\end{eqnarray*}$$

We define g^m as

$$\begin{eqnarray}&&{g}^{m}(\nu ,\mu )=\displaystyle \sum _{l=| m| }^{L}{\beta }_{l}{p}_{l}^{m}(\mu ){g}_{l}^{m}(\nu ).\end{eqnarray} \tag{ 16 }$$

Singular eigenfunctions ${\phi }^{m}(\nu ,\mu )$ are thus obtained as

$$\begin{eqnarray*}&&{\phi }^{m}(\nu ,\mu )=\displaystyle \frac{\varpi \nu }{2}{\mathcal{P}}\displaystyle \frac{{g}^{m}(\nu ,\mu )}{\nu -\mu }+{\lambda }^{m}(\nu ){\left(1-{\mu }^{2}\right)}^{-| m| }\delta (\nu -\mu ),\end{eqnarray*}$$

where ${\mathcal{P}}$ denotes the Cauchy principal value. Here the separation constant ν has discrete values $\pm {\nu }_{j}^{m}$ (${\nu }_{j}^{m}\gt 1$, $j=0,1,...,{M}^{m}-1$) and the continuous spectrum between −1 and 1. The number M^m of discrete eigenvalues depends on ϖ and ${\beta }_{l}$. The function ${\lambda }^{m}(\nu )$ is given by

$$\begin{eqnarray*}&&{\lambda }^{m}(\nu )=1-\displaystyle \frac{\varpi \nu }{2}\;{\mathcal{P}}{\int }_{-1}^{1}\displaystyle \frac{{g}^{m}(\nu ,\mu )}{\nu -\mu }{(1-{\mu }^{2})}^{| m| }\;{\rm{d}}\mu .\end{eqnarray*}$$

Discrete eigenvalues satisfy

$$\begin{eqnarray*}&&{\Lambda }^{m}({\nu }_{j}^{m})=0,\end{eqnarray*}$$

where for $w\in {\mathbb{C}}$

$$\begin{eqnarray*}&&{\Lambda }^{m}(w)=1-\displaystyle \frac{\varpi w}{2}{\displaystyle \int }_{-1}^{1}\displaystyle \frac{{g}^{m}(w,\mu )}{w-\mu }{(1-{\mu }^{2})}^{| m| }\;{\rm{d}}\mu .\end{eqnarray*}$$

By using ${P}_{l}^{-m}={P}_{l}^{m}{(-1)}^{m}(l-m)!/(l+m)!$, we can readily check that ${g}_{l}^{m}(\nu )$ in (15) satisfy ${g}_{l}^{-m}(\nu )={(-1)}^{m}{g}_{l}^{m}(\nu )$. This implies ${\phi }^{-m}={\phi }^{m}$. Singular eigenfunctions ${\phi }^{m}(\nu ,\mu )$ satisfy [4, 45, 47]

$$\begin{eqnarray*}&&{\displaystyle \int }_{-1}^{1}\mu {\phi }^{m}(\nu ,\mu ){\phi }^{m}({\nu }^{\prime },\mu )\;{\rm{d}}\mu ={{\mathcal{N}}}^{m}(\nu ){\delta }_{\nu \nu \prime },\end{eqnarray*}$$

where the Kronecker delta ${\delta }_{\nu \nu \prime }$ is replaced by the Dirac delta $\delta (\nu -{\nu }^{\prime })$ if $\nu ,\nu \prime$ are in the continuous spectrum. The normalization factor ${{\mathcal{N}}}^{m}(\nu )$ is given by

$$\begin{eqnarray}{{\mathcal{N}}}^{m}(\nu )=\left\{\begin{array}{ll}\displaystyle \frac{1}{2}{({\nu }_{j}^{m})}^{2}g({\nu }_{j}^{m},{\nu }_{j}^{m}){\left.\displaystyle \frac{{\rm{d}}{\Lambda }^{m}(w)}{{\rm{d}}w}\right|}_{w={\nu }_{j}^{m}}, & \nu ={\nu }_{j}^{m},\\ \nu {\Lambda }^{m+}(\nu ){\Lambda }^{m-}(\nu ){(1-{\nu }^{2})}^{-| m| }, & \nu \in (-1,1),\end{array}\right.\end{eqnarray} \tag{ 17 }$$

where ${\Lambda }^{m\pm }(\nu )={\mathrm{lim}}_{\epsilon \to {0}^{+}}\;{\Lambda }^{m}(\nu \pm {\rm{i}}\epsilon )$.

We can numerically obtain the discrete eigenvalues ${\nu }_{j}^{m}$ as eigenvalues of a tridiagonal matrix $B(m)$ below. For l_B ($\geqslant L$) and m ($-L\leqslant m\leqslant L$), the matrix $B(m)$ is given by

$$\begin{eqnarray}B(m)=\left(\begin{array}{ccccc}0 & {b}_{| m| +1} & 0 & & \\ {b}_{| m| +1} & 0 & {b}_{| m| +2} & & \\ 0 & {b}_{| m| +2} & 0 & \ddots & \\ & & \ddots & \ddots & {b}_{{l}_{B}}\\ & & & {b}_{{l}_{B}} & 0\end{array}\right),\end{eqnarray} \tag{ 18 }$$

where ${b}_{l}(m)=\sqrt{({l}^{2}-{m}^{2})/({h}_{l}{h}_{l-1})}$. The matrix $B(m)$ has $({l}_{B}-| m| +1)/2$ or $({l}_{B}-| m| )/2$ positive eigenvalues for ${l}_{B}-| m| +1$ even or odd, respectively. To see how $B(m)$ is obtained, we first prove the following proposition.

([15]).

Proposition 2.4 Discrete eigenvalues are zeros of ${g}_{l}^{m}$ as $l\to \infty$.

Proof. We define

$$\begin{eqnarray*}&&{q}_{l}^{m}(w)=\displaystyle \frac{1}{2}{\displaystyle \int }_{-1}^{1}\displaystyle \frac{{p}_{l}^{m}(\mu )}{w-\mu }{(1-{\mu }^{2})}^{| m| }\;{\rm{d}}\mu ,\quad w\notin [-1,1].\end{eqnarray*}$$

For $\nu \notin [-1,1]$, the three-term recurrence relation of ${p}_{l}^{m}$ implies

$$\begin{eqnarray}\begin{array}{rcl}\sqrt{{(l+1)}^{2}-{m}^{2}}{q}_{l+1}^{m}(\nu ) & = & (2l+1)\nu {q}_{l}^{m}(\nu )-\sqrt{{l}^{2}-{m}^{2}}{q}_{l-1}^{m}(\nu )\\ & & -{(\mathrm{sgn}(m))}^{m}\displaystyle \frac{\sqrt{(2| m| )!}}{(2| m| -1)!!}{\delta }_{l,| m| }.\end{array}\end{eqnarray} \tag{ 19 }$$

By subtracting (6) multiplied by ${q}_{l}^{m}(\nu )$ on both sides from (19) multiplied by ${g}_{l}^{m}(\nu )$ on both sides, we obtain

$$\begin{eqnarray*}&&\begin{array}{l}\sqrt{{(l+1)}^{2}-{m}^{2}}\left({q}_{l+1}^{m}(\nu ){g}_{l}^{m}(\nu )-{q}_{l}^{m}(\nu ){g}_{l+1}^{m}(\nu )\right)\\ \qquad =(2l+1)\nu {q}_{l}^{m}(\nu ){g}_{l}^{m}(\nu )-\nu {h}_{l}{q}_{l}^{m}(\nu ){g}_{l}^{m}(\nu )\\ \qquad \;\;\;-\sqrt{{l}^{2}-{m}^{2}}\left({q}_{l-1}^{m}(\nu ){g}_{l}^{m}(\nu )-{q}_{l}^{m}(\nu ){g}_{l-1}^{m}(\nu )\right)-{\delta }_{l,| m| }.\end{array}\end{eqnarray*}$$

Suppose ${l}_{B}\geqslant L$. By taking the summation ${\displaystyle \sum }_{l=| m| }^{{l}_{B}}$ we obtain

$$\begin{eqnarray*}&&\begin{array}{l}\sqrt{{({l}_{B}+1)}^{2}-{m}^{2}}\left[{q}_{{l}_{B}+1}^{m}(\nu ){g}_{{l}_{B}}^{m}(\nu )-{q}_{{l}_{B}}^{m}(\nu ){g}_{{l}_{B}+1}^{m}(\nu )\right]\\ \qquad =\displaystyle \sum _{l=| m| }^{{l}_{B}}\left((2l+1)\nu -\nu {h}_{l}\right){q}_{l}^{m}(\nu ){g}_{l}^{m}(\nu )-1.\end{array}\end{eqnarray*}$$

Noting that ${\Lambda }^{m}(\nu )=1-\varpi \nu {\displaystyle \sum }_{l=| m| }^{L}{\beta }_{l}{g}_{l}^{m}(\nu ){q}_{l}^{m}(\nu )$, we obtain (the Christoffel–Darboux formula)

$$\begin{eqnarray}&&{\Lambda }^{m}(\nu )=\sqrt{{({l}_{B}+1)}^{2}-{m}^{2}}\left[{q}_{{l}_{B}}^{m}(\nu ){g}_{{l}_{B}+1}^{m}(\nu )-{q}_{{l}_{B}+1}^{m}(\nu ){g}_{{l}_{B}}^{m}(\nu )\right].\end{eqnarray} \tag{ 20 }$$

Next we subtract (19) multiplied by ${p}_{l}^{m}(\nu )$ on both sides from (9) multiplied by ${q}_{l}^{m}(\nu )$ on both sides. By summing the resulting expression over l from $| m|$ to l_B, we have

$$\begin{eqnarray}&&1=\sqrt{{({l}_{B}+1)}^{2}-{m}^{2}}\left({p}_{{l}_{B}+1}^{m}(\mu ){q}_{{l}_{B}}^{m}(\nu )-{p}_{{l}_{B}}^{m}(\nu ){q}_{{l}_{B}+1}^{m}(\nu )\right).\end{eqnarray} \tag{ 21 }$$

Similarly we subtract (6) multiplied by ${p}_{l}^{m}(\nu )$ on both sides from (9) multiplied by ${g}_{l}^{m}(\nu )$ on both sides, and take the sum over l from $| m|$ to l_B. We obtain

$$\begin{eqnarray}&&\varpi \nu {g}^{m}(\nu ,\nu )=\sqrt{{({l}_{B}+1)}^{2}-{m}^{2}}\left({p}_{{l}_{B}+1}^{m}(\mu ){g}_{{l}_{B}}^{m}(\nu )-{p}_{{l}_{B}}^{m}(\nu ){g}_{{l}_{B}+1}^{m}(\nu )\right).\end{eqnarray} \tag{ 22 }$$

Using (20), (21), and (22), we obtain

$$\begin{eqnarray*}\begin{array}{rcl}{p}_{{l}_{B}+1}^{m}(\nu ){\Lambda }^{m}(\nu ) & = & \sqrt{{({l}_{B}+1)}^{2}-{m}^{2}}\left[{p}_{{l}_{B}+1}^{m}(\nu ){q}_{{l}_{B}}^{m}(\nu ){g}_{{l}_{B}+1}^{m}(\nu )-{p}_{{l}_{B}+1}^{m}(\nu ){q}_{{l}_{B}+1}^{m}(\nu ){g}_{{l}_{B}}^{m}(\nu )\right]\\ & = & {g}_{{l}_{B}+1}^{m}(\nu )+\sqrt{{({l}_{B}+1)}^{2}-{m}^{2}}\left[{p}_{{l}_{B}}^{m}(\nu ){g}_{{l}_{B}+1}^{m}(\nu )-{p}_{{l}_{B}+1}^{m}(\nu ){g}_{{l}_{B}}^{m}(\nu )\right]{q}_{{l}_{B}+1}^{m}(\nu )\\ & = & {g}_{{l}_{B}+1}^{m}(\nu )-\varpi \nu {g}^{m}(\nu ,\nu ){q}_{{l}_{B}+1}^{m}(\nu ).\end{array}\end{eqnarray*}$$

We note that ${\mathrm{lim}}_{l\to \infty }\;{q}_{l}^{m}(w)/{p}_{l}^{m}(w)={\mathrm{lim}}_{l\to \infty }\;{Q}_{l}^{m}(w)/{P}_{l}^{m}(w)=0$ ($w\notin [-1,1]$), where ${Q}_{l}^{m}$ is the associated Legendre polynomial of the second kind. Therefore we obtain

$$\begin{eqnarray*}&&{\Lambda }^{m}(\nu )=\underset{{l}_{B}\to \infty }{\mathrm{lim}}\;\displaystyle \frac{{g}_{{l}_{B}+1}^{m}(\nu )}{{p}_{{l}_{B}+1}^{m}(\nu )}.\end{eqnarray*}$$

Thus the proof is completed. □

Let us recall that the recurrence relation (6) for ${g}_{l}^{m}(\nu )$ is derived for an eigenvalue ν in (11) and rewrite (6) as

$$\begin{eqnarray*}&&\sqrt{\displaystyle \frac{{l}^{2}-{m}^{2}}{{h}_{l}{h}_{l-1}}}\sqrt{{h}_{l-1}}{g}_{l-1}^{m}(\nu )+\sqrt{\displaystyle \frac{{(l+1)}^{2}-{m}^{2}}{{h}_{l}{h}_{l+1}}}\sqrt{{h}_{l+1}}{g}_{l+1}^{m}(\nu )=\nu \sqrt{{h}_{l}}{g}_{l}^{m}(\nu ).\end{eqnarray*}$$

Hence eigenvalues of $B(m)$ are zeros of ${g}_{{l}_{B}+1}^{m}$. Together with proposition 2.4, we see that discrete eigenvalues ${\nu }_{j}^{m}$ can be computed as eigenvalues of $B(m)$ for sufficiently large l_B. More sophisticated ways of obtaining discrete eigenvalues are discussed in [17].

The tridiagonal matrix $B(m)$ can be alternatively obtained as follows. Let us write ${\Phi }_{\nu }^{m}(\hat{{\bf{s}}})$ as

$$\begin{eqnarray*}&&{\Phi }_{\nu }^{m}(\hat{{\bf{s}}})=\displaystyle \sum _{l\prime =0}^{\infty }\displaystyle \sum _{m\prime =-l\prime }^{{l}^{\prime }}{c}_{l\prime m\prime }^{m}(\nu ){Y}_{l\prime m\prime }(\hat{{\bf{s}}}),\end{eqnarray*}$$

where ${c}_{l\prime m\prime }^{m}(\nu )\in {\mathbb{C}}$. Then (13) can be rewritten as

$$\begin{eqnarray}&&\left(1-\displaystyle \frac{\mu }{\nu }\right)\displaystyle \sum _{l\prime m\prime }{c}_{l\prime m\prime }{Y}_{l\prime m\prime }(\hat{{\bf{s}}})=\varpi \displaystyle \sum _{l\prime m\prime }\displaystyle \frac{{\beta }_{{l}^{\prime }}}{2l\prime +1}{c}_{l\prime m\prime }{Y}_{l\prime m\prime }(\hat{{\bf{s}}}).\end{eqnarray} \tag{ 23 }$$

Thus for $| m| \leqslant L$ we have

$$\begin{eqnarray*}&&{c}_{{lm}}-\displaystyle \frac{1}{\nu }\displaystyle \sum _{l\prime m\prime }{c}_{l\prime m\prime }{\displaystyle \int }_{{{\mathbb{S}}}^{2}}\mu {Y}_{l\prime m\prime }(\hat{{\bf{s}}}){Y}_{{lm}}^{*}(\hat{{\bf{s}}})\;{\rm{d}}\hat{{\bf{s}}}=\displaystyle \frac{\varpi {\beta }_{l}}{2l+1}{\chi }_{[0,L]}(l){c}_{{lm}}.\end{eqnarray*}$$

Using the orthogonality relation for associated Legendre polynomials: ${\displaystyle \int }_{-1}^{1}{P}_{l}^{m}(\mu ){P}_{{l}^{\prime }}^{m}\;{\rm{d}}\mu ={\delta }_{{{ll}}^{\prime }}2(l+m)!/[(2l+1)(l-m)!]$, we obtain

$$\begin{eqnarray*}&&\begin{array}{l}\sqrt{\displaystyle \frac{(2l+1)(2l\prime +1)}{{h}_{l}{h}_{{l}^{\prime }}}}\left(\sqrt{\displaystyle \frac{{l}^{2}-{m}^{2}}{4{l}^{2}-1}}{\delta }_{l\prime ,l-1}+\sqrt{\displaystyle \frac{{(l+1)}^{2}-{m}^{2}}{4{(l+1)}^{2}-1}}{\delta }_{l\prime ,l+1}\right){c}_{l\prime m}\sqrt{\displaystyle \frac{{h}_{{l}^{\prime }}}{2l\prime +1}}\\ \qquad =\left({b}_{l}(m){\delta }_{l\prime ,l-1}+{b}_{{l}^{\prime }}(m){\delta }_{l\prime ,l+1}\right){c}_{l\prime m}\sqrt{\displaystyle \frac{{h}_{{l}^{\prime }}}{2l\prime +1}}=\nu {c}_{{lm}}\sqrt{\displaystyle \frac{{h}_{l}}{2l+1}}.\end{array}\end{eqnarray*}$$

The above equation forms an eigenvalue problem for $B(m)$, and c_lm are given in terms of eigenvectors of $B(m)$.

2.3. Singular eigenfunctions for three dimensions

(Rotated reference frames).

Definition 2.5 Let $\hat{{\bf{k}}}\in {{\mathbb{C}}}^{3}$ be a unit vector such that $\hat{{\bf{k}}}\cdot \hat{{\bf{k}}}=1$. We define an invertible linear operator ${{\mathcal{R}}}_{\hat{{\bf{k}}}}\;:{\mathbb{C}}\mapsto {\mathbb{C}}$. For a function ${f}_{1}(\hat{{\bf{s}}})\in {\mathbb{C}}$ ($\hat{{\bf{s}}}\in {{\mathbb{S}}}^{2}$), ${{\mathcal{R}}}_{\hat{{\bf{k}}}}{f}_{1}(\hat{{\bf{s}}})$ is the value of ${f}_{1}(\hat{{\bf{s}}})$ where $\hat{{\bf{s}}}$ is measured in the rotated reference frame whose z-axis lies in the direction of $\hat{{\bf{k}}}$.

Suppose that ${f}_{1}(\hat{{\bf{s}}})\in {\mathbb{C}}$ is given by spherical harmonics:

$$\begin{eqnarray*}&&{f}_{1}(\hat{{\bf{s}}})=\displaystyle \sum _{l=0}^{\infty }\displaystyle \sum _{m=-l}^{l}{f}_{{lm}}{Y}_{{lm}}(\hat{{\bf{s}}}),\end{eqnarray*}$$

where ${f}_{{lm}}\in {\mathbb{C}}$. Then we have [8, 24, 43]

$$\begin{eqnarray*}\begin{array}{rcl}{{\mathcal{R}}}_{\hat{{\bf{k}}}}{f}_{1}(\hat{{\bf{s}}}) & = & \displaystyle \sum _{l=0}^{\infty }\displaystyle \sum _{m=-l}^{l}{f}_{{lm}}\displaystyle \sum _{m\prime =-l}^{l}{D}_{m\prime m}^{l}({\varphi }_{\hat{{\bf{k}}}},{\theta }_{\hat{{\bf{k}}}},0){Y}_{{{lm}}^{\prime }}(\hat{{\bf{s}}})\\ & = & \displaystyle \sum _{l=0}^{\infty }\displaystyle \sum _{m=-l}^{l}{f}_{{lm}}\displaystyle \sum _{m\prime =-l}^{l}{{\rm{e}}}^{-{\rm{i}}m\prime {\varphi }_{\hat{{\bf{k}}}}}{d}_{m\prime m}^{l}({\theta }_{\hat{{\bf{k}}}}){Y}_{{{lm}}^{\prime }}(\hat{{\bf{s}}}),\end{array}\end{eqnarray*}$$

where ${\theta }_{\hat{{\bf{k}}}}$ and ${\varphi }_{\hat{{\bf{k}}}}$ are the polar and azimuthal angles of $\hat{{\bf{k}}}$ in the laboratory frame. Here ${D}_{m\prime m}^{l}$ and ${d}_{m\prime m}^{l}$ are Wignerʼs D-matrices and d-matrices. Moreover we obtain

$$\begin{eqnarray*}\begin{array}{rcl}{{\mathcal{R}}}_{\hat{{\bf{k}}}}^{-1}{f}_{1}(\hat{{\bf{s}}}) & = & \displaystyle \sum _{l=0}^{\infty }\displaystyle \sum _{m=-l}^{l}{f}_{{lm}}\displaystyle \sum _{m\prime =-l}^{l}{D}_{m\prime m}^{l}(0,-{\theta }_{\hat{{\bf{k}}}},-{\varphi }_{\hat{{\bf{k}}}}){Y}_{{{lm}}^{\prime }}(\hat{{\bf{s}}})\\ & = & \displaystyle \sum _{l=0}^{\infty }\displaystyle \sum _{m=-l}^{l}{f}_{{lm}}\displaystyle \sum _{m\prime =-l}^{l}{{\rm{e}}}^{{\rm{i}}m{\varphi }_{\hat{{\bf{k}}}}}{d}_{{{mm}}^{\prime }}^{l}({\theta }_{\hat{{\bf{k}}}}){Y}_{{{lm}}^{\prime }}(\hat{{\bf{s}}}).\end{array}\end{eqnarray*}$$

We can directly show ${{\mathcal{R}}}_{\hat{{\bf{k}}}}^{-1}{{\mathcal{R}}}_{\hat{{\bf{k}}}}{f}_{1}(\hat{{\bf{s}}})={f}_{1}(\hat{{\bf{s}}})$ by using $\displaystyle \sum _{m\prime =-l}^{l}{d}_{m\prime m}^{l}({\theta }_{\hat{{\bf{k}}}}){d}_{m\prime m\prime \prime }^{l}({\theta }_{\hat{{\bf{k}}}})={\delta }_{{{mm}}^{\prime\prime }}$. We have for ${f}_{1}(\hat{{\bf{s}}}),{f}_{2}(\hat{{\bf{s}}})\in {\mathbb{C}}$,

$$\begin{eqnarray*}&&{{\mathcal{R}}}_{\hat{{\bf{k}}}}{f}_{1}(\hat{{\bf{s}}}){f}_{2}(\hat{{\bf{s}}})=\left({{\mathcal{R}}}_{\hat{{\bf{k}}}}{f}_{1}(\hat{{\bf{s}}})\right)\left({{\mathcal{R}}}_{\hat{{\bf{k}}}}{f}_{1}(\hat{{\bf{s}}})\right),\qquad {\displaystyle \int }_{{{\mathbb{S}}}^{2}}{{\mathcal{R}}}_{\hat{{\bf{k}}}}{f}_{1}(\hat{{\bf{s}}})\;{\rm{d}}\hat{{\bf{s}}}={\displaystyle \int }_{{{\mathbb{S}}}^{2}}{f}_{1}(\hat{{\bf{s}}})\;{\rm{d}}\hat{{\bf{s}}}.\end{eqnarray*}$$

Example 2.6. For any function ${f}_{1}(\hat{{\bf{s}}})$ and the unit vector $\hat{{\bf{z}}}$ in the positive direction on the z-axis, we have ${{\mathcal{R}}}_{\hat{{\bf{z}}}}{f}_{1}(\hat{{\bf{s}}})={f}_{1}(\hat{{\bf{s}}})$.

Example 2.7.  ${{\mathcal{R}}}_{\hat{{\bf{k}}}}\hat{{\bf{s}}}\cdot \hat{{\bf{s}}}\prime =\hat{{\bf{s}}}\cdot \hat{{\bf{s}}}\prime$ for $\hat{{\bf{s}}},\hat{{\bf{s}}}\prime \in {{\mathbb{S}}}^{2}$.

Example 2.8.  ${{\mathcal{R}}}_{\hat{{\bf{k}}}}\mu =\sqrt{\frac{4\pi }{3}}{{\mathcal{R}}}_{\hat{{\bf{k}}}}{Y}_{10}(\hat{{\bf{s}}})=\displaystyle \sum _{m\prime =-1}^{1}{{\rm{e}}}^{-{\rm{i}}m\prime {\varphi }_{\hat{{\bf{k}}}}}{d}_{m\prime 0}^{1}({\theta }_{\hat{{\bf{k}}}}){Y}_{1{m}^{\prime }}(\hat{{\bf{s}}})=\hat{{\bf{s}}}\cdot \hat{{\bf{k}}}$.

(Plane wave decomposition).

Definition 2.9 Complex unit vectors $\hat{{\bf{k}}}(\nu ,{\bf{q}})\in {{\mathbb{C}}}^{3}$ ($\nu \in {\mathbb{R}}$, ${\bf{q}}\in {{\mathbb{R}}}^{2}$) are given by

$$\begin{eqnarray*}&&\hat{{\bf{k}}}(\nu ,{\bf{q}})=\left(\begin{array}{c}{\rm{i}}\nu {\bf{q}}\\ {\hat{k}}_{z}(\nu q)\end{array}\right),\end{eqnarray*}$$

where $q=| {\bf{q}}|$ and

$$\begin{eqnarray*}&&{\hat{k}}_{z}(\nu q)=\sqrt{1+{(\nu q)}^{2}}.\end{eqnarray*}$$

Example 2.10. For $\nu \in {\mathbb{R}}$, ${\bf{q}}\in {{\mathbb{R}}}^{2}$, we obtain

$$\begin{eqnarray}&&{{\mathcal{R}}}_{\hat{{\bf{k}}}(-\nu ,{\bf{q}})}\mu ={\hat{k}}_{z}(\nu q)\mu -{\rm{i}}\nu q\sqrt{1-{\mu }^{2}}\mathrm{cos}(\varphi -{\varphi }_{{\bf{q}}}),\end{eqnarray} \tag{ 24 }$$

$$\begin{eqnarray}&&{{\mathcal{R}}}_{\hat{{\bf{k}}}(-\nu ,{\bf{q}})}^{-1}\mu =\sqrt{\displaystyle \frac{4\pi }{3}}{{\mathcal{R}}}_{\hat{{\bf{k}}}(-\nu ,{\bf{q}})}^{-1}{Y}_{10}(\hat{{\bf{s}}})={\hat{k}}_{z}(\nu q)\mu -{\rm{i}}| \nu q| \sqrt{1-{\mu }^{2}}\mathrm{cos}\varphi .\end{eqnarray} \tag{ 25 }$$

([41, 48]).

Definition 2.11 We define

$$\begin{eqnarray*}&&\mathrm{cos}[i\tau (\nu q)]=\mathrm{cos}{\theta }_{\hat{{\bf{k}}}(\nu ,{\bf{q}})},\qquad \mathrm{sin}[i\tau (| \nu q| )]=\mathrm{sin}{\theta }_{\hat{{\bf{k}}}(\nu ,{\bf{q}})}.\end{eqnarray*}$$

Since Wignerʼs d-matrices ${d}_{{{mm}}^{\prime }}^{l}(\theta )$ are given in terms of $\mathrm{cos}\theta$, we also write

$$\begin{eqnarray*}&&{d}_{{{mm}}^{\prime }}^{l}[i\tau (\nu q)]={d}_{{{mm}}^{\prime }}^{l}({\theta }_{\hat{{\bf{k}}}(\nu ,{\bf{q}})}).\end{eqnarray*}$$

To compute Wignerʼs d-matrices, we take square roots such that $0\leqslant \mathrm{arg}(\sqrt{z})\lt \pi$ for all $z\in {\mathbb{C}}$ [41, 48]. We have

$$\begin{eqnarray*}&&\mathrm{cos}{\varphi }_{\hat{{\bf{k}}}(\nu ,{\bf{q}})}=\displaystyle \frac{\hat{{\bf{x}}}\cdot ({\rm{i}}\nu {\bf{q}})}{\sqrt{({\rm{i}}\nu {\bf{q}})\cdot ({\rm{i}}\nu {\bf{q}})}}=\displaystyle \frac{\nu }{| \nu | }\mathrm{cos}{\varphi }_{{\bf{q}}},\qquad \mathrm{sin}{\varphi }_{\hat{{\bf{k}}}(\nu ,{\bf{q}})}=\displaystyle \frac{\hat{{\bf{y}}}\cdot ({\rm{i}}\nu {\bf{q}})}{\sqrt{({\rm{i}}\nu {\bf{q}})\cdot ({\rm{i}}\nu {\bf{q}})}}=\displaystyle \frac{\nu }{| \nu | }\mathrm{sin}{\varphi }_{{\bf{q}}},\end{eqnarray*}$$

and

$$\begin{eqnarray*}&&\mathrm{cos}{\theta }_{\hat{{\bf{k}}}(\nu ,{\bf{q}})}=\hat{{\bf{z}}}\cdot \hat{{\bf{k}}}={\hat{k}}_{z}(\nu q),\qquad \mathrm{sin}{\theta }_{\hat{{\bf{k}}}(\nu ,{\bf{q}})}=\sqrt{1-{\mathrm{cos}}^{2}{\theta }_{\hat{{\bf{k}}}(\nu ,{\bf{q}})}}=\sqrt{({\rm{i}}\nu {\bf{q}})\cdot ({\rm{i}}\nu {\bf{q}})}={\rm{i}}| \nu q| .\end{eqnarray*}$$

In particular we obtain

$$\begin{eqnarray*}{\varphi }_{\hat{{\bf{k}}}(\nu ,{\bf{q}})}=\left\{\begin{array}{ll}{\varphi }_{{\bf{q}}}, & \mathrm{for}\ \nu \gt 0,\\ {\varphi }_{{\bf{q}}}+\pi , & \mathrm{for}\ \nu \lt 0,\end{array}\right.\end{eqnarray*}$$

where ${\varphi }_{{\bf{q}}}$ is the polar angle of ${\bf{q}}$.

Let us consider

$$\begin{eqnarray}&&\hat{{\bf{s}}}\cdot \nabla I({\bf{r}},\hat{{\bf{s}}})+I({\bf{r}},\hat{{\bf{s}}})=\varpi {\displaystyle \int }_{{{\mathbb{S}}}^{2}}p(\hat{{\bf{s}}},\hat{{\bf{s}}}\prime )I({\bf{r}},\hat{{\bf{s}}}\prime )\;{\rm{d}}\hat{{\bf{s}}}\prime ,\end{eqnarray} \tag{ 26 }$$

where ${\bf{r}}\in {{\mathbb{R}}}^{3}$, $\hat{{\bf{s}}}\in {{\mathbb{S}}}^{2}$. Solutions to the above equation are given by a superposition of eigenmodes

$$\begin{eqnarray}&&I({\bf{r}},\hat{{\bf{s}}})={{\mathcal{R}}}_{\hat{{\bf{k}}}}{\Phi }_{\nu }^{m}(\hat{{\bf{s}}}){{\rm{e}}}^{-\hat{{\bf{k}}}\cdot {\bf{r}}/\nu },\end{eqnarray} \tag{ 27 }$$

where $\hat{{\bf{k}}}=\hat{{\bf{k}}}(\nu ,{\bf{q}})$. To see this we substitute the separated solution (27) in the above homogeneous three-dimensional radiative transport equation (26) and obtain

$$\begin{eqnarray}&&\left(1-\displaystyle \frac{{{\mathcal{R}}}_{\hat{{\bf{k}}}}\mu }{\nu }\right){{\mathcal{R}}}_{\hat{{\bf{k}}}}{\Phi }_{\nu }^{m}(\hat{{\bf{s}}})=\varpi {\displaystyle \int }_{{{\mathbb{S}}}^{2}}p(\hat{{\bf{s}}},\hat{{\bf{s}}}\prime ){{\mathcal{R}}}_{\hat{{\bf{k}}}}{\Phi }_{\nu }^{m}(\hat{{\bf{s}}}\prime )\;{\rm{d}}\hat{{\bf{s}}}\prime .\end{eqnarray} \tag{ 28 }$$

The right-hand side can be written as

$$\begin{eqnarray}&&\varpi {\displaystyle \int }_{{{\mathbb{S}}}^{2}}p(\hat{{\bf{s}}},\hat{{\bf{s}}}\prime ){{\mathcal{R}}}_{\hat{{\bf{k}}}}{\Phi }_{\nu }^{m}(\hat{{\bf{s}}}\prime )\;{\rm{d}}\hat{{\bf{s}}}\prime =\varpi {\displaystyle \int }_{{{\mathbb{S}}}^{2}}p({{\mathcal{R}}}_{\hat{{\bf{k}}}}\hat{{\bf{s}}},{{\mathcal{R}}}_{\hat{{\bf{k}}}}\hat{{\bf{s}}}\prime ){{\mathcal{R}}}_{\hat{{\bf{k}}}}{\Phi }_{\nu }^{m}(\hat{{\bf{s}}}\prime )\;{\rm{d}}\hat{{\bf{s}}}\prime .\end{eqnarray} \tag{ 29 }$$

That is,

$$\begin{eqnarray*}&&{{\mathcal{R}}}_{\hat{{\bf{k}}}}\left(1-\displaystyle \frac{\mu }{\nu }\right){\Phi }_{\nu }^{m}(\hat{{\bf{s}}})={{\mathcal{R}}}_{\hat{{\bf{k}}}}\varpi {\displaystyle \int }_{{{\mathbb{S}}}^{2}}p(\hat{{\bf{s}}},\hat{{\bf{s}}}\prime ){\Phi }_{\nu }^{m}(\hat{{\bf{s}}}\prime )\;{\rm{d}}\hat{{\bf{s}}}\prime .\end{eqnarray*}$$

Thus the three-dimensional equation (28) reduces to the one-dimensional equation (13). Recall that ${\Phi }_{\nu }^{m}(\hat{{\bf{s}}})$ given in (12) is constructed so that (13) obtained from (10) and (11) is satisfied. We have

$$\begin{eqnarray}&&{{\mathcal{R}}}_{\hat{{\bf{k}}}}{\phi }^{m}(\nu ,\mu )=\displaystyle \frac{\varpi \nu }{2}{\mathcal{P}}\displaystyle \frac{{g}^{m}(\nu ,\hat{{\bf{s}}}\cdot \hat{{\bf{k}}})}{\nu -\hat{{\bf{s}}}\cdot \hat{{\bf{k}}}}+{\lambda }^{m}(\nu ){\left(1-{\nu }^{2}\right)}^{-| m| }\delta (\nu -\hat{{\bf{s}}}\cdot \hat{{\bf{k}}}).\end{eqnarray} \tag{ 30 }$$

Proposition 2.12. The following orthogonality relation holds.

$$\begin{eqnarray*}&&{\displaystyle \int }_{{{\mathbb{S}}}^{2}}\mu \left({{\mathcal{R}}}_{\hat{{\bf{k}}}(\nu ,{\bf{q}})}{\Phi }_{\nu }^{m}(\hat{{\bf{s}}})\right)\left({{\mathcal{R}}}_{\hat{{\bf{k}}}({\nu }^{\prime },{\bf{q}})}{\Phi }_{{\nu }^{\prime }}^{{m\prime }^{*}}(\hat{{\bf{s}}})\right)\;{\rm{d}}\hat{{\bf{s}}}=2\pi {\hat{k}}_{z}(\nu q){\mathcal{N}}(\nu ){\delta }_{\nu \nu \prime }{\delta }_{{{mm}}^{\prime }}.\end{eqnarray*}$$

Proof. The full-range orthogonality is obtained in [42] through the Greenʼs function. Here we give a direct proof.

We perform separation of variables to the homogeneous equation by assuming the form (27). By substituting the separated solution into the radiative transport equation (26), we obtain

$$\begin{eqnarray*}&&\left(1-\displaystyle \frac{{{\mathcal{R}}}_{\hat{{\bf{k}}}}\mu }{\nu }\right){{\mathcal{R}}}_{\hat{{\bf{k}}}}{\Phi }_{\nu }^{m}(\hat{{\bf{s}}})=\varpi \displaystyle \sum _{l=0}^{L}\displaystyle \sum _{m=-l}^{l}\displaystyle \frac{{\beta }_{l}}{2l+1}{Y}_{{lm}}(\hat{{\bf{s}}}){\displaystyle \int }_{{{\mathbb{S}}}^{2}}{Y}_{{lm}}^{*}(\hat{{\bf{s}}}\prime ){{\mathcal{R}}}_{\hat{{\bf{k}}}}{\Phi }_{\nu }^{m}(\hat{{\bf{s}}}\prime )\;{\rm{d}}\hat{{\bf{s}}}\prime ,\end{eqnarray*}$$

For fixed ${\bf{q}}$, we consider $({m}_{1},{\nu }_{1}^{{m}_{1}})$ and $({m}_{2},{\nu }_{2}^{{m}_{2}})$. Let us write ${\hat{{\bf{k}}}}_{1}=\hat{{\bf{k}}}({\nu }_{1}^{{m}_{1}},{\bf{q}})$, ${\hat{{\bf{k}}}}_{2}=\hat{{\bf{k}}}({\nu }_{2}^{{m}_{2}},{\bf{q}})$. We write the following two equations.

$$\begin{eqnarray*}&&\begin{array}{l}\left({{\mathcal{R}}}_{{\hat{{\bf{k}}}}_{2}}{\Phi }_{{\nu }_{2}}^{{m}_{2}}(\hat{{\bf{s}}})\right){{\mathcal{R}}}_{{\hat{{\bf{k}}}}_{1}}\left(1-\displaystyle \frac{\mu }{{\nu }_{1}}\right){\Phi }_{{\nu }_{1}}^{{m}_{1}}(\hat{{\bf{s}}})\\ \qquad =\varpi \displaystyle \sum _{l=0}^{L}\displaystyle \sum _{m=-l}^{l}\displaystyle \frac{{\beta }_{l}}{2l+1}{Y}_{{lm}}(\hat{{\bf{s}}})\left({{\mathcal{R}}}_{{\hat{{\bf{k}}}}_{2}}{\Phi }_{{\nu }_{2}}^{{m}_{2}}(\hat{{\bf{s}}})\right){\displaystyle \int }_{{{\mathbb{S}}}^{2}}{Y}_{{lm}}^{*}(\hat{{\bf{s}}}\prime )\left({{\mathcal{R}}}_{{\hat{{\bf{k}}}}_{1}}{\Phi }_{{\nu }_{1}}^{{m}_{1}}(\hat{{\bf{s}}}\prime )\right)\;{\rm{d}}\hat{{\bf{s}}}\prime ,\\ \left({{\mathcal{R}}}_{{\hat{{\bf{k}}}}_{1}}{\Phi }_{{\nu }_{1}}^{{m}_{1}}(\hat{{\bf{s}}})\right){{\mathcal{R}}}_{{\hat{{\bf{k}}}}_{2}}\left(1-\displaystyle \frac{\mu }{{\nu }_{2}}\right){\Phi }_{{\nu }_{2}}^{{m}_{2}}(\hat{{\bf{s}}})\\ \qquad =\varpi \displaystyle \sum _{l=0}^{L}\displaystyle \sum _{m=-l}^{l}\displaystyle \frac{{\beta }_{l}}{2l+1}{Y}_{{lm}}^{*}(\hat{{\bf{s}}})\left({{\mathcal{R}}}_{{\hat{{\bf{k}}}}_{1}}{\Phi }_{{\nu }_{1}}^{{m}_{1}}(\hat{{\bf{s}}})\right){\displaystyle \int }_{{{\mathbb{S}}}^{2}}{Y}_{{lm}}(\hat{{\bf{s}}}\prime )\left({{\mathcal{R}}}_{{\hat{{\bf{k}}}}_{2}}{\Phi }_{{\nu }_{2}}^{{m}_{2}}(\hat{{\bf{s}}}\prime )\right)\;{\rm{d}}\hat{{\bf{s}}}\prime .\end{array}\end{eqnarray*}$$

We note (24). By subtraction and integration over $\hat{{\boldsymbol{}}s}$ we have

$$\begin{eqnarray*}&&{\displaystyle \int }_{{{\mathbb{S}}}^{2}}\left(\displaystyle \frac{{{\mathcal{R}}}_{{\hat{{\bf{k}}}}_{2}}\mu }{{\nu }_{2}}-\displaystyle \frac{{{\mathcal{R}}}_{{\hat{{\bf{k}}}}_{1}}\mu }{{\nu }_{1}}\right)\left({{\mathcal{R}}}_{{\hat{{\bf{k}}}}_{1}}{\Phi }_{{\nu }_{1}}^{{m}_{1}}(\hat{{\boldsymbol{}}{\bf{s}}})\right)\left({{\mathcal{R}}}_{{\hat{{\bf{k}}}}_{2}}{\Phi }_{{\nu }_{2}}^{{m}_{2}}(\hat{{\boldsymbol{}}{\bf{s}}})\right)\;{\rm{d}}\hat{{\bf{s}}}\end{eqnarray*}$$

$$\begin{eqnarray*}&&\qquad =\;\left(\displaystyle \frac{{\hat{k}}_{z}({\nu }_{2}q)}{{\nu }_{2}}-\displaystyle \frac{{\hat{k}}_{z}({\nu }_{1}q)}{{\nu }_{1}}\right){\displaystyle \int }_{{{\mathbb{S}}}^{2}}\mu \left({{\mathcal{R}}}_{{\hat{{\bf{k}}}}_{1}}{\Phi }_{{\nu }_{1}}^{{m}_{1}}(\hat{{\bf{s}}})\right)\left({{\mathcal{R}}}_{{\hat{{\bf{k}}}}_{2}}{\Phi }_{{\nu }_{2}}^{{m}_{2}}(\hat{{\bf{s}}})\right)\;{\rm{d}}\hat{{\bf{s}}}\end{eqnarray*}$$

$$\begin{eqnarray*}&&\qquad =\;0.\end{eqnarray*}$$

Therefore,

$$\begin{eqnarray}&&{\displaystyle \int }_{{{\mathbb{S}}}^{2}}\mu \left({{\mathcal{R}}}_{{\hat{{\bf{k}}}}_{1}}{\Phi }_{{\nu }_{1}}^{{m}_{1}}(\hat{{\bf{s}}})\right)\left({{\mathcal{R}}}_{{\hat{{\bf{k}}}}_{2}}{\Phi }_{{\nu }_{2}}^{{m}_{2}}(\hat{{\bf{s}}})\right)\;{\rm{d}}\hat{{\bf{s}}}=0\quad {\rm{for}}\ {\nu }_{1}\ne {\nu }_{2}.\end{eqnarray} \tag{ 31 }$$

Suppose $\nu ={\nu }_{1}={\nu }_{2}$, $\hat{{\bf{k}}}={\hat{{\bf{k}}}}_{1}={\hat{{\bf{k}}}}_{2}$, ${m}_{1}\ne {m}_{2}$. In this case we have

$$\begin{eqnarray}&&\begin{array}{l}{\displaystyle \int }_{{{\mathbb{S}}}^{2}}\mu \left({{\mathcal{R}}}_{\hat{{\bf{k}}}}{\Phi }_{\nu }^{{m}_{1}}(\hat{{\bf{s}}})\right)\left({{\mathcal{R}}}_{\hat{{\bf{k}}}}{\Phi }_{\nu }^{{m}_{2}}(\hat{{\bf{s}}})\right)\;{\rm{d}}\hat{{\boldsymbol{}}{\bf{s}}}={\displaystyle \int }_{{{\mathbb{S}}}^{2}}\mu {\Phi }_{\nu }^{{m}_{1}}(\hat{{\boldsymbol{}}{\bf{s}}}){\Phi }_{\nu }^{{m}_{2}}(\hat{{\boldsymbol{}}{\bf{s}}})\;{\rm{d}}\hat{{\boldsymbol{}}{\bf{s}}}\\ \qquad ={\displaystyle \int }_{0}^{2\pi }{{\rm{e}}}^{{\rm{i}}({m}_{1}+{m}_{2})\varphi }\;{\rm{d}}\varphi \ {\displaystyle \int }_{-1}^{1}\mu {\phi }^{{m}_{1}}(\nu ,\mu ){\phi }^{{m}_{2}}(\nu ,\mu ){\left(1-{\mu }^{2}\right)}^{(| {m}_{1}| +| {m}_{2}| )/2}\;{\rm{d}}\mu \\ \qquad \propto {\delta }_{{m}_{1},-{m}_{2}}.\end{array}\end{eqnarray} \tag{ 32 }$$

We note that ${\bf{q}}$

$$\begin{eqnarray*}&&{\Phi }_{\nu }^{-m}(\hat{{\bf{s}}})={\Phi }_{\nu }^{m}{(\hat{{\bf{s}}})}^{*}.\end{eqnarray*}$$

Using (31) and (32), for arbitrary $\nu ,{\nu }^{\prime },m,m\prime$, we have

$$\begin{eqnarray*}&&{\displaystyle \int }_{{{\mathbb{S}}}^{2}}\mu \left({{\mathcal{R}}}_{\hat{{\bf{k}}}}{\Phi }_{\nu }^{m}(\hat{{\bf{s}}})\right)\left({{\mathcal{R}}}_{\hat{{\bf{k}}}\prime }{\Phi }_{{\nu }^{\prime }}^{{m}^{\prime }}{(\hat{{\bf{s}}})}^{*}\right)\;{\rm{d}}\hat{{\bf{s}}}\propto {\delta }_{\nu \nu \prime }{\delta }_{{{mm}}^{\prime }}.\end{eqnarray*}$$

If $\nu =\nu \prime$, $m=m\prime$, we have

$$\begin{eqnarray*}\begin{array}{rcl}{\displaystyle \int }_{{{\mathbb{S}}}^{2}}\mu \left({{\mathcal{R}}}_{\hat{{\bf{k}}}}{\Phi }_{\nu }^{m}(\hat{{\boldsymbol{}}{\bf{s}}})\right)\left({{\mathcal{R}}}_{\hat{{\bf{k}}}}{\Phi }_{\nu }^{m}{(\hat{{\boldsymbol{}}{\bf{s}}})}^{*}\right)\;{\rm{d}}\hat{{\boldsymbol{}}{\bf{s}}} & = & {\displaystyle \int }_{{{\mathbb{S}}}^{2}}\left({{\mathcal{R}}}_{\hat{{\bf{k}}}}^{-1}\mu \right){\Phi }_{\nu }^{m}(\hat{{\boldsymbol{}}{\bf{s}}}){\Phi }_{\nu }^{m}{(\hat{{\bf{s}}})}^{*}\;{\rm{d}}\hat{{\boldsymbol{}}s}\\ & = & {\displaystyle \int }_{{{\mathbb{S}}}^{2}}\left({{\mathcal{R}}}_{\hat{{\bf{k}}}}^{-1}\mu \right){\left[{\phi }^{m}(\nu ,\mu )\right]}^{2}{\left(1-{\mu }^{2}\right)}^{| m| }\;{\rm{d}}\hat{{\boldsymbol{}}{\bf{s}}}.\end{array}\end{eqnarray*}$$

Hence,

$$\begin{eqnarray*}\begin{array}{rcl}{\displaystyle \int }_{{{\mathbb{S}}}^{2}}\mu \left({{\mathcal{R}}}_{\hat{{\bf{k}}}}{\Phi }_{\nu }^{m}(\hat{{\bf{s}}})\right)\left({{\mathcal{R}}}_{\hat{{\bf{k}}}}{\Phi }_{\nu }^{m}{(\hat{{\bf{s}}})}^{*}\right)\;{\rm{d}}\hat{{\bf{s}}} & = & 2\pi {\hat{k}}_{z}(\nu q){\displaystyle \int }_{-1}^{1}\mu {\left[{\phi }^{m}(\nu ,\mu )\right]}^{2}{\left(1-{\mu }^{2}\right)}^{| m| }\;{\rm{d}}\mu \\ & = & 2\pi {\hat{k}}_{z}(\nu q){{\mathcal{N}}}^{m}(\nu ),\end{array}\end{eqnarray*}$$

where the normalization factor ${{\mathcal{N}}}^{m}(\nu )$ is given in (17). Thus we obtain the full-range orthogonality relation. □

2.4. Method of rotated reference frames

The method of rotated reference frames does not rely on singular eigenfunctions ${\Phi }_{\nu }^{{m}^{\prime }}(\hat{{\bf{s}}})$ and uses the expansion (5), in which ${c}_{{lm}}^{{m}^{\prime }}(\nu )$ are unknown coefficients that can be fully numerically computed as eigenvectors of $B(m\prime )$. The method is summarized in appendix A. We describe below how the matrix $B(m\prime )$ appears in this method.

We plug (5) into (28):

$$\begin{eqnarray*}&&\left(1-\displaystyle \frac{{{\mathcal{R}}}_{\hat{{\bf{k}}}}\mu }{\nu }\right)\displaystyle \sum _{{lm}}{c}_{{lm}}^{{m}^{\prime }}{{\mathcal{R}}}_{\hat{{\bf{k}}}}{Y}_{{lm}}(\hat{{\bf{s}}})=\varpi {\displaystyle \int }_{{{\mathbb{S}}}^{2}}p({{\mathcal{R}}}_{\hat{{\bf{k}}}}\hat{{\bf{s}}},{{\mathcal{R}}}_{\hat{{\bf{k}}}}\hat{{\bf{s}}}\prime )\displaystyle \sum _{{lm}}{c}_{{lm}}^{{m}^{\prime }}{{\mathcal{R}}}_{\hat{{\bf{k}}}}{Y}_{{lm}}(\hat{{\bf{s}}}\prime )\;{\rm{d}}\hat{{\bf{s}}}\prime .\end{eqnarray*}$$

By operating ${{\mathcal{R}}}_{\hat{{\bf{k}}}}^{-1}$, the above equation reduces to (23), from which the matrix $B(m\prime )$ is derived.

To show how the F_N method can be extended to three dimensions, we will consider the half-space geometry in which a homogeneous random medium with optical parameter ϖ exists only in the lower half $z\lt 0$. By the Placzek lemma [5] we can consider the following radiative transport equation in ${{\mathbb{R}}}^{3}$ instead of (1).

$$\begin{eqnarray*}\left\{\begin{array}{ll}\hat{{\bf{s}}}\cdot \nabla \psi ({\bf{r}},\hat{{\bf{s}}})+\psi ({\bf{r}},\hat{{\bf{s}}})=\varpi {\displaystyle \int }_{{{\mathbb{S}}}^{2}}p(\hat{{\bf{s}}},\hat{{\bf{s}}}\prime )\psi ({\bf{r}},\hat{{\bf{s}}}\prime )\;{\rm{d}}\hat{{\bf{s}}}\prime & \\ \qquad +{\chi }_{(0,\infty )}(z)S({\bf{r}},\hat{{\bf{s}}})+\mu I({\bf{r}},\hat{{\bf{s}}})\delta (z), & z\in (-\infty ,\infty ),\\ \psi ({\bf{r}},\hat{{\bf{s}}})\to 0, & | z| \to \infty ,\end{array}\right.\end{eqnarray*}$$

where ${\chi }_{(0,\infty )}(z)=1$ for $z\gt 0$ and $=\;0$ otherwise. We have the jump condition

$$\begin{eqnarray*}&&\psi ({\boldsymbol{\rho }},{0}^{+},\hat{{\bf{s}}})-\psi ({\boldsymbol{\rho }},{0}^{-},\hat{{\bf{s}}})=I({\boldsymbol{\rho }},0,\hat{{\bf{s}}}).\end{eqnarray*}$$

Since $I({\boldsymbol{\rho }},0,\hat{{\bf{s}}})$ is given by only eigenmodes with positive eigenvalues and $\psi ({\boldsymbol{\rho }},{0}^{-},\hat{{\bf{s}}})$ is given by only eigenmodes with negative eigenvalues, we see that $\psi ({\boldsymbol{\rho }},{0}^{-},\hat{{\bf{s}}})=0$. Therefore we obtain the relation

$$\begin{eqnarray*}\psi ({\bf{r}},\hat{{\bf{s}}})=\left\{\begin{array}{ll}I({\bf{r}},\hat{{\bf{s}}}), & z\gt 0,\\ 0, & z\lt 0.\end{array}\right.\end{eqnarray*}$$

Let us introduce the Greenʼs function $G({\bf{r}},\hat{{\bf{s}}};{{\bf{r}}}_{0},{\hat{{\bf{s}}}}_{0})$ for the infinite medium as

$$\begin{eqnarray*}\left\{\begin{array}{ll}\hat{{\bf{s}}}\cdot \nabla G({\bf{r}},\hat{{\bf{s}}})+G({\bf{r}},\hat{{\bf{s}}})=\varpi {\displaystyle \int }_{{{\mathbb{S}}}^{2}}p(\hat{{\bf{s}}},\hat{{\bf{s}}}\prime )G({\bf{r}},\hat{{\bf{s}}}\prime )\;{\rm{d}}\hat{{\bf{s}}}\prime & \\ \qquad +\delta ({\bf{r}}-{{\bf{r}}}_{0})\delta (\hat{{\bf{s}}}-{\hat{{\bf{s}}}}_{0}), & z\in (-\infty ,\infty ),\\ G({\bf{r}},\hat{{\bf{s}}})\to 0, & | z| \to \infty .\end{array}\right.\end{eqnarray*}$$

Thus we obtain

$$\begin{eqnarray*}\begin{array}{rcl}\psi ({\bf{r}},\hat{{\bf{s}}}) & = & {\displaystyle \int }_{{{\mathbb{S}}}^{2}}{\displaystyle \int }_{{{\mathbb{R}}}^{2}}G({\bf{r}},\hat{{\bf{s}}};{{\boldsymbol{\rho }}}^{\prime },0,\hat{{\bf{s}}}\prime )\mu \prime I({{\boldsymbol{\rho }}}^{\prime },0,\hat{{\bf{s}}}\prime )\;d{\boldsymbol{\rho }}\prime {\rm{d}}\hat{{\bf{s}}}\prime \\ & & +{\displaystyle \int }_{{{\mathbb{S}}}^{2}}{\displaystyle \int }_{0}^{\infty }{\displaystyle \int }_{{{\mathbb{R}}}^{2}}G({\bf{r}},\hat{{\bf{s}}};{{\boldsymbol{\rho }}}^{\prime },z\prime ,\hat{{\bf{s}}}\prime )S({{\boldsymbol{\rho }}}^{\prime },z\prime ,\hat{{\bf{s}}}\prime )\;d{\boldsymbol{\rho }}\prime {\rm{d}}z\prime {\rm{d}}\hat{{\bf{s}}}\prime ,\quad {\bf{r}}\in {{\mathbb{R}}}^{3},\ \hat{{\bf{s}}}\in {{\mathbb{S}}}^{2}.\end{array}\end{eqnarray*}$$

The Greenʼs function is obtained as [42]

$$\begin{eqnarray*}&&G({\bf{r}},\hat{{\bf{s}}};{{\bf{r}}}_{0},{\hat{{\bf{s}}}}_{0})=\displaystyle \frac{1}{{(2\pi )}^{2}}{\displaystyle \int }_{{{\mathbb{R}}}^{2}}{{\rm{e}}}^{-{\rm{i}}{\bf{q}}\cdot ({\boldsymbol{\rho }}-{{\boldsymbol{\rho }}}_{0})}\tilde{G}({\bf{q}};z,\hat{{\bf{s}}};{z}_{0},{\hat{{\bf{s}}}}_{0})\;{\rm{d}}{\bf{q}}.\end{eqnarray*}$$

Here,

$$\begin{eqnarray*}\begin{array}{rcl}\tilde{G}({\bf{q}};z,\hat{{\bf{s}}};{z}_{0},{\hat{{\bf{s}}}}_{0}) & = & \displaystyle \sum _{m=-L}^{L}\left\{\displaystyle \sum _{j=0}^{{M}^{m}-1}\displaystyle \frac{1}{2\pi {\hat{k}}_{z}({\nu }_{j}^{m}q){\mathcal{N}}({\nu }_{j}^{m})}\right.\\ & & \times {{\mathcal{R}}}_{\hat{{\bf{k}}}(\pm {\nu }_{j}^{m},{\bf{q}})}{\Phi }_{j\pm }^{m}(\hat{{\bf{s}}}){\Phi }_{j\pm }^{m*}({\hat{{\bf{s}}}}_{0}){{\rm{e}}}^{-{\hat{k}}_{z}({\nu }_{j}^{m}q)| z-{z}_{0}| /{\nu }_{j}^{m}}\\ & & +\left.{\displaystyle \int }_{0}^{1}\displaystyle \frac{1}{2\pi {\hat{k}}_{z}(\nu q){\mathcal{N}}(\nu )}{{\mathcal{R}}}_{\hat{{\bf{k}}}(\pm \nu ,{\bf{q}})}{\Phi }_{\pm \nu }^{m}(\hat{{\bf{s}}}){\Phi }_{\pm \nu }^{m*}({\hat{{\bf{s}}}}_{0}){{\rm{e}}}^{-{\hat{k}}_{z}(\nu q)| z-{z}_{0}| /\nu }\;{\rm{d}}\nu \right\},\end{array}\end{eqnarray*}$$

where upper signs are chosen for $z\gt {z}_{0}$ and lower signs are chosen for $z\lt {z}_{0}$. By letting $z\to {0}^{+}$ we obtain

$$\begin{eqnarray*}\begin{array}{rcl}I({\boldsymbol{\rho }},0,\hat{{\bf{s}}}) & = & {\displaystyle \int }_{{{\mathbb{S}}}^{2}}{\displaystyle \int }_{{{\mathbb{R}}}^{2}}G({\boldsymbol{\rho }},{0}^{+},\hat{{\bf{s}}};{{\boldsymbol{\rho }}}^{\prime },0,\hat{{\bf{s}}}\prime )\mu \prime I({{\boldsymbol{\rho }}}^{\prime },0,\hat{{\bf{s}}}\prime )\;d{\boldsymbol{\rho }}\prime {\rm{d}}\hat{{\bf{s}}}\prime \\ & & +{\displaystyle \int }_{{{\mathbb{S}}}^{2}}{\displaystyle \int }_{0}^{\infty }{\displaystyle \int }_{{{\mathbb{R}}}^{2}}G({\boldsymbol{\rho }},0,\hat{{\bf{s}}};{{\boldsymbol{\rho }}}^{\prime },z\prime ,\hat{{\bf{s}}}\prime )S({{\boldsymbol{\rho }}}^{\prime },z\prime ,\hat{{\bf{s}}}\prime )\;d{\boldsymbol{\rho }}\prime {\rm{d}}z\prime {\rm{d}}\hat{{\bf{s}}}\prime ,\end{array}\end{eqnarray*}$$

where $\hat{{\bf{s}}}\in {{\mathbb{S}}}^{2}$. We have

$$\begin{eqnarray}\begin{array}{rcl}\tilde{I}({\bf{q}},0,\hat{{\bf{s}}}) & = & {\displaystyle \int }_{{{\mathbb{S}}}^{2}}\tilde{G}({\bf{q}};{0}^{+},\hat{{\bf{s}}};0,\hat{{\bf{s}}}\prime ){\mu }^{\prime }\tilde{I}({\bf{q}},0,\hat{{\bf{s}}}\prime )\;{\rm{d}}\hat{{\bf{s}}}\prime \\ & & +{\displaystyle \int }_{{{\mathbb{S}}}^{2}}{\displaystyle \int }_{0}^{\infty }\tilde{G}({\bf{q}};0,\hat{{\bf{s}}};z\prime ,\hat{{\bf{s}}}\prime )\tilde{S}({\bf{q}},z\prime ,\hat{{\bf{s}}}\prime )\;{\rm{d}}z\prime {\rm{d}}\hat{{\bf{s}}}\prime .\end{array}\end{eqnarray} \tag{ 33 }$$

Definition 3.1. Let ${\xi }^{m}$ denote the positive eigenvalues, i.e., ${\xi }^{m}={\nu }_{j}^{m}$ ($j=0,1,...,{M}^{m}-1$) or ${\xi }^{m}=\nu \in (0,1)$. We drop the superscript and write $\xi ={\xi }^{m}$ if there is no confusion.

If we multiply (33) by $\mu {{\mathcal{R}}}_{\hat{{\bf{k}}}(-\xi ,{\bf{q}})}{\Phi }_{-\xi }^{{m\prime }^{*}}(\hat{{\bf{s}}})$ with some m' and $\xi ={\xi }^{{m}^{\prime }}\gt 0$, and integrate over ${{\mathbb{S}}}^{2}$, we obtain

$$\begin{eqnarray*}&&\begin{array}{l}{\displaystyle \int }_{{{\mathbb{S}}}^{2}}\mu \left({{\mathcal{R}}}_{\hat{{\bf{k}}}(-\xi ,{\bf{q}})}{\Phi }_{-\xi }^{{m\prime }^{*}}(\hat{{\boldsymbol{}}{\bf{s}}})\right)\tilde{I}({\bf{q}},0,\hat{{\bf{s}}})\;{\rm{d}}\hat{{\boldsymbol{}}{\bf{s}}}\\ \qquad ={\displaystyle \int }_{{{\mathbb{S}}}^{2}}{\displaystyle \int }_{0}^{\infty }\left({{\mathcal{R}}}_{\hat{{\bf{k}}}(-\xi ,{\bf{q}})}{\Phi }_{-\xi }^{{m\prime }^{*}}(\hat{{\boldsymbol{}}{\bf{s}}}\prime )\right){{\rm{e}}}^{-{\hat{k}}_{z}(\xi q)z\prime /\xi }\tilde{S}({\bf{q}},z\prime ,\hat{{\bf{s}}}\prime )\;{\rm{d}}z\prime {\rm{d}}\hat{{\boldsymbol{}}{\bf{s}}}\prime .\end{array}\end{eqnarray*}$$

Hence we can write the above equation as

$$\begin{eqnarray}&&\begin{array}{l}{\displaystyle \int }_{{{\mathbb{S}}}_{+}^{2}}\mu \left({{\mathcal{R}}}_{\hat{{\bf{k}}}(-\xi ,{\bf{q}})}{\Phi }_{-\xi }^{{m\prime }^{*}}(-\hat{{\bf{s}}})\right)\tilde{I}({\bf{q}},0,-\hat{{\bf{s}}})\;{\rm{d}}\hat{{\boldsymbol{}}{\bf{s}}}={\displaystyle \int }_{{{\mathbb{S}}}_{+}^{2}}\mu \left({{\mathcal{R}}}_{\hat{{\bf{k}}}(-\xi ,{\bf{q}})}{\Phi }_{-\xi }^{{m\prime }^{*}}(\hat{{\boldsymbol{}}{\bf{s}}})\right)\tilde{f}({\bf{q}},\hat{{\bf{s}}})\;{\rm{d}}\hat{{\boldsymbol{}}{\bf{s}}}\\ \qquad -{\displaystyle \int }_{{{\mathbb{S}}}^{2}}{\displaystyle \int }_{0}^{\infty }\left({{\mathcal{R}}}_{\hat{{\bf{k}}}(-\xi ,{\bf{q}})}{\Phi }_{-\xi }^{{m\prime }^{*}}(\hat{{\boldsymbol{}}{\bf{s}}}\prime )\right){{\rm{e}}}^{-{\hat{k}}_{z}(\xi q)z\prime /\xi }\tilde{S}({\bf{q}},z\prime ,\hat{{\boldsymbol{}}{\bf{s}}}\prime )\;{\rm{d}}z\prime {\rm{d}}\hat{{\boldsymbol{}}{\bf{s}}}\prime .\end{array}\end{eqnarray} \tag{ 34 }$$

By the expansion in (3), we obtain the following key F_N equation

$$\begin{eqnarray}&&\displaystyle \sum _{m=-{l}_{\mathrm{max}}}^{{l}_{\mathrm{max}}}\displaystyle \sum _{l=| m| ,| m| +2,...}{A}_{{lm}}^{{m}^{\prime }}(\xi ,{\bf{q}}){c}_{{lm}}({\bf{q}})={K}^{{m}^{\prime }}(\xi ,{\bf{q}}),\end{eqnarray} \tag{ 35 }$$

where $-L\leqslant m\prime \leqslant L$. Here,

$$\begin{eqnarray*}\begin{array}{rcl}{A}_{{lm}}^{{m}^{\prime }}(\xi ,{\bf{q}}) & = & {\displaystyle \int }_{{{\mathbb{S}}}_{+}^{2}}\mu {Y}_{{lm}}(\hat{{\bf{s}}}){{\mathcal{R}}}_{\hat{{\bf{k}}}(-\xi ,{\bf{q}})}{\Phi }_{-\xi }^{{m\prime }^{*}}(-\hat{{\bf{s}}})\;{\rm{d}}\hat{{\bf{s}}},\\ {K}^{{m}^{\prime }}(\xi ,{\bf{q}}) & = & {\displaystyle \int }_{{{\mathbb{S}}}_{+}^{2}}\mu \tilde{f}({\bf{q}},\hat{{\bf{s}}}){{\mathcal{R}}}_{\hat{{\bf{k}}}(-\xi ,{\bf{q}})}{\Phi }_{-\xi }^{{m\prime }^{*}}(\hat{{\bf{s}}})\;{\rm{d}}\hat{{\bf{s}}}\\ & & -{\displaystyle \int }_{{{\mathbb{S}}}^{2}}{\displaystyle \int }_{0}^{\infty }{{\rm{e}}}^{-{\hat{k}}_{z}(\xi q)z\prime /\xi }\tilde{S}({\bf{q}},z\prime ,\hat{{\bf{s}}}\prime ){{\mathcal{R}}}_{\hat{{\bf{k}}}(-\xi ,{\bf{q}})}{\Phi }_{-\xi }^{{m\prime }^{*}}(\hat{{\bf{s}}}\prime )\;{\rm{d}}z\prime {\rm{d}}\hat{{\bf{s}}}\prime .\end{array}\end{eqnarray*}$$

Remark 3.2. In the above proof we used the Greenʼs function in the free space to derive (34). This approach is similar to the C_N method [2, 23]. If the Greenʼs function for the half space is used, we can explicitly give $\tilde{I}({\bf{q}},0,-\hat{{\bf{s}}})$ without relying on (3) and (35) [52]. However, the half-space Greenʼs function in three dimensions is not yet known.

We obtain

$$\begin{eqnarray*}&&{A}_{{lm}}^{{m}^{\prime }}(\xi ,{\bf{q}})={A}_{{lm}}^{{m}^{\prime }}(\xi ,q){{\rm{e}}}^{{\rm{i}}m{\varphi }_{{\bf{q}}}},\end{eqnarray*}$$

where

$$\begin{eqnarray*}\begin{array}{rcl}{A}_{{lm}}^{{m}^{\prime }}(\xi ,q) & = & {(-1)}^{m}{\hat{k}}_{z}(\xi q)\sqrt{\displaystyle \frac{\pi }{2l+1}}{d}_{{{mm}}^{\prime }}^{l}[{\rm{i}}\tau (\xi q)]\left(\sqrt{{(l+1)}^{2}-m{\prime }^{2}}{g}_{l+1}^{{m}^{\prime }}(\xi )\right.\\ & & +\left.\sqrt{{l}^{2}-m{\prime }^{2}}{g}_{l-1}^{{m}^{\prime }}(\xi )\right)\\ & & -{\rm{i}}\displaystyle \frac{| \xi q| }{2}\sqrt{\displaystyle \frac{\pi }{2l+1}}{(-1)}^{m}\displaystyle \sum _{m\prime \prime =-l}^{l}{d}_{{{mm}}^{\prime\prime }}^{l}[i\tau (\xi q)]\\ & & \times \left[{\delta }_{m\prime \prime ,m\prime -1}\left(\sqrt{(l-m\prime \prime )(l-m\prime )}{g}_{l-1}^{{m}^{\prime }}(\xi )-\sqrt{(l+m\prime +1)(l+m\prime )}{g}_{l+1}^{{m}^{\prime }}(\xi )\right)\right.\\ & & +\left.{\delta }_{m\prime \prime ,m\prime +1}\left(\sqrt{(l-m\prime +1)(l-m\prime )}{g}_{l+1}^{{m}^{\prime }}(\xi )-\sqrt{(l+m\prime \prime )(l+m\prime )}{g}_{l-1}^{{m}^{\prime }}(\xi )\right)\right]\\ & & +\displaystyle \frac{\varpi \xi }{2}{(-1)}^{l}\sqrt{\displaystyle \frac{2l+1}{4\pi }\displaystyle \frac{(l-m)!}{(l+m)!}}{[\mathrm{sgn}(m\prime )]}^{{m}^{\prime }}\displaystyle \frac{\sqrt{(2| m\prime | )!}}{(2| m\prime | -1)!!}\end{array}\end{eqnarray*}$$

$$\begin{eqnarray}\begin{array}{rcl} & & \times \displaystyle \sum _{m\prime \prime =-| m\prime | }^{| m\prime | }{(-1)}^{{m}^{\prime\prime }}\sqrt{\displaystyle \frac{(| m\prime | -m\prime \prime )!}{(| m\prime | +m\prime \prime )!}}\\ & & \times {d}_{m\prime \prime ,-m\prime }^{| m\prime | }[{\rm{i}}\tau (\xi q)]{\displaystyle \int }_{{{\mathbb{S}}}_{+}^{2}}\displaystyle \frac{{g}^{{m}^{\prime }}\left(-\xi ,{\hat{k}}_{z}(\xi q)\mu -{\rm{i}}\xi q\sqrt{1-{\mu }^{2}}\mathrm{cos}\varphi \right)}{\xi +{\hat{k}}_{z}(\xi q)\mu -{\rm{i}}\xi q\sqrt{1-{\mu }^{2}}\mathrm{cos}\varphi }\\ & & \times \mu {P}_{| m\prime | }^{{m}^{\prime\prime }}(\mu ){P}_{l}^{m}(\mu ){{\rm{e}}}^{{\rm{i}}(m+m\prime \prime )\varphi }\;{\rm{d}}\hat{{\bf{s}}}.\end{array}\end{eqnarray} \tag{ 36 }$$

If ${K}^{{m}^{\prime }}(\xi ,{\bf{q}})$ is independent of ${\varphi }_{{\bf{q}}}$ and ${K}^{-m\prime }={K}^{{m}^{\prime }}$, then

$$\begin{eqnarray*}&&{c}_{{lm}}({\bf{q}})={c}_{{lm}}(q){{\rm{e}}}^{-{\rm{i}}m{\varphi }_{{\bf{q}}}},\qquad {c}_{l,-m}(q)={(-1)}^{m}{c}_{{lm}}(q).\end{eqnarray*}$$

Here the coefficients ${c}_{{lm}}(q)$ are solutions to

$$\begin{eqnarray}&&\begin{array}{l}\displaystyle \sum _{m=0}^{{l}_{\mathrm{max}}}\displaystyle \sum _{\alpha =0}^{\left\lfloor ({l}_{\mathrm{max}}-m)/2\right\rfloor}\left[{A}_{m+2\alpha ,m}^{{m}^{\prime }}(\xi ,q)+(1-{\delta }_{m0}){(-1)}^{m}{A}_{m+2\alpha ,-m}^{{m}^{\prime }}(\xi ,q)\right]{c}_{m+2\alpha ,m}(q)\\ \qquad ={K}^{{m}^{\prime }}(\xi ,{\bf{q}}),\end{array}\end{eqnarray} \tag{ 37 }$$

where ${A}_{m+2\alpha ,m}^{{m}^{\prime }}(\xi ,q)$ are given in (36).

The rest of the section is devoted to the calculations of (36) and (37).

First, ${A}_{{lm}}^{{m}^{\prime }}(\xi ,{\bf{q}})$ are computed as follow. We begin by noting that

$$\begin{eqnarray}\begin{array}{rcl}{A}_{{lm}}^{{m}^{\prime }}(\xi ,{\bf{q}}) & = & {\displaystyle \int }_{{{\mathbb{S}}}^{2}}\mu {Y}_{{lm}}(\hat{{\bf{s}}}){{\mathcal{R}}}_{\hat{{\bf{k}}}(-\xi ,{\bf{q}})}{\Phi }_{-\xi }^{{m\prime }^{*}}(-\hat{{\bf{s}}})\;{\rm{d}}\hat{{\bf{s}}}-{\displaystyle \int }_{{{\mathbb{S}}}_{-}^{2}}\mu {Y}_{{lm}}(\hat{{\bf{s}}}){{\mathcal{R}}}_{\hat{{\bf{k}}}(-\xi ,{\bf{q}})}{\Phi }_{-\xi }^{{m\prime }^{*}}(-\hat{{\bf{s}}})\;{\rm{d}}\hat{{\bf{s}}}\\ & = & {\displaystyle \int }_{{{\mathbb{S}}}^{2}}\left({{\mathcal{R}}}_{\hat{{\bf{k}}}(-\xi ,{\bf{q}})}^{-1}\mu {Y}_{{lm}}(\hat{{\bf{s}}})\right){\Phi }_{-\xi }^{{m\prime }^{*}}(-\hat{{\bf{s}}})\;{\rm{d}}\hat{{\bf{s}}}+{\displaystyle \int }_{{{\mathbb{S}}}_{+}^{2}}\mu {Y}_{{lm}}(-\hat{{\bf{s}}}){{\mathcal{R}}}_{\hat{{\bf{k}}}(-\xi ,{\bf{q}})}{\Phi }_{-\xi }^{{m\prime }^{*}}(\hat{{\bf{s}}})\;{\rm{d}}\hat{{\bf{s}}}.\end{array}\end{eqnarray} \tag{ 38 }$$

We obtain the first term on the right-hand side of (38) as

$$\begin{eqnarray*}\begin{array}{rcl}{\rm{(1st}}\;{\rm{term)}} & = & {\displaystyle \int }_{{{\mathbb{S}}}^{2}}{\Phi }_{-\xi }^{{m\prime }^{*}}(-\hat{{\bf{s}}}){{\mathcal{R}}}_{\hat{{\bf{k}}}(-\xi ,{\bf{q}})}^{-1}\mu {Y}_{{lm}}(\hat{{\bf{s}}})\;{\rm{d}}\hat{{\bf{s}}}\\ & = & \displaystyle \sum _{m\prime \prime =-l}^{l}{{\rm{e}}}^{{\rm{i}}m{\varphi }_{\hat{{\bf{k}}}}}{d}_{{{mm}}^{\prime\prime }}^{l}({\theta }_{\hat{{\bf{k}}}}){\displaystyle \int }_{{{\mathbb{S}}}^{2}}{\Phi }_{-\xi }^{{m\prime }^{*}}(-\hat{{\bf{s}}})\left({{\mathcal{R}}}_{\hat{{\bf{k}}}(-\xi ,{\bf{q}})}^{-1}\mu \right){Y}_{{{lm}}^{\prime\prime }}(\hat{{\bf{s}}})\;{\rm{d}}\hat{{\bf{s}}}\\ & = & {(-1)}^{l}\displaystyle \sum _{m\prime \prime =-l}^{l}{{\rm{e}}}^{{\rm{i}}m{\varphi }_{\hat{{\bf{k}}}}}{d}_{{{mm}}^{\prime\prime }}^{l}({\theta }_{\hat{{\bf{k}}}}){\displaystyle \int }_{{{\mathbb{S}}}^{2}}{\Phi }_{-\xi }^{{m\prime }^{*}}(\hat{{\bf{s}}})\left(-{\hat{k}}_{z}(\xi q)\mu \right.\\ & & +\left.{\rm{i}}| \xi q| \sqrt{1-{\mu }^{2}}\mathrm{cos}\varphi \right){Y}_{{{lm}}^{\prime\prime }}(\hat{{\bf{s}}})\;{\rm{d}}\hat{{\bf{s}}}.\end{array}\end{eqnarray*}$$

Here,

$$\begin{eqnarray*}&&\begin{array}{l}{\displaystyle \int }_{{{\mathbb{S}}}^{2}}{\Phi }_{-\xi }^{{m\prime }^{*}}(\hat{{\bf{s}}})\mu {Y}_{{{lm}}^{\prime\prime }}(\hat{{\bf{s}}})\;{\rm{d}}\hat{{\bf{s}}}\\ \qquad =\sqrt{\displaystyle \frac{{(l+1)}^{2}-m\prime {\prime }^{2}}{4{(l+1)}^{2}-1}}{\displaystyle \int }_{{{\mathbb{S}}}^{2}}{\Phi }_{-\xi }^{{m\prime }^{*}}(\hat{{\bf{s}}}){Y}_{l+1,m\prime \prime }(\hat{{\bf{s}}})\;{\rm{d}}\hat{{\bf{s}}}\\ \qquad \;\;\;+\sqrt{\displaystyle \frac{{l}^{2}-m\prime {\prime }^{2}}{4{l}^{2}-1}}{\displaystyle \int }_{{{\mathbb{S}}}^{2}}{\Phi }_{-\xi }^{{m\prime }^{*}}(\hat{{\bf{s}}}){Y}_{l-1,m\prime \prime }(\hat{{\bf{s}}})\;{\rm{d}}\hat{{\bf{s}}}\end{array}\end{eqnarray*}$$

$$\begin{eqnarray*}&&\begin{array}{l}\qquad ={\delta }_{m\prime m\prime \prime }{(-1)}^{{m}^{\prime }}\sqrt{\displaystyle \frac{\pi }{2l+1}}\left(\sqrt{{(l+1)}^{2}-m{\prime }^{2}}{g}_{l+1}^{{m}^{\prime }}(-\xi )+\sqrt{{l}^{2}-m{\prime }^{2}}{g}_{l-1}^{{m}^{\prime }}(-\xi )\right),\end{array}\end{eqnarray*}$$

where we used $\mu {P}_{l}^{{m}^{\prime }}(\mu )=[(l+m\prime ){P}_{l-1}^{{m}^{\prime }}(\mu )+(l-m\prime +1){P}_{l+1}^{{m}^{\prime }}]/(2l+1)$. We also have

$$\begin{eqnarray*}&&\begin{array}{l}{\displaystyle \int }_{{{\mathbb{S}}}^{2}}{\Phi }_{-\xi }^{{m\prime }^{*}}(\hat{{\bf{s}}})\sqrt{1-{\mu }^{2}}\mathrm{cos}\varphi {Y}_{{{lm}}^{\prime\prime }}(\hat{{\bf{s}}})\;{\rm{d}}\hat{{\bf{s}}}=\sqrt{\displaystyle \frac{(2l+1)\pi }{4}\displaystyle \frac{(l-m\prime \prime )!}{(l+m\prime \prime )!}}\left({\delta }_{m\prime \prime ,m\prime -1}+{\delta }_{m\prime \prime ,m\prime +1}\right)\\ \qquad \times {\displaystyle \int }_{-1}^{1}{\phi }^{{m}^{\prime }}(-\xi ,\mu ){(1-{\mu }^{2})}^{(| m\prime | +1)/2}{P}_{l}^{{m}^{\prime\prime }}(\mu )\;{\rm{d}}\mu .\end{array}\end{eqnarray*}$$

Using $\sqrt{1-{\mu }^{2}}{P}_{l}^{m\prime -1}(\mu )=[{P}_{l-1}^{{m}^{\prime }}(\mu )-{P}_{l+1}^{{m}^{\prime }}(\mu )]/(2l+1)$, $\sqrt{1-{\mu }^{2}}{P}_{l}^{m\prime +1}(\mu )=(l-m\prime )\mu {P}_{l}^{{m}^{\prime }}(\mu )-(l+m\prime ){P}_{l-1}^{{m}^{\prime }}(\mu )$, we obtain

$$\begin{eqnarray*}&&\begin{array}{l}{\displaystyle \int }_{{{\mathbb{S}}}^{2}}{\Phi }_{-\xi }^{{m\prime }^{*}}(\hat{{\bf{s}}})\sqrt{1-{\mu }^{2}}\mathrm{cos}\varphi {Y}_{{{lm}}^{\prime\prime }}(\hat{{\bf{s}}})\;{\rm{d}}\hat{{\bf{s}}}\\ \qquad =\displaystyle \frac{1}{2}\sqrt{\displaystyle \frac{\pi }{2l+1}}{(-1)}^{l+1}\left[{\delta }_{m\prime \prime ,m\prime -1}\left(\sqrt{(l-m\prime +1)(l-m\prime )}{g}_{l-1}^{{m}^{\prime }}(\xi )\right.\right.\\ \qquad -\left.\sqrt{(l+m\prime +1)(l+m\prime )}{g}_{l+1}^{{m}^{\prime }}(\xi )\right)\\ \qquad +\left.{\delta }_{m\prime \prime ,m\prime +1}\left(\sqrt{(l-m\prime +1)(l-m\prime )}{g}_{l+1}^{{m}^{\prime }}(\xi )-\sqrt{(l+m\prime +1)(l+m\prime )}{g}_{l-1}^{{m}^{\prime }}(\xi )\right)\right].\end{array}\end{eqnarray*}$$

Therefore,

$$\begin{eqnarray*}&&\begin{array}{l}{\rm{(1st}}\;{\rm{term)}}\\ \qquad =-{(-1)}^{l+m\prime }{\hat{k}}_{z}(\xi q)\sqrt{\displaystyle \frac{\pi }{2l+1}}{{\rm{e}}}^{{\rm{i}}m{\varphi }_{\hat{{\bf{k}}}}}{d}_{{{mm}}^{\prime }}^{l}({\theta }_{\hat{{\bf{k}}}})\left(\sqrt{{(l+1)}^{2}-m{\prime }^{2}}{g}_{l+1}^{{m}^{\prime }}(-\xi )\right.\\ \qquad +\left.\sqrt{{l}^{2}-m{\prime }^{2}}{g}_{l-1}^{{m}^{\prime }}(-\xi )\right)-{\rm{i}}\displaystyle \frac{| \xi q| }{2}\sqrt{\displaystyle \frac{\pi }{2l+1}}\displaystyle \sum _{m\prime \prime =-l}^{l}{{\rm{e}}}^{{\rm{i}}m{\varphi }_{\hat{{\bf{k}}}}}{d}_{{{mm}}^{\prime\prime }}^{l}({\theta }_{\hat{{\bf{k}}}})\\ \qquad \times \left[{\delta }_{m\prime \prime ,m\prime -1}\left(\sqrt{(l-m\prime +1)(l-m\prime )}{g}_{l-1}^{{m}^{\prime }}(\xi )-\sqrt{(l+m\prime +1)(l+m\prime )}{g}_{l+1}^{{m}^{\prime }}(\xi )\right)\right.\\ \qquad +\left.{\delta }_{m\prime \prime ,m\prime +1}\left(\sqrt{(l-m\prime +1)(l-m\prime )}{g}_{l+1}^{{m}^{\prime }}(\xi )-\sqrt{(l+m\prime +1)(l+m\prime )}{g}_{l-1}^{{m}^{\prime }}(\xi )\right)\right].\end{array}\end{eqnarray*}$$

We will use

$$\begin{eqnarray*}\begin{array}{rcl}{(1-{\mu }^{2})}^{| m\prime | /2}{{\rm{e}}}^{-{\rm{i}}m\prime \varphi } & = & \displaystyle \frac{{(-1)}^{| m\prime | }}{(2| m\prime | -1)!!}{P}_{| m\prime | }^{| m\prime | }(\mu ){{\rm{e}}}^{-{\rm{i}}m\prime \varphi }\\ & = & {[\mathrm{sgn}(m\prime )]}^{{m}^{\prime }}\displaystyle \frac{\sqrt{4\pi (2| m\prime | +1)!}}{(2| m\prime | +1)!!}{Y}_{| m\prime | ,-m\prime }(\hat{{\bf{s}}}).\end{array}\end{eqnarray*}$$

The second term on the right-hand side of (38) is calculated as

$$\begin{eqnarray*}\begin{array}{rcl}{\rm{(2nd}}\;{\rm{term)}} & = & {\displaystyle \int }_{{{\mathbb{S}}}_{+}^{2}}\left({{\mathcal{R}}}_{\hat{{\bf{k}}}(-\xi ,{\bf{q}})}{\Phi }_{-\xi }^{{m\prime }^{*}}(\hat{{\bf{s}}})\right)\mu {Y}_{{lm}}(-\hat{{\bf{s}}})\;{\rm{d}}\hat{{\bf{s}}}\\ & = & \displaystyle \frac{\varpi \xi }{2}{\displaystyle \int }_{{{\mathbb{S}}}_{+}^{2}}\displaystyle \frac{{g}^{{m}^{\prime }}\left(-\xi ,{\hat{k}}_{z}(\xi q)\mu -{\rm{i}}\xi q\sqrt{1-{\mu }^{2}}\mathrm{cos}(\varphi -{\varphi }_{{\bf{q}}})\right)}{\xi +{\hat{k}}_{z}(\xi q)\mu -{\rm{i}}\xi q\sqrt{1-{\mu }^{2}}\mathrm{cos}(\varphi -{\varphi }_{{\bf{q}}})}\end{array}\end{eqnarray*}$$

$$\begin{eqnarray*}\begin{array}{rcl} & & \times {[\mathrm{sgn}(m\prime )]}^{{m}^{\prime }}\displaystyle \frac{\sqrt{4\pi (2| m\prime | +1)!}}{(2| m\prime | +1)!!}\\ & & \times \displaystyle \sum _{m\prime \prime =-| m\prime | }^{| m\prime | }{{\rm{e}}}^{-{\rm{i}}m\prime \prime {\varphi }_{\hat{{\bf{k}}}}}{d}_{m\prime \prime ,-m\prime }^{| m\prime | }({\theta }_{\hat{{\bf{k}}}})\mu {Y}_{| m\prime | m\prime \prime }(\hat{{\bf{s}}}){Y}_{{lm}}(-\hat{{\bf{s}}})\;{\rm{d}}\hat{{\bf{s}}}\\ & = & \displaystyle \frac{\varpi \xi }{2}{(-1)}^{l}\sqrt{\displaystyle \frac{2l+1}{4\pi }\displaystyle \frac{(l-m)!}{(l+m)!}}{[\mathrm{sgn}(m\prime )]}^{{m}^{\prime }}\displaystyle \frac{\sqrt{(2| m\prime | )!}}{(2| m\prime | -1)!!}\\ & & \times \displaystyle \sum _{m\prime \prime =-| m\prime | }^{| m\prime | }\sqrt{\displaystyle \frac{(| m\prime | -m\prime \prime )!}{(| m\prime | +m\prime \prime )!}}{{\rm{e}}}^{-{\rm{i}}m\prime \prime {\varphi }_{\hat{{\bf{k}}}}}{d}_{m\prime \prime ,-m\prime }^{| m\prime | }({\theta }_{\hat{{\bf{k}}}})\\ & & \times {\displaystyle \int }_{{{\mathbb{S}}}_{+}^{2}}\left({{\mathcal{R}}}_{\hat{{\bf{k}}}(-\xi ,{\bf{q}})}\displaystyle \frac{{g}^{{m}^{\prime }}(-\xi ,\mu )}{\xi +\mu }\right)\mu {P}_{| m\prime | }^{{m}^{\prime\prime }}(\mu ){P}_{l}^{m}(\mu ){{\rm{e}}}^{{\rm{i}}(m+m\prime \prime )\varphi }\;{\rm{d}}\hat{{\bf{s}}}.\end{array}\end{eqnarray*}$$

We note that the relation ${g}^{m}(-\xi ,\mu )={g}^{m}(\xi ,-\mu )$ implies ${{\mathcal{R}}}_{\hat{{\bf{k}}}}{\phi }^{m}(-\xi ,\mu )={{\mathcal{R}}}_{\hat{{\bf{k}}}}{\phi }^{m}(\xi ,-\mu )$ for a fixed $\hat{{\bf{k}}}$. Thus (36) is obtained.

Next, (37) is obtained as follows. Using ${p}_{l}^{-m}(\mu )={(-1)}^{m}{p}_{l}^{m}(\mu )$, ${g}_{l}^{-m}(\xi )={(-1)}^{m}{g}_{l}^{m}(\xi )$, and ${g}^{m}(\xi ,\mu )={g}^{-m}(\xi ,\mu )$, we can show that

$$\begin{eqnarray*}&&{A}_{l,-m}^{-m\prime }(\xi ,{\bf{q}}){{\rm{e}}}^{{\rm{i}}m{\varphi }_{{\bf{q}}}}={(-1)}^{m}{A}_{{lm}}^{{m}^{\prime }}(\xi ,{\bf{q}}){{\rm{e}}}^{-{\rm{i}}m{\varphi }_{{\bf{q}}}}.\end{eqnarray*}$$

Since we assume ${K}^{-m\prime }(\xi ,{\bf{q}})={K}^{{m}^{\prime }}(\xi ,{\bf{q}})$, we have

$$\begin{eqnarray*}\begin{array}{rcl}\displaystyle \sum _{m=-{l}_{\mathrm{max}}}^{{l}_{\mathrm{max}}}\displaystyle \sum _{l}{A}_{{lm}}^{{m}^{\prime }}(\xi ,{\bf{q}}){c}_{{lm}}({\bf{q}}) & = & \displaystyle \sum _{m=-{l}_{\mathrm{max}}}^{{l}_{\mathrm{max}}}\displaystyle \sum _{l}{A}_{{lm}}^{-m\prime }(\xi ,{\bf{q}}){c}_{{lm}}({\bf{q}})\\ & = & \displaystyle \sum _{m=-{l}_{\mathrm{max}}}^{{l}_{\mathrm{max}}}\displaystyle \sum _{l}{A}_{l,-m}^{-m\prime }(\xi ,{\bf{q}}){c}_{l,-m}({\bf{q}})\\ & = & \displaystyle \sum _{m=-{l}_{\mathrm{max}}}^{{l}_{\mathrm{max}}}\displaystyle \sum _{l}{A}_{{lm}}^{{m}^{\prime }}(\xi ,{\bf{q}}){(-1)}^{m}{{\rm{e}}}^{-2{\rm{i}}m{\varphi }_{{\bf{q}}}}{c}_{l,-m}({\bf{q}}).\end{array}\end{eqnarray*}$$

This implies

$$\begin{eqnarray*}&&{c}_{l,-m}({\bf{q}})={(-1)}^{m}{{\rm{e}}}^{2{\rm{i}}m{\varphi }_{{\bf{q}}}}{c}_{{lm}}({\bf{q}}).\end{eqnarray*}$$

Moreover since we assume that ${K}^{{m}^{\prime }}(\xi ,{\bf{q}})$ is independent of ${\varphi }_{{\bf{q}}}$, we have

$$\begin{eqnarray*}&&\displaystyle \sum _{m=-{l}_{\mathrm{max}}}^{{l}_{\mathrm{max}}}\displaystyle \sum _{l}{A}_{{lm}}^{{m}^{\prime }}(\xi ,{\bf{q}}){c}_{{lm}}({\bf{q}})=\displaystyle \sum _{m=-{l}_{\mathrm{max}}}^{{l}_{\mathrm{max}}}\displaystyle \sum _{l}{A}_{{lm}}^{{m}^{\prime }}(\xi ,q){{\rm{e}}}^{{\rm{i}}m{\varphi }_{{\bf{q}}}}{c}_{{lm}}({\bf{q}}).\end{eqnarray*}$$

This implies

$$\begin{eqnarray*}&&{c}_{{lm}}({\bf{q}})={c}_{{lm}}(q){{\rm{e}}}^{-{\rm{i}}m{\varphi }_{{\bf{q}}}}.\end{eqnarray*}$$

Therefore we obtain

$$\begin{eqnarray*}&&{c}_{l,-m}(q)={(-1)}^{m}{c}_{{lm}}(q).\end{eqnarray*}$$

By using this relation in (35), we obtain (37).

Let us consider a structured illumination in the half space:

$$\begin{eqnarray*}\left\{\begin{array}{ll}\hat{{\bf{s}}}\cdot \nabla I({\bf{r}},\hat{{\bf{s}}})+I({\bf{r}},\hat{{\bf{s}}})=\varpi {\displaystyle \int }_{{{\mathbb{S}}}^{2}}p(\hat{{\bf{s}}},\hat{{\bf{s}}}\prime )I({\bf{r}},\hat{{\bf{s}}}\prime )\;{\rm{d}}\hat{{\bf{s}}}\prime , & z\gt 0,\\ I({\bf{r}},\hat{{\bf{s}}})=f({\boldsymbol{\rho }},\hat{{\bf{s}}}), & z=0,\mu \in (0,1],\\ I({\bf{r}},\hat{{\bf{s}}})\to 0, & z\to \infty .\end{array}\right.\end{eqnarray*}$$

Here the incoming boundary value f is given by

$$\begin{eqnarray*}&&f({\boldsymbol{\rho }},\hat{{\bf{s}}})={I}_{0}\left[1+{A}_{0}\mathrm{cos}({{\bf{q}}}_{0}\cdot {\boldsymbol{\rho }}+{B}_{0})\right]\delta (\hat{{\bf{s}}}-{\hat{{\bf{s}}}}_{0}),\qquad {\hat{{\bf{s}}}}_{0}\in {{\mathbb{S}}}_{+}^{2},\end{eqnarray*}$$

where I₀ is the amplitude, A₀ is the modulation depth, and B₀ is the phase of the source. It is enough if we consider [40]

$$\begin{eqnarray}&&f({\boldsymbol{\rho }},\hat{{\bf{s}}})={{\rm{e}}}^{-{\rm{i}}{{\bf{q}}}_{0}\cdot {\boldsymbol{\rho }}}\delta (\hat{{\bf{s}}}-{\hat{{\bf{s}}}}_{0}),\quad {\hat{{\bf{s}}}}_{0}\in {{\mathbb{S}}}_{+}^{2},\end{eqnarray} \tag{ 39 }$$

where ${\hat{{\bf{s}}}}_{0}$ has the azimuthal angle ${\varphi }_{0}$ and the cosine of the polar angle ${\mu }_{0}$. By collision expansion we can write I as

$$\begin{eqnarray*}&&I({\bf{r}},\hat{{\bf{s}}})={I}_{b}({\bf{r}},\hat{{\bf{s}}})+{I}_{s}({\bf{r}},\hat{{\bf{s}}}),\end{eqnarray*}$$

where I_b is the ballistic term and I_s is the scattered part. They satisfy

$$\begin{eqnarray*}\left\{\begin{array}{ll}\hat{{\bf{s}}}\cdot \nabla {I}_{b}({\bf{r}},\hat{{\bf{s}}})+{I}_{b}({\bf{r}},\hat{{\bf{s}}})=0, & z\gt 0,\\ {I}_{b}({\bf{r}},\hat{{\bf{s}}})=f({\boldsymbol{\rho }},\hat{{\bf{s}}}), & z=0,\mu \in (0,1],\end{array}\right.\end{eqnarray*}$$

and

$$\begin{eqnarray*}\left\{\begin{array}{ll}\hat{{\bf{s}}}\cdot \nabla {I}_{s}({\bf{r}},\hat{{\bf{s}}})+{I}_{s}({\bf{r}},\hat{{\bf{s}}})=\varpi {\displaystyle \int }_{{{\mathbb{S}}}^{2}}p(\hat{{\bf{s}}},\hat{{\bf{s}}}\prime ){I}_{s}({\bf{r}},\hat{{\bf{s}}}\prime )\;{\rm{d}}\hat{{\bf{s}}}\prime +S({\bf{r}},\hat{{\bf{s}}}), & z\gt 0,\\ {I}_{s}({\bf{r}},\hat{{\bf{s}}})=0, & z=0,\mu \in (0,1],\end{array}\right.\end{eqnarray*}$$

where

$$\begin{eqnarray*}&&S({\bf{r}},\hat{{\bf{s}}})=\varpi {\displaystyle \int }_{{{\mathbb{S}}}^{2}}p(\hat{{\bf{s}}},\hat{{\bf{s}}}\prime ){I}_{b}({\bf{r}},\hat{{\bf{s}}}\prime )\;{\rm{d}}\hat{{\bf{s}}}\prime .\end{eqnarray*}$$

We also assume ${I}_{b},{I}_{s}\to 0$ as $z\to \infty$. Let us put

$$\begin{eqnarray*}&&{\hat{{\bf{s}}}}_{0}=\hat{{\bf{z}}}.\end{eqnarray*}$$

We obtain

$$\begin{eqnarray*}&&{I}_{{\rm{b}}}({\bf{r}},\hat{{\bf{s}}})={{\rm{e}}}^{-{\rm{i}}{{\bf{q}}}_{0}\cdot {\boldsymbol{\rho }}}{{\rm{e}}}^{-z}\delta (\hat{{\bf{s}}}-\hat{{\bf{z}}}).\end{eqnarray*}$$

We have

$$\begin{eqnarray*}&&S({\bf{r}},\hat{{\bf{s}}})=\displaystyle \frac{\varpi }{4\pi }{{\rm{e}}}^{-{\rm{i}}{{\bf{q}}}_{0}\cdot {\boldsymbol{\rho }}}{{\rm{e}}}^{-z}\displaystyle \sum _{l=0}^{L}{\beta }_{l}{P}_{l}(\mu ),\qquad \tilde{S}({\bf{q}},z,\hat{{\bf{s}}})=\pi \varpi \delta ({\bf{q}}-{{\bf{q}}}_{0}){{\rm{e}}}^{-z}\displaystyle \sum _{l=0}^{L}{\beta }_{l}{P}_{l}(\mu ).\end{eqnarray*}$$

Furthermore we assume that ${{\bf{q}}}_{0}$ is parallel to the x-axis:

$$\begin{eqnarray}&&{{\bf{q}}}_{0}={q}_{0}\hat{{\bf{x}}}.\end{eqnarray} \tag{ 40 }$$

We obtain

$$\begin{eqnarray*}\begin{array}{rcl}{K}^{{m}^{\prime }}(\xi ,{\bf{q}}) & = & -{\displaystyle \int }_{{{\mathbb{S}}}^{2}}{\displaystyle \int }_{0}^{\infty }{{\rm{e}}}^{-{\hat{k}}_{z}(\xi q)z\prime /\xi }\tilde{S}({\bf{q}},z\prime ,\hat{{\bf{s}}}\prime ){{\mathcal{R}}}_{\hat{{\bf{k}}}(-\xi ,{\bf{q}})}{\Phi }_{-\xi }^{{m\prime }^{*}}(\hat{{\bf{s}}}\prime )\;{\rm{d}}z\prime {\rm{d}}\hat{{\bf{s}}}\prime \\ & = & \displaystyle \frac{-2{\pi }^{3/2}\varpi \xi }{\xi +{\hat{k}}_{z}(\xi q)}\delta ({\bf{q}}-{{\bf{q}}}_{0})\displaystyle \sum _{l=0}^{L}\displaystyle \frac{{\beta }_{l}}{\sqrt{2l+1}}{\displaystyle \int }_{{{\mathbb{S}}}^{2}}{Y}_{l0}(\hat{{\bf{s}}}\prime ){{\mathcal{R}}}_{\hat{{\bf{k}}}(-\xi ,{\bf{q}})}{\Phi }_{-\xi }^{{m\prime }^{*}}(\hat{{\bf{s}}}\prime )\;{\rm{d}}\hat{{\bf{s}}}\prime .\end{array}\end{eqnarray*}$$

We note that

$$\begin{eqnarray*}\begin{array}{rcl}{\displaystyle \int }_{{{\mathbb{S}}}^{2}}{Y}_{l0}(\hat{{\bf{s}}}){{\mathcal{R}}}_{\hat{{\bf{k}}}(-\xi ,{\bf{q}})}{\Phi }_{-\xi }^{{m\prime }^{*}}(\hat{{\bf{s}}})\;{\rm{d}}\hat{{\bf{s}}} & = & {\displaystyle \int }_{{{\mathbb{S}}}^{2}}\left({{\mathcal{R}}}_{\hat{{\bf{k}}}(-\xi ,{\bf{q}})}^{-1}{Y}_{l0}(\hat{{\bf{s}}})\right){\Phi }_{-\xi }^{{m\prime }^{*}}(\hat{{\bf{s}}})\;{\rm{d}}\hat{{\bf{s}}}\\ & = & \displaystyle \sum _{m\prime \prime =-l}^{l}{d}_{0{m}^{\prime\prime }}^{l}({\theta }_{\hat{{\bf{k}}}}){\displaystyle \int }_{{{\mathbb{S}}}^{2}}{\Phi }_{-\xi }^{{m\prime }^{*}}(\hat{{\bf{s}}}){Y}_{{{lm}}^{\prime\prime }}(\hat{{\bf{s}}})\;{\rm{d}}\hat{{\bf{s}}}\\ & = & \displaystyle \sum _{m\prime \prime =-l}^{l}{d}_{0{m}^{\prime\prime }}^{l}({\theta }_{\hat{{\bf{k}}}})\sqrt{(2l+1)\pi }{(-1)}^{{m}^{\prime }}{\delta }_{m\prime m\prime \prime }{g}_{l}^{{m}^{\prime }}(-\xi )\\ & = & {\chi }_{[0,l]}(| m\prime | ){d}_{0{m}^{\prime }}^{l}({\theta }_{\hat{{\bf{k}}}})\sqrt{(2l+1)\pi }{(-1)}^{l}{g}_{l}^{{m}^{\prime }}(\xi ).\end{array}\end{eqnarray*}$$

Therefore,

$$\begin{eqnarray*}\begin{array}{rcl}{K}^{{m}^{\prime }}(\xi ,{{\bf{q}}}_{0}) & = & {\stackrel{\check{} }{K}}^{{m}^{\prime }}(\xi ,{\bf{q}})\delta ({\bf{q}}-{{\bf{q}}}_{0}),\\ {\stackrel{\check{} }{K}}^{{m}^{\prime }}(\xi ,{{\bf{q}}}_{0}) & = & \displaystyle \frac{-2{\pi }^{2}\varpi \xi }{\xi +{\hat{k}}_{z}(\xi {q}_{0})}\displaystyle \sum _{l=| m\prime | }^{L}{(-1)}^{l}{\beta }_{l}{d}_{0{m}^{\prime }}^{l}({\theta }_{\hat{{\bf{k}}}}){g}_{l}^{{m}^{\prime }}(\xi ).\end{array}\end{eqnarray*}$$

This implies that ${c}_{{lm}}({\bf{q}})$ have the form

$$\begin{eqnarray*}&&{c}_{{lm}}({\bf{q}})={\stackrel{\check{} }{c}}_{{lm}}({{\bf{q}}}_{0})\delta ({\bf{q}}-{{\bf{q}}}_{0}).\end{eqnarray*}$$

Since ${\stackrel{\check{} }{K}}^{{m}^{\prime }}(\xi ,{{\bf{q}}}_{0})$ is independent of ${\varphi }_{{{\bf{q}}}_{0}}$ and ${\stackrel{\check{} }{K}}^{-m\prime }={\stackrel{\check{} }{K}}^{{m}^{\prime }}$, we can write the key F_N equation as

$$\begin{eqnarray}&&\begin{array}{l}\displaystyle \sum _{m=0}^{{l}_{\mathrm{max}}}\displaystyle \sum _{\alpha =0}^{\left\lfloor ({l}_{\mathrm{max}}-m)/2\right\rfloor}\left[{A}_{m+2\alpha ,m}^{{m}^{\prime }}(\xi ,{q}_{0})+(1-{\delta }_{m0}){(-1)}^{m}{A}_{m+2\alpha ,-m}^{{m}^{\prime }}(\xi ,{q}_{0})\right]{\stackrel{\check{} }{c}}_{m+2\alpha ,m}({q}_{0})\\ \qquad ={\stackrel{\check{} }{K}}^{{m}^{\prime }}(\xi ,{{\bf{q}}}_{0}).\end{array}\end{eqnarray} \tag{ 41 }$$

The number of columns of the matrix ${\{A({q}_{0})\}}_{{\xi }^{{m}^{\prime }},{lm}}={A}_{{lm}}^{{m}^{\prime }}(\xi ,{q}_{0})$ is ${N}_{\mathrm{tot}}$, where

$$\begin{eqnarray*}{N}_{\mathrm{tot}}=\displaystyle \sum _{m=0}^{{l}_{\mathrm{max}}}{N}_{\mathrm{col}}^{m}=\left\{\begin{array}{ll}\displaystyle \frac{{({l}_{\mathrm{max}}+2)}^{2}}{4}, & {l}_{\mathrm{max}}\mathrm{even},\\ \displaystyle \frac{({l}_{\mathrm{max}}+1)({l}_{\mathrm{max}}+3)}{4}, & {l}_{\mathrm{max}}\mathrm{odd},\end{array}\right.\end{eqnarray*}$$

where ${N}_{\mathrm{col}}^{m}=\left\lfloor ({l}_{\mathrm{max}}-m)/2\right\rfloor+1$. We choose the number of rows so that $A({q}_{0})$ becomes square. For this purpose, different collocation schemes have been proposed [14, 16, 20, 46]. Here we take, in addition to discrete eigenvalues ${\xi }_{j}={\nu }_{j-1}^{{m}^{\prime }}$ ($j=1,...,{M}^{{m}^{\prime }}$), ${N}_{\mathrm{col}}^{{m}^{\prime }}-{M}^{{m}^{\prime }}$ points according to

$$\begin{eqnarray}&&{\xi }_{j}=\mathrm{cos}\left(\displaystyle \frac{\pi }{2}\displaystyle \frac{j-{M}^{{m}^{\prime }}}{{N}_{\mathrm{col}}^{{m}^{\prime }}-{M}^{{m}^{\prime }}+1}\right),\quad j={M}^{{m}^{\prime }}+1,...,{N}_{\mathrm{col}}^{{m}^{\prime }}.\end{eqnarray} \tag{ 42 }$$

The number of components of the vector $\{\stackrel{\check{} }{{\bf{K}}}({q}_{0})\}{}_{{\xi }^{{m}^{\prime }}}={\stackrel{\check{} }{K}}^{{m}^{\prime }}(\xi ,{{\bf{q}}}_{0})$ is ${N}_{\mathrm{tot}}$.

The hemispheric flux ${J}_{+({\boldsymbol{\rho }};{{\bf{q}}}_{0})}$ exiting the boundary is

$$\begin{eqnarray*}\begin{array}{rcl}{J}_{+}({\boldsymbol{\rho }};{{\bf{q}}}_{0}) & = & {\displaystyle \int }_{0}^{2\pi }{\displaystyle \int }_{0}^{1}\mu I({\boldsymbol{\rho }},0,-\hat{{\bf{s}}})\;{\rm{d}}\mu {\rm{d}}\varphi \\ & \approx & \displaystyle \frac{1}{4{\pi }^{3/2}}{{\rm{e}}}^{-{\rm{i}}{{\bf{q}}}_{0}\cdot {\boldsymbol{\rho }}}\displaystyle \sum _{l=0,2,...}\sqrt{2l+1}{\stackrel{\check{} }{c}}_{l0}({{\bf{q}}}_{0}){\displaystyle \int }_{0}^{1}\mu {P}_{l}(\mu )\;{\rm{d}}\mu .\end{array}\end{eqnarray*}$$

Here for even l

$$\begin{eqnarray*}&&{\displaystyle \int }_{0}^{1}\mu {P}_{l}(\mu )\;{\rm{d}}\mu =\displaystyle \frac{-{(-1)}^{l/2}(l-1)!!}{l!!(l-1)(l+2)}=\displaystyle \frac{-{(-1)}^{l/2}l!}{{2}^{l}(l-1)(l+2){\left[\left(\displaystyle \frac{l}{2}\right)!\right]}^{2}}.\end{eqnarray*}$$

Therefore we obtain

$$\begin{eqnarray}&&{J}_{+}({\boldsymbol{\rho }};{{\bf{q}}}_{0})\approx \displaystyle \frac{1}{4{\pi }^{3/2}}{{\rm{e}}}^{-{\rm{i}}{{\bf{q}}}_{0}\cdot {\boldsymbol{\rho }}}\displaystyle \sum _{l=0,2,...}\displaystyle \frac{\sqrt{2l+1}{(-1)}^{1+l/2}l!}{{2}^{l}(l-1)(l+2){\left[\left(\displaystyle \frac{l}{2}\right)!\right]}^{2}}{\stackrel{\check{} }{c}}_{l0}({q}_{0}).\end{eqnarray} \tag{ 43 }$$

Let us express the absolute value as

$$\begin{eqnarray}&&{J}_{+}({q}_{0})=\left|{J}_{+}({\boldsymbol{\rho }};{{\bf{q}}}_{0})\right|.\end{eqnarray} \tag{ 44 }$$

The algorithm of the three-dimensional F_N method can be summarized as follows.

Step 1. The integral over μ in (36) is done using the Golub–Welsch algorithm [21] of the Gauss–Legendre quadrature with points ${\mu }_{n}$ and weights w_n ($n=1,2,...,{N}_{\mu }$). The integral over $\varphi$ in (36) is computed using the trapezoid rule with points ${\varphi }_{j}=2\pi j/{N}_{\varphi }$ ($j=0,1,...,{N}_{\varphi }$). We use eigenvalues of the matrix $B(m\prime )$ in (18) for ${\xi }_{j}^{{m}^{\prime }}$ corresponding to discrete eigenvalues and use (42) for ${\xi }_{j}^{{m}^{\prime }}$ corresponding to the continuous spectrum. We calculate ${P}_{l}^{m}({\mu }_{n})$ and ${g}_{l}^{m}({\xi }_{j}^{{m}^{\prime }})$ with recurrence relations. The polynomials ${g}_{l}^{m}(\xi )$ are evaluated according to [17, 18]. That is, when ξ is a discrete eigenvalue, we obtain ${g}_{l}^{m}(\xi )$ starting with a large degree using backward recursion. For ξ in the continuous spectrum, we begin with the initial term and successively obtain ${g}_{l}^{m}(\xi )$ using the three-term recurrence relation (6).

Step 2. The analytically continued Wigner d-matrices are computed using the recurrence relation. See appendix B.

Step 3. We compute the double integrals in (36). In the function ${g}^{{m}^{\prime }}$, we compute ${p}_{l}^{m}(\mu )$ by using the recurrence relation (9). The computation time for each double integral grows as ${N}_{\mu }{N}_{\varphi }$.

Step 4. The coefficients ${\stackrel{\check{} }{c}}_{{lm}}({q}_{0})$ are obtained from the linear system (41) with the ${N}_{\mathrm{tot}}\times {N}_{\mathrm{tot}}$ matrix $A({q}_{0})$ and the vector $\stackrel{\check{} }{{\bf{K}}}({q}_{0})$ of length ${N}_{\mathrm{tot}}$.

Step 5. Once ${\stackrel{\check{} }{c}}_{{lm}}({q}_{0})$ are obtained, ${J}_{+({\boldsymbol{\rho }};{{\bf{q}}}_{0})}$ is immediately calculated by using (43).

Remark 4.1. The computation time is dominated by the integral in (36), which does not exist in the method of rotated reference frames (appendix A). For a given ${{\bf{q}}}_{0}$, the computation time for the double integrals grows as $O({l}_{\mathrm{max}}^{5}{N}_{\mu }{N}_{\phi })$ whereas the computation time of ${J}_{+}({q}_{0})$ scales as $O({l}_{\mathrm{max}}^{5})$ in the method of rotated reference frames.

For numerical calculation, let us set the absorption and scattering coefficients to

$$\begin{eqnarray*}&&{\mu }_{a}=0.05,\qquad {\mu }_{s}=100.\end{eqnarray*}$$

We set the scattering asymmetry parameter to ${\rm{g}}=0.9$ and ${\rm{g}}=0.01$ (almost isotropic). Although the unit of length has been $1/{\mu }_{t}$, we take the transport mean free path ${{\ell }}^{*}=1/({\mu }_{t}-{\mu }_{s}{\rm{g}})$ to be the unit of length in the figures.

In figures 2 and 3, ${J}_{+}({q}_{0})$ in (44) is plotted as a function of the spatial frequency q₀. The F_N result is compared with Monte Carlo simulation and the method of rotated reference frames. In Monte Carlo simulation 10⁸ particles were used. To obtain Monte Carlo simulation for structured illumination, Fourier transform was performed to results from Monte Carlo simulation for the delta-function source [35]. Monte Carlo simulation assumed the Henyey–Greenstein model for the scattering phase function. The method of rotated reference frames for structured illumination [33, 35] is summarized in appendix A.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The scattering asymmetry parameter ${\rm{g}}=0.01$ in figure 2 and ${\rm{g}}=0.9$ in figure 3. We set $L={l}_{\mathrm{max}}$. For both the F_N method and the method of rotated reference frames we consider ${l}_{\mathrm{max}}=9$ and 25. In figure 2, the three methods agree reasonably well for ${l}_{\mathrm{max}}=9$. When we increase ${l}_{\mathrm{max}}$ aiming at more accuracy, however, ${J}_{+}$ from the method of rotated reference frames becomes unstable. Note that in this case scattering is almost isotropic and discrete eigenvalues are rather close to 1. Hence we have $\nu -\mu {\hat{k}}_{z}(\nu {q}_{0})\lt 0$ (see (30)) for relatively small q₀. In figure 3, the result from the 3D F_N method has a jump near ${q}_{0}=3.7$ for ${l}_{\mathrm{max}}=9$ because this ${l}_{\mathrm{max}}$ is not sufficiently large in this case. A smooth curve is obtained if large enough ${l}_{\mathrm{max}}$ is used as shown in the right panel of figure 3 for ${l}_{\mathrm{max}}=25$.

The F_N method is similar to the method of rotated reference frames in the sense that spherical-harmonic expansion is used. However, in the F_N method, there is no need of expanding singular eigenfunctions. The extension of the F_N method in the half space to the slab geometry is straight forward. In the slab geometry, in addition to conditions such as (4) for one plane at z = 0, we have another set of conditions that corresponds to the other plane. Once the specific intensity on the boundary is obtained, it is also possible to compute the specific intensity inside the medium for the half space geometry and the slab geometry [51].

The author learned the F_N method at the 23rd International Conference on Transport Theory (September 2013, Santa Fe, New Mexico). He is grateful to C E Siewert and J C Schotland for fruitful discussion on collision expansion. Monte Carlo simulation was performed using the Monte Carlo solver MC developed by V A Markel.

In this section we solve (1) with the method of rotated reference frames [33, 35]. We consider structured illumination and assume the source term (39) with (40).

We write the eigenvector of the matrix $B(M)$ in (18) corresponding to the eigenvalue ν as $| {y}_{\nu }\rangle$ ($\langle {y}_{\nu }| {y}_{\nu }\rangle =1$). Note that ν and $| {y}_{\nu }\rangle$ depend on M. In the method of rotated reference frames, we write the specific intensity as a superposition of ${I}^{(+)}$ and ${I}^{(-)}$ [41, 48], where

$$\begin{eqnarray*}\begin{array}{rcl}{I}_{M}^{(+)}({\bf{r}},\hat{{\bf{s}}}) & = & {{\rm{e}}}^{{\rm{i}}{\bf{q}}\cdot {\boldsymbol{\rho }}-{\hat{k}}_{z}(\nu q)z/\nu }\displaystyle \sum _{l=0}^{{l}_{\mathrm{max}}}\sqrt{\displaystyle \frac{2l+1}{{h}_{l}}}\displaystyle \sum _{m=-l}^{l}{Y}_{{lm}}(\hat{{\bf{s}}}){(-1)}^{m}{{\rm{e}}}^{-{\rm{i}}m{\varphi }_{{\bf{q}}}}\langle l| {y}_{\nu }\rangle {d}_{{mM}}^{l}[{\rm{i}}\tau (\nu q)],\\ {I}_{M}^{(-)}({\bf{r}},\hat{{\bf{s}}}) & = & {{\rm{e}}}^{{\rm{i}}{\bf{q}}\cdot {\boldsymbol{\rho }}+{\hat{k}}_{z}(\nu q)z/\nu }\displaystyle \sum _{l=0}^{{l}_{\mathrm{max}}}\sqrt{\displaystyle \frac{2l+1}{{h}_{l}}}\displaystyle \sum _{m=-l}^{l}{Y}_{{lm}}(-\hat{{\bf{s}}}){{\rm{e}}}^{-{\rm{i}}m{\varphi }_{{\bf{q}}}}\langle l| {y}_{\nu }\rangle {d}_{m,-M}^{l}[{\rm{i}}\tau (\nu q)],\end{array}\end{eqnarray*}$$

In the half space ${{\mathbb{R}}}_{+}^{3}$, the specific intensity is given by

$$\begin{eqnarray*}&&I({\bf{r}},\hat{{\bf{s}}})\approx \displaystyle \frac{1}{{(2\pi )}^{2}}\displaystyle \sum _{M=-L}^{L}\displaystyle \sum _{\nu }{\displaystyle \int }_{{{\mathbb{R}}}^{2}}{F}_{M}^{(+)}{I}_{M}^{(+)}({\bf{r}},\hat{{\bf{s}}})\;{\rm{d}}{\bf{q}},\end{eqnarray*}$$

where ${\displaystyle \sum }_{\nu }$ stands for the sum over all positive eigenvalues of $B(M)$. From the boundary conditions we obtain

$$\begin{eqnarray*}&&{F}_{M}^{(+)}={f}_{M}^{(+)}(q){(2\pi )}^{2}\delta ({q}_{x}+{q}_{0})\delta ({q}_{y}),\qquad {f}_{-M}^{(+)}(q)={(-1)}^{M}{f}_{M}^{(+)}(q).\end{eqnarray*}$$

Here ${f}_{M}^{(+)}(q)$ are solutions to

$$\begin{eqnarray*}&&{\mathcal{M}}(q){f}_{M}^{(+)}(q)={v}^{(+)},\quad M\geqslant 0,\end{eqnarray*}$$

where

$$\begin{eqnarray*}&&{\{{\mathcal{M}}(q)\}}_{{lm},\nu }=\displaystyle \sum _{l\prime =0}^{{l}_{\mathrm{max}}}\sqrt{\displaystyle \frac{2l\prime +1}{{h}_{{l}^{\prime }}}}{{\mathcal{B}}}_{{{ll}}^{\prime }}^{m}\langle l\prime | {y}_{\nu }\rangle \left({d}_{{mM}}^{{l}^{\prime }}[i\tau (\nu q)]+(1-{\delta }_{M0}){(-1)}^{M}{d}_{m,-M}^{{l}^{\prime }}[i\tau (\nu q)]\right),\end{eqnarray*}$$

and

$$\begin{eqnarray*}&&{\{{v}^{(+)}\}}_{{lm}}={\delta }_{m0}\displaystyle \sum _{l\prime =0}^{{l}_{\mathrm{max}}}{{\mathcal{B}}}_{{{ll}}^{\prime }}^{0}\sqrt{\displaystyle \frac{2l\prime +1}{4\pi }}.\end{eqnarray*}$$

Here,

$$\begin{eqnarray*}&&{{\mathcal{B}}}_{{{ll}}^{\prime }}^{m}=\displaystyle \frac{1}{2}\sqrt{\displaystyle \frac{(2l+1)(2l\prime +1)(l-m)!(l\prime -m)!}{(l+m)!(l\prime +m)!}}{\displaystyle \int }_{0}^{1}{P}_{l}^{m}(\mu ){P}_{{l}^{\prime }}^{m}(\mu )\;{\rm{d}}\mu .\end{eqnarray*}$$

That is,

$$\begin{eqnarray*}&&I({\bf{r}},\hat{{\bf{s}}})\approx \displaystyle \frac{1}{2\pi }\displaystyle \sum _{l=0}^{{l}_{\mathrm{max}}}\displaystyle \sum _{m=0}^{l}{i}^{m}\sqrt{\displaystyle \frac{2l+1}{{h}_{l}}}\left[{Y}_{{lm}}(\hat{{\bf{s}}})+(1-{\delta }_{m0}){Y}_{{lm}}^{*}(\hat{{\bf{s}}})\right]{K}_{{lm}}({\boldsymbol{\rho }},z),\end{eqnarray*}$$

where

$$\begin{eqnarray*}\begin{array}{rcl}{K}_{{lm}}({\boldsymbol{\rho }},z) & = & 2\pi {(-{\rm{i}})}^{m}{{\rm{e}}}^{{\rm{i}}{\bf{q}}\cdot {\boldsymbol{\rho }}}\displaystyle \sum _{M\geqslant 0}\displaystyle \sum _{\nu }\langle l| {y}_{\nu }\rangle {{\rm{e}}}^{-{\hat{k}}_{z}(\nu {q}_{0})z/\nu }{f}_{M}^{(+)}({q}_{0})\\ & & \times \left[{d}_{{mM}}^{l}[{\rm{i}}\tau ({q}_{0}\nu )]+(1-{\delta }_{M0}){(-1)}^{M}{d}_{m,-M}^{l}[{\rm{i}}\tau ({q}_{0}\nu )]\right].\end{array}\end{eqnarray*}$$

The hemispheric flux is obtained as

$$\begin{eqnarray}\begin{array}{rcl}{J}_{+}({\boldsymbol{\rho }}) & = & {\displaystyle \int }_{0}^{2\pi }{\displaystyle \int }_{-1}^{0}(\hat{{\bf{s}}}\cdot \hat{{\bf{z}}})I({\boldsymbol{\rho }},0,\hat{{\bf{s}}})\;{\rm{d}}\mu {\rm{d}}\varphi \\ & = & {\displaystyle \int }_{0}^{2\pi }{\displaystyle \int }_{-1}^{1}(\hat{{\bf{s}}}\cdot \hat{{\bf{z}}})I({\boldsymbol{\rho }},0,\hat{{\bf{s}}})\;{\rm{d}}\mu {\rm{d}}\varphi -{\displaystyle \int }_{0}^{2\pi }{\displaystyle \int }_{0}^{1}(\hat{{\bf{s}}}\cdot \hat{{\bf{z}}})I({\boldsymbol{\rho }},0,\hat{{\bf{s}}})\;{\rm{d}}\mu {\rm{d}}\varphi \\ & = & \displaystyle \frac{1}{\sqrt{\pi {h}_{1}}}{K}_{10}({\boldsymbol{\rho }},0)-{{\rm{e}}}^{-{\rm{i}}{q}_{0}x}{\mu }_{0}{\chi }_{[\mathrm{0,1}]}({\mu }_{0}),\end{array}\end{eqnarray} \tag{ A.1 }$$

where we used $I({\boldsymbol{\rho }},0,\hat{{\bf{s}}})={{\rm{e}}}^{-{\rm{i}}{q}_{0}x}\delta (\hat{{\bf{s}}}-{\hat{{\bf{s}}}}_{0})$ for $\mu \gt 0$. The expression (A.1) is used for figures 2 and 3.

To compute the analytically continued Wigner d-matrices we use a pyramid scheme with recurrence relations [3]. We begin with ${d}_{00}^{0}[{\rm{i}}\tau (x)](=1)$, ${d}_{00}^{1}[{\rm{i}}\tau (x)]$, ${d}_{1-1}^{1}[{\rm{i}}\tau (x)]$, ${d}_{10}^{1}[{\rm{i}}\tau (x)]$, and ${d}_{11}^{1}[{\rm{i}}\tau (x)]$:

$$\begin{eqnarray*}&&{d}_{00}^{1}=\sqrt{1+{x}^{2}},\quad {d}_{1-1}^{1}=\displaystyle \frac{1-\sqrt{1+{x}^{2}}}{2},\quad {d}_{10}^{1}=-{\rm{i}}\displaystyle \frac{x}{\sqrt{2}},\quad {d}_{11}^{1}=\displaystyle \frac{1+\sqrt{1+{x}^{2}}}{2}.\end{eqnarray*}$$

Let us we increase l iteratively up to ${l}_{\mathrm{max}}$. For each value of l, we first compute ${d}_{{{mm}}^{\prime }}^{l}[{\rm{i}}\tau (x)]$ ($m=0,...,l-2;m\prime =-m,...,m$) according to

$$\begin{eqnarray*}\begin{array}{rcl}{d}_{{{mm}}^{\prime }}^{l} & = & \displaystyle \frac{l(2l-1)}{\sqrt{({l}^{2}-{m}^{2})({l}^{2}-m{\prime }^{2})}}\\ & & \times \left[\left({d}_{00}^{1}-\displaystyle \frac{{mm}\prime }{l(l-1)}\right){d}_{{{mm}}^{\prime }}^{l-1}-\displaystyle \frac{\sqrt{\left[{(l-1)}^{2}-{m}^{2}\right]\left[{(l-1)}^{2}-m{\prime }^{2}\right]}}{(l-1)(2l-1)}{d}_{{{mm}}^{\prime }}^{l-2}\right].\end{array}\end{eqnarray*}$$

We obtain ${d}_{{ll}}^{l}[{\rm{i}}\tau (x)]$ and ${d}_{l-1,l-1}^{l}[{\rm{i}}\tau (x)]$ as

$$\begin{eqnarray*}&&{d}_{{ll}}^{l}={d}_{11}^{1}{d}_{l-1,l-1}^{l-1},\qquad {d}_{l-1,l-1}^{l}=({{ld}}_{00}^{1}-l+1){d}_{l-1,l-1}^{l-1},\end{eqnarray*}$$

and ${d}_{{{lm}}^{\prime }}^{l}[\tau (x)]$ $(m\prime =l-1,...,-l)$ as

$$\begin{eqnarray*}&&{d}_{{{lm}}^{\prime }}^{l}=-{\rm{i}}\sqrt{\displaystyle \frac{l+m\prime +1}{l-m\prime }}\sqrt{\left|\displaystyle \frac{{d}_{1-1}^{1}}{{d}_{11}^{1}}\right|}{d}_{l,m\prime +1}^{l}.\end{eqnarray*}$$

With the relation

$$\begin{eqnarray*}&&{d}_{l-1,m\prime }^{l}=-{\rm{i}}\displaystyle \frac{{{ld}}_{00}^{1}-m\prime }{{{ld}}_{00}^{1}-m\prime -1}\sqrt{\displaystyle \frac{l+m\prime +1}{l-m\prime }}\sqrt{\left|\displaystyle \frac{{d}_{1-1}^{1}}{{d}_{11}^{1}}\right|}{d}_{l-1,m\prime +1}^{l},\end{eqnarray*}$$

we have ${d}_{l-1,m\prime }^{l}[{\rm{i}}\tau (x)]$ ($m\prime =l-2,...,1-l$). Other functions ${d}_{{{mm}}^{\prime }}^{l}[{\rm{i}}\tau (x)]$ are obtained by using the symmetry properties

$$\begin{eqnarray*}&&{d}_{{{mm}}^{\prime }}^{l}={d}_{-m\prime ,-m}^{l}={(-1)}^{m+m\prime }{d}_{-m,-m\prime }^{l}={(-1)}^{m+m\prime }{d}_{m\prime m}^{l}.\end{eqnarray*}$$