Journal of Computational and Theoretical Transport

ISSN: 2332-4309 (Print) 2332-4325 (Online) Journal homepage: http://www.tandfonline.com/loi/ltty21

How to Construct Three-Dimensional Transport
Theory Using Rotated Reference Frames
Manabu MachidaTo cite this article: Manabu Machida (2016) How to Construct Three-Dimensional Transport
Theory Using Rotated Reference Frames, Journal of Computational and Theoretical Transport,
45:7, 594-609, DOI: 10.1080/23324309.2016.1239639
To link to this article: http://dx.doi.org/10.1080/23324309.2016.1239639

Published online: 01 Nov 2016.

Submit your article to this journal

View related articles

View Crossmark data

Full Terms & Conditions of access and use can be found at
http://www.tandfonline.com/action/journalInformation?journalCode=ltty21
Download by: [Manabu Machida]

Date: 01 November 2016, At: 16:08

JOURNAL OF COMPUTATIONAL AND THEORETICAL TRANSPORT
, VOL. , NO. , –
http://dx.doi.org/./..

How to Construct Three-Dimensional Transport Theory
Using Rotated Reference Frames
Manabu Machida
Institute for Medical Photonics Research, Hamamatsu University School of Medicine, Hamamatsu,
Shizuoka, Japan

ABSTRACTKEYWORDS

In linear transport theory, 3D equations reduce to 1D equations bymeans of rotated reference frames. In this paper, we illustrate how
the technique works and 3D transport theories are obtained.

Linear transport theory;
rotated reference frames;
singular eigenfunctions

1. Introduction
In linear transport theory, 3D equations reduce to 1D equations with rotated reference frames when the angular flux has the structure of separation of variables.
Although the technique can be used for anisotropic scattering in the presence ofboundaries, calculation becomes complicated in such general cases and sometimes
the essence is buried in straightforward but lengthy and tedious calculations. In
this paper, aiming at elucidating how rotated reference frames help construct linear
transport theory in three dimensions, we will consider the case of isotropic scattering in an infinite medium.
To the best of the author’s knowledge, rotated reference frames were first introduced in transport theory by Dede (1964) in the context of the PN method, andthen Kobayashi (1977) expanded Dede’s calculation. Forty years after Dede’s finding, Markel (2004) devised an efficient numerical algorithm of computing solutions
to the 3D transport equation by reinventing rotated reference frames. The method
is called the method of rotated reference frames. A lot of numerical calculations(Panasyuk et al., 2006; Xu and Patterson, 2006a,b; Machida et al., 2010; Liemert
and Kienle, 2011a,b,c, 2012a,b,c,d,e, 2013a,b,c, 2015) done during the last decade
have proved the usefulness and efficiency of this method. It was then found thatthe technique of rotated reference frames is not merely for the particular numericalmethod but is a tool to build bridges between 3D transport theory and 1D transport theory. Case’s method (Case, 1960; Mika, 1961; McCormick and Kuščer, 1966;
Case and Zweifel, 1967) was extended to three dimensions (Machida, 2014), andthe FN method (Siewert, 1978; Siewert and Benoist, 1979) was also extended to three
CONTACT Manabu Machidamachida@hama-med.ac.jp
Institute for Medical Photonics Research, Hamamatsu
University School of Medicine, Hamamatsu, Shizuoka, - Japan.
©  Taylor & Francis Group, LLC

JOURNAL OF COMPUTATIONAL AND THEORETICAL TRANSPORT

595

dimensions (Machida, 2015). Recently, the angular flux of the 3D transport equation
with anisotropic scattering was computed using the Fourier transform by making
use of rotated reference frames (Machida, 2016). The technique of rotated reference
frames is also applied to optical tomography (Schotland and Markel, 2007; Machida
et al., 2016; Machida, in press).
There have been many attempts to construct multi-dimensional transport theory
especially since Case’s method (1960) of singular eigenfunctions appeared. In addition to trials of direct extension of the singular-eigenfunction approach (Williams,
1967; Kaper, 1969; Case and Hazeltine, 1970; Gibbs, 1969; Cannon, 1973; Pomraning, 1996), there were approaches based on the integral transform (Bareiss and
Abu-Shumays, 1967; Williams, 1968; Schreiner et al., 1969; Garrettson and Leonard,
1970; Leonard, 1971; Lam and Leonard, 1971, 1973; Williams, 1982). However,
the efforts of extending the singular-eigenfunction approach ended up with formal
complicated calculations and could not arrive at useful formulas.
Let us consider the following transport equation in a 3D infinite medium:ˆ · ∇ + 1 ψ (r, )ˆ =ˆ dˆ + S(r, ),ˆψ (r, )(1)
4π S2
ˆ ∈ S2 , and  ∈ (0, 1) is the albedo for single scattering. Here,
where r ∈ R3 , ˆ is the angular flux and S(r, )ˆ is the source term. The unit vector ˆ has theψ (r, )polar angle θ and azimuthal angle ϕ. We let μ denote the cosine of θ , i.e., μ = cos θ .
ˆ we introduce an operator R
Let k̂ ∈ C3 be a unit vector. For a function f (),k̂ˆ is the value of f in which ˆ is measured in the reference frame
such that Rk̂ f ()ˆ = μ, we have
whose z-axis lies in the direction of k̂. For example, when f ()ˆRk̂ μ = k̂ · .
Such reference frames which are rotated in directions of given unit vectors k̂ arecalled rotated reference frames. We define spherical harmonics byˆ = 2l + 1 (l − m)! Pm (μ)eimϕ ,Ylm ()4π (l + m)! l
where Plm (μ) are associated Legendre polynomials. We have the following calculation for spherical harmonics:ˆ =Rk̂ Ylm ()lm =−lˆDlm m ϕk̂ , θk̂ , 0 Ylm (),where θk̂ and ϕk̂ are the polar and azimuthal angles of k̂ in the laboratory
frame, and Dlm m are Wigner’s D-matrices (Varshalovich et al., 1988). The function
Dlm m (α, β, γ ), where α, β, γ are called the Euler angles, expresses rotation of thereference frame. If we set α = ϕk̂ , β = θk̂ , and γ = 0, we can rotate the laboratory frame to the reference frame whose z-axis lies in the direction of k̂. That is,
by Dlm m (ϕk̂ , θk̂ , 0), the reference frame is rotated about the original z-axis by ϕk̂

and then about the new x-axis by θk̂ . As a result, the new z-axis coincides with k̂. If

596

M. MACHIDA

γ is non-zero, we further rotate the reference frame about the resulting z-axis by theangle γ . We note that if we rotate the obtained rotated reference frame about the xaxis by −θk̂ and then rotate it on the x–y plane by −ϕk̂ , the reference frame returns
to the laboratory frame. The rotated reference frame also comes back to the laboratory frame by the following successive rotations: the rotation about the new z-axis
or the vector k̂ by π , the rotation about the flipped x-axis by θk̂ , and the rotation on
the x–y plane by π − ϕk̂ . From the above consideration, we haveˆ =Ylm ()lm =−l=lm =−lˆDlm m 0, −θk̂ , −ϕk̂ Rk̂ Ylm ()ˆDlm m π , θk̂ , π − ϕk̂ Rk̂ Ylm ().


l
−imγ
l
, where dm
It is possible to write Dlm m (α, β, γ ) = e−im α dm
 m (β )e m (β ) arecalled Wigner’s d-matrices. We haves!(s + m1 + m2 )! 1 − cos β m1 /2 1 + cos β m2 /2ldm m (β ) = ξm m(s + m1 )!(s + m2 )!
2
2

×Ps(m1 ,m2 ) (cos β ),where ξm m = 1 for m ≥ m , (−1)m−m for m < m , m1 = |m − m|, m2 = |m + m|,
s = l − (m1 + m2 )/2, and the Jacobi polynomials Ps(m1 ,m2 ) (μ) are given byPs(m1 ,m2 ) (μ) =(s + m1 )!
(s + m2 )!
1 (μ − 1) j (μ + 1)s− j .
s
2 j=0 (s − j)!(m1 + j)! j!(s + m2 − j)!
s





l
m −m ll
We note that dm
d−m ,−m (β ) = (−1)m −m dmm
 m (β ) = (−1)
 (β ) and0(β ) = 1,
d00

−1
1
d10
(β ) = √ sin β,
21d00(β ) = cos β,1(β ) =
d1±1

1 ± cos β.
2

2. Case’s method in three dimensions
For a unit vector k̂, we assume the following separated solutions with separation
parameter ν:ˆ k̂) =ψν (r, ;
where the unknown function12π−k̂·r/νˆ,ν (; k̂)eˆ

ν (; k̂) is normalized as



ˆ

S2

ˆ

ν (; k̂) d  = 1.

ˆ k̂) into the homogeneous equation
We plug the above ψν (r, ;ˆ d ,ˆˆˆ
ψ (r, ) · ∇ + 1 ψ (r, ) =
4π S2

(2)

JOURNAL OF COMPUTATIONAL AND THEORETICAL TRANSPORT

597

and obtain
1−
Let us expressˆk̂ · νˆν (; k̂) =


.
2

(3)

ˆ

ν (; k̂) asˆν (; k̂) = Rk̂ φ(ν, μ),

where

 1
−1

φ(ν, μ) dμ = 1.

We will see below that φ(ν, μ) is independent of ϕ. Then, Equation (3) is writtenasR μ
1 − k̂
Rk̂ φ(ν, μ) = .
(4)
ν2
Since the right-hand side of Equation (4) is a scalar, by operating R−1 , Equation (4)k̂reduces to

μ

1−
(5)
φ(ν, μ) = .
ν
2The above Equation (5) is the equation appearing in 1D transport theory (Case,
1960; Case and Zweifel, 1967). This is different from the pseudo-problem approach(Williams, 1967, 1968, 1982) and the coefficients in the corresponding 1D transport
equation remains constant. That is, in one dimension, Equation (5) is derived fromthe following equation:∂ 1
μ + 1 ψ (z, μ) =
ψ (z, μ) dμ,
∂z
2 −1
where
ψ (z, μ) = φ(ν, μ)e−z/ν .
Therefore, it turns out that φ(ν, μ) are singular eigenfunctions (Case, 1960), which
are obtained as1νP+ λ(ν )δ(ν − μ),φ(ν, μ) =
2
ν−μ
where

 1
ν
1
λ(ν ) = 1 −P
dμ = 1 −  ν tanh−1 (ν ).
2−1 ν − μThe separation constant ν takes values in (−1, 1) in addition to ±ν0 , where ν0 > 1
is the positive root of (ν ) such that (ν0 ) = 0. Here, the function (w) is definedfor w ∈ C \ [−1, 1] as

1
w 1
dμ.
(6)
(w) = 1 −2 −1 w − μ598

M. MACHIDA

Thus, ν0 is given as the positive solution to the transcendental equation
1 −  ν0 tanh−1

1
ν0

= 0.

When  is near 1, which is typical for light propagating in biological tissue, ν0 is
approximately calculated as Case and Zweifel (1967)
ν0 ≈ √

1.
3(1 −  )Now, we return to eigenmodes (2) in three dimensions. We obtainˆ k̂) = R φ(ν, μ)e−k̂·r/ν = φ(ν, k̂ · )eˆ −k̂·r/ν .
ψν (r, ;k̂

(7)

The angular flux in Equation (1) is given by a superposition of eigenmodesˆ k̂).
ψν (r, ;So far, k̂ ∈ C3 has been arbitrary as long as k̂ · k̂ = 1 is satisfied. Let us write k̂ =(k̂x , k̂y , k̂z )T . We see that k̂ · k̂ = 1 cannot be achieved if k̂x = k̂y = k̂z = 0.
Therefore, at least one component of k̂ must have a non-zero real part. Hereafter,
we consider k̂ of the following form, so that the Green’s function (10) is given as
plane-wave decomposition:k̂(ν, q) =

−iνq

,

k̂z (νq)

where q ∈ R2 , q = |q|, and k̂z (νq) = 1+(νq)2 . Below, the Green’s function (10)is given as a superposition of all possible k̂ and ν. It turns out that the above form ofk̂(ν, q) is enough to construct a complete set (Panasyuk et al., 2006; Kim, 2004). The
present formulation is similar to the pseudo-problem approach (Williams, 1967,
1968, 1982) in the sense that Fourier transform is taken for the 2D position vector ρ. However, as is seen below, the structure of the Green’s function (10) becomesmuch simpler thanks to the use of rotated reference frames. The angles θk̂ and ϕk̂
are introduced as⎞⎛sin θk̂ cos ϕk̂⎟⎜k̂(ν, q) = ⎝ sin θk̂ sin ϕk̂ ⎠ .
cos θk̂
1 − cos2 θk̂ = i|νq| by√z for any
choosing the branch cut on the positive real axis, i.e., we have 0 ≤z ∈ C. We can write the vector q as

We see that cos θk̂ = k̂z (νq). Moreover, we have sin θk̂ =q=q

cos ϕqsin ϕq.

JOURNAL OF COMPUTATIONAL AND THEORETICAL TRANSPORT599

We have ϕk̂ = ϕq + π (ν > 0) and ϕk̂ = ϕq (ν < 0). We can rotate the reference
frame using k̂. For this k̂, we have
4π
ˆR
Y10 ()Rk̂(ν,q) μ =3 k̂(ν,q)
1
4π  1ˆ=D  (ϕ , θ , 0)Y1m ()3 m =−1 m 0 k̂ k̂1
(1 − m )! m1
P1 (μ)eim (ϕ−ϕk̂ )dm=
 0 (θ )k̂(1 + m )!
m =−1= cos θk̂ cos θ − 1 − cos2 θk̂ sin θ cos(ϕ − ϕk̂ )= k̂z (νq)μ − iνq 1 − μ2 cos(ϕ − ϕq ),1(θk̂ ) =where we used P1−1 (μ) = 12 sin θ , P1 (μ) = cos θ , P11 (μ) = − sin θ , d00
11√
cos θk̂ , and d±10 (θk̂ ) = ± 2 sin θk̂ . Similarly, we obtain

R−1 μ =k̂(ν,q)

=


4π −1
ˆ
R
Y10 ()3 k̂(ν,q)1
4π  1ˆD  (0, −θk̂ , −ϕk̂ )Y1m ()3 m =−1 m 0

1
4π  1ˆ=D  (π , θk̂ , π − ϕk̂ )Y1m ()3 m =−1 m 01(1 − m )! m1
P1 (μ)eim ϕ(−1)m dm= 0 (θ )k̂(1 + m )!
m =−1= cos θk̂ cos θ + 1 − cos2 θk̂ sin θ cos ϕ= k̂z (νq)μ − i|νq| 1 − μ2 cos ϕ.
Theorem 2.1 (Orthogonality relations (Machida, 2014, 2015)). For 3D singulareigenfunctions, we have
ˆ = 2π k̂z (νq)N (ν )δνν  .
μ Rk̂(ν,q) φ(ν, μ) Rk̂(ν  ,q) φ(ν  , μ) d 
S2Here, the Kronecker delta δνν  is read as the Dirac delta function δ(ν − ν  ) when ν, ν 
are in the continuous spectrum (−1, 1). The normalization factor N (ν ) is from 1Dtransport theory and given by
⎧


 3
⎨
ν ν 2−1 − ν12 ,ν = ±ν0 ,2N (ν ) =⎩ ν 1 −  ν tanh−1 (ν ) 2 +  πν 2 ,ν ∈ (−1, 1).
2

600

M. MACHIDA

Proof. For fixed q, we consider two eigenvalues ν1 and ν2 . Correspondingly, we writek̂1 = k̂(ν1 , q) and k̂2 = k̂(ν2 , q). We multiply Equation (3) for ν1 by Rk̂2 φ(ν2 , μ)and multiply Equation (3) for ν2 by Rk̂1 φ(ν1 , μ). We haveˆk̂1 · R φ(ν2 , μ),Rk̂1 φ(ν1 , μ) =Rk̂2 φ(ν2 , μ) 1 −
ν1
2 k̂2


ˆ
k̂2 · R φ(ν1 , μ).
Rk̂2 φ(ν2 , μ) =Rk̂1 φ(ν1 , μ) 1 −ν2
2 k̂1
By integrating both sides and subtracting the latter from former, we obtain
Rk̂2 μ Rk̂1 μ ˆ = 0.
Rk̂1 φ(ν1 , μ) Rk̂2 φ(ν2 , μ) d 
−
ν2
ν1
S2
Therefore,


S2




ˆ = 0,
μ Rk̂1 φ(ν1 , μ) Rk̂2 φ(ν2 , μ) d 

ν 1 = ν2 .

When ν1 = ν2 = ν, we have

 


2
ˆ
ˆ
μ Rk̂ φ(ν, μ) d  =R−1 μ φ(ν, μ)2 d 
S2

k̂

S2

= 2π k̂z (νq)
The proof is completed by noticing N (ν ) =1967).

1

 1−1

μφ(ν, μ)2 dμ.

2−1 μφ(ν, μ) dμ (Case and Zweifel,

By using 3D singular eigenfunctions, let us compute the Green’s function. Whenthe source term in Equation (1) is given byˆ = δ(r)δ(ˆ −ˆ 0 ),
S(r, )(8)

the angular flux becomes the Green’s function:ˆˆ ˆ 0 ) = ψ (r, ).
G(r, ;(9)ˆ ˆ 0 ) in terms of theWith some coefficients a± (q), Aν (q), we can write G(r, ;eigenmodes ψν in Equation (7) as
 


1
ˆ
ˆ
ˆ k̂(ν0 , q)aG(r, ; 0 ) =
(q)ψ
r,
;
+
ν
0
(2π )2 R2ˆ k̂(−ν0 , q)
+ a− (q)ψ−ν0 r, ;
+

 1
−1

Aν (q)ψν ˆ
r, ; k̂(ν, q) dν dq.

JOURNAL OF COMPUTATIONAL AND THEORETICAL TRANSPORT

601

ˆ 
ˆ 0 ) → 0 as
Since the Green’s function vanishes at infinity, particularly G(r, ;ˆ ˆ 0 ) as|z| → ∞, we can write G(r, ;ˆ ˆ 0)G(r, ;⎧ ⎪ˆ k̂) + 1 Aν (q)ψν (r, ;ˆ k̂) dν dq, z > 0,
⎨ (2π1 )2 R2 a+ (q)ψν0 (r, ;
0
=




⎪ −1
ˆ k̂) + 0 Aν (q)ψν (r, ;ˆ k̂) dν dq, z < 0,
⎩
(q)ψ
(r,
;
a
2
−
−ν
2
0
−1
(2π ) R
where coefficients a± (q), Aν (q) are determined below. We let ρ ∈ R2 be the position
vector in the plane perpendicular to the z-axis, i.e., r = (ρ, z) and ρ = (x, y). The
jump condition is written asˆ −ˆ 0 ).
ˆ ˆ 0 ) − G(ρ, 0− , ;ˆ ˆ 0 ) = 1 δ(ρ)δ(G(ρ, 0+ , ;μLet us multiply e−iq·ρ on both sides of the above-mentioned jump condition andintegrate both sides over ρ. As a result, we obtain
a+ (q) Rk̂(ν0 ,q) φ(ν0 , μ) + a− (q) Rk̂(−ν0 ,q) φ(−ν0 , μ)
 11 ˆˆ 0 ).
−Aν (q) Rk̂(ν,q) φ(ν, μ) dν = δ(+μ−1
By using the orthogonality relations in Theorem 2.1, we can determine the coefficients as
1
a± (q) =Rk̂(±ν0 ,q) φ(±ν0 , μ0 ),2π k̂z (ν0 q)N (±ν0 )
1Aν (q) =We note that

2π k̂z (νq)N (ν )Rk̂(ν,q) φ(ν, μ0 ).

ˆ .
Rk̂(ν,q) φ(ν, μ) = φ ν, k̂(ν, q) · 

Finally, the Green’s function is obtained as Machida (2014)ˆ ˆ 0) =G(r, ;×+1(2π )31k̂z (ν0 q)N (ν0 )
 1
1
0

R2

eiq·ρ

 


ˆ φ ±ν0 , k̂(±ν0 , q) · ˆ 0 e∓k̂z (ν0 q)z/ν0
φ ±ν0 , k̂(±ν0 , q) · k̂z (νq)N (ν)

 ∓k̂z (νq)z/νˆˆdν dq,
φ ±ν, k̂(±ν, q) ·  φ ±ν, k̂(±ν, q) · 0 e
(10)where upper signs are used for z > 0 and lower signs are chosen for z < 0. We canreadily extend Equation (10) to the case of anisotropic scattering (Machida, 2014).
The right-hand side of Equation (10) has the form of plane-wave decomposition.
The fact implies that Rk̂ φ are complete if φ form a complete set in the 1D transport

602

M. MACHIDA

theory (Panasyuk et al., 2006; Kim, 2004). The structure of the right-hand side ofEquation (10) is similar to Williams (1967) and Garrettson and Leonard (1970) in
the sense that Fourier transforms are used in the (x, y) directions and the dependence of the integrand on z is exponentially decaying.
1
Let us consider the energy density U (ρ, z) for an isotropic point source 4π
δ(r),
i.e.,
 
1
ˆ 
ˆ 0 ) d dˆ ˆ 0.
U (ρ, z) =G(r, ;4π S2 S2
It is worth mentioning that the singular eigenfunctions do not appear in U (ρ, z).
We have"!
 1 ∓k̂z (νq)z/νe∓k̂z (ν0 q)z/ν0e1
iq·ρU (ρ, z) =+dν dq.
e
8π 2 R20 k̂z (νq)N (ν )k̂z (ν0 q)N (ν0 )Let us suppose z > 0. We note that for ν > 0e−k̂z (νq)z/ν
2π
dq =
eiq·ρ
ν
R2
k̂z (νq) ∞J0 (qρ )

0

 ∞

e−z

√

(1/ν 2 )+q2

(1/ν 2 ) + q2

q dq

ρ 
√
u
2
u e−(z/ν) u +1 √du
ν
u2 + 1
0
√
2π
1
2
2
= 2
e− (ρ/ν) +(z/ν) ,22ν(ρ/ν ) + (z/ν )

2π
= 2
ν

J0

where J0 is the Bessel function of the zeroth order and we used u = νq. Thus, we
obtain −r/ν0 1 −r/ν1ee
U (ρ, z) =
+
dν ,
(11)
4π r ν0 N (ν0 )0 νN (ν )where r = ρ 2 + z2 . In some special cases such as the above U (ρ, z), it is possible
to arrive at the same expression without introducing the 3D singular eigenfunctions.
Indeed, in Ganapol and Kornreich (1995), Ganapol and Kornreich obtained exactly
the same formula (11) using the Fourier transform and pseudo-problem. In Ganapolˆ ẑ) d ˆ and the case of an isotropic line source are
and Kornreich (1995), S2 G(r, ;also calculated.
3. Ganapol’s Fourier-transform approach
Alternative expressions of the Green’s function (10) can be obtained with the Fouriertransform. We obtain (Ishimaru, 1978)ˆ ˆ 0 ) = 1 e−r δ ˆ − r δ(ˆ −
ˆ 0)G(r, ;
r2
r

−1

1 − k tan−1 (k)

ik·r
+
dk,
(12)
e
ˆ
ˆ 0)
2(2π )4 R3(1 + ik · )(1
+ ik · JOURNAL OF COMPUTATIONAL AND THEORETICAL TRANSPORT

603

where r = |r| and k = |k|. The extension of Equation (12) to anisotropic scattering is also possible by using rotated reference frames (Machida, 2016). Cassell andWilliams proposed an alternative expression of the Green’s function (12) by separating the once-collided term in addition to the uncollided term. This expression is more useful for benchmarking purposes (Cassell and Williams, 2000). In
one dimension, Ganapol has developed an alternative Fourier-transform approach(Ganapol, 2000, 2015), which is different from the conventional derivation thatyields the 1D version of Equation (12). The new formula is potentially more suitable for numerical calculation. Here, we compute the 3D Green’s function using
Ganapol’s approach.
In this section, we use the source term (8). The angular flux or the Green’s function is then symmetric about the azimuthal angle. Using the Fourier transform of
the Green’s function

ˆ
ˆ
ˆ 
ˆ 0 ) dr,G̃(k, ; 0 ) =e−ik·r G(r, ;R3

in the Fourier space, the transport equation is written as



ˆ 
ˆ 0) = ˆ G̃(k, ;ˆ ˆ 0 ) d
ˆ + δ(ˆ −ˆ 0 ).
1 + ik · G̃(k, ;4π S2Note thatˆ ˆ 0) =G(r, ;1(2π )3


R3

(13)

ˆ 
ˆ 0 )eik·r dk,
G̃(k, ;ˆ and rhas the structure of separation of variables in the sense that G̃ depends on 
ik·rexists only in e . Keeping in mind the notation (9), let us define



ˆ 
ˆ 0 ) d .
ˆRk̂ Pl (μ) G̃(k, ;(14)ψ̃l (k) =S2

We will use a new variableiz= .
k
Noting the recurrence relation(2l + 1)μPl (μ) = (l + 1)Pl+1 (μ) + lPl−1 (μ),
we obtain
zhl ψ̃l (k) − (l + 1)ψ̃l+1 (k) − l ψ̃l−1 (k) = zSl (k̂),
wherehl = 2l + 1 −  δl0 ,andˆ 0 ).
Sl (k̂) = (2l + 1) Rk̂ Pl (μ0 ) = (2l + 1)Pl (k̂ · 

(15)

604

M. MACHIDA

Chandrasekhar polynomials of the first and second kinds are defined aszhl gl (z) − (l + 1)gl+1 (z) − lgl−1 (z) = 0,g0 (z) = 1,

g1 (z) = z(1 −  ),andzhl ρl (z) − (l + 1)ρl+1 (z) − lρl−1 (z) = 0,

ρ0 (z) = 0,

ρ1 (z) = z.

We can express ψ̃l asψ̃l = a(z)gl (z) + b(z)ρl (z) + (1 − δl0 ) zl

αl, j (z)S j .

(16)

j=1

By setting l = 0 in Equation (16), we first notice thata(z) = ψ̃0 .
By plugging Equation (16) into Equation (15), we have− b(z)ρ1 (z) + zα1,1 (z)S1 = z.
Suppose l > 0. Let us impose
zhl αl, j − (l + 1)αl+1, j − lαl−1, j = 0.

(17)

By substituting Equation (16) for ψ̃l in Equation (15), we obtainzhl αl,l Sl − (l + 1) αl+1,l Sl + αl+1,l+1 Sl+1 = Sl .
The left-hand side of the above equation can be rewritten as
LHS = −(l + 1)αl+1,l+1 Sl+1 + lαl−1,l Sl .
Hence, we can put1
αl−1,l = ,l

αl+1,l+1 = 0.

(18)

Thus, we findb(z) = −1.
To find αl, j (z), let us plug the expression αl, j = u j gl + v j ρl into αl,l = 0 and
αl−1,l = 1/l. We obtain
ul =

ρl
,
l[gl−1 ρl − gl ρl−1 ]

vl =

−gl
.
l[gl−1 ρl − gl ρl−1 ]We subtract ρl [zhl gl − (l + 1)gl+1 − lgl−1 ] = 0 from gl [zhl ρl − (l + 1)ρl+1 −lρl−1 ] = 0, and obtain(l + 1) gl (z)ρl+1 (z) − gl+1 (z)ρl (z) = l gl−1 (z)ρl (z) − gl (z)ρl−1 (z)= (l − 1) gl−2 (z)ρl−1 (z) − gl−1 (z)ρl−2 (z)= g0 (z)ρ1 (z) − g1 (z)ρ0 (z) = z.

JOURNAL OF COMPUTATIONAL AND THEORETICAL TRANSPORT

605

Thus, we obtainzαl, j (z) = ρ j (z)gl (z) − g j (z)ρl (z).
Finally, Equation (16) becomesψ̃l = gl (z)ψ̃0 − χl (k),

(19)

where
l



ρl (z)g j (z) − gl (z)ρ j (z) S j (k̂) + ρl (z)χl (k) = (1 − δl0 )
j=1

=

l




ρl (z)g j (z) − gl (z)ρ j (z) S j (k̂).

(20)

j=0

To find the initial term ψ̃0 , we return to Equation (13). We obtain1ˆ ˆ 0 ) =  ψ̃0 (k) +ˆ −ˆ 0 ).
G̃(k, ;δ(ˆ
ˆ
4π 1 + ik · 
1 + ik · Thus, we have

 1
Pl (μ)1ˆ 0 ).
dμ +ψ̃l (k) = ψ̃0 (k)Pl (k̂ · 
ˆ02−1 1 + ikμ
1 + ik · 

Setting l = 0, we have
[1 −  L0 (z)] ψ̃0 (k) =where
1
Ll (z) =2

 1

Pl (μ)zdμ =2−1 1 + ikμ1
ˆ0
1 + ik · , 1

Pl (μ)
dμ = zQl (z).
−1 z − μHere, Ql (z) is the Legendre function of the second kind which has a branch cut on[−1, 1]. We obtainψ̃0 (k) =z1,ˆ0(z) z − k̂ · (21)

where we used
1 −  L0 (z) = 1 −  zQ0 (z) = 1 −1
z z + 1ln= 1 −  z tanh−12z−1z= (z).

The function (z) is defined in Equation (6). We can calculate ψ̃l (k) using Equations (19) and (21).
Equation (14) implies that the Fourier transform of the angular flux is given byˆ =ψ̃ (k, )∞
2l + 1l=04π

ψ̃l (k) Rk̂ Pl (μ).

(22)

606

M. MACHIDA

When the above equation is rewritten using Equation (19), the dependence of ψ̃0
becomes evident asˆ ψ̃0 (k) − T (z, ),ˆˆ = φ (z, )ψ̃ (k, )
k̂
k̂
where
ˆ =
φk̂ (z, )∞
2l + 1l=0ˆ =Tk̂ (z, )4π∞
2l + 1l=04π

gl (z) Rk̂ Pl (μ),χl (k) Rk̂ Pl (μ).

ˆ
Now, we return to the real space from the Fourier space by inverting ψ̃ (k, ).
From Equations (19), (21), and (22), we obtain

1
ˆ dk
ˆ
eik·r ψ̃ (k, )ψ (r, ) =(2π )3 R3∞ 
2l + 1 
1
ik·r
ˆ
=
e
Y
(
)
R
l0
k̂
(2π )3 R3
4π
l=0


× gl (z)ψ̃0 (k) − χl (k) dk.
(23)
Moreover from (23), we have
∞ 
l 

1
2l + 1ˆ =ˆψ (r, )Ylm ()(2π )34π
l=0 m=−lleik·r e−imϕk̂ dm0(θk̂ ) gl (z)ψ̃0 (k) − χl (k) dk
×
R3


∞
l
1 
ˆ
=
Ylm ()eik·r e−imϕk̂ κlm (k) dk,
3(2π )3R(24)

l=0 m=−lwhereκlm (k) =!
"gl (z)
2l + 1 lzd (θ )− χl (k) .
ˆ0
4π m0 k̂ (z) z − k̂ · 

We can obtain the Green’s function (24) with this approach also in the case ofanisotropic scattering (Machida, 2016).
4. Concluding remarks
Using the simple case of isotropic scattering in a 3D infinite medium, we haveseen how the angular flux is obtained with rotated reference frames. In one dimension, the solution by the singular-eigenfunction approach can be derived from the

JOURNAL OF COMPUTATIONAL AND THEORETICAL TRANSPORT

607

Fourier-transform approach (Ganapol, 2000, 2015). It is an interesting future problem to show the equivalence of Equations (10) and (24).

References
Bareiss, E. H., Abu-Shumays, I. K. (1967). On the structure of the isotropic transport operatorsin three independent space variables. ANL-7328.
Cannon, C. J. (1973). An exact solution to the multi-dimensional line transfer equation. J. Quant.
Spec. Rad. Trans. 13:1011–1021.
Case, K. M. (1960). Elementary solutions of the transport equation and their applications. Ann.
Phys. 9:1–23.
Case, K. M., Hazeltine, R. D. (1970). Three-dimensional linear transport theory. J. Math. Phys.
11:1126–1135.
Case, K. M., Zweifel, P. F. (1967). Linear Transport Theory. Massachusetts: Addison-Wesley.
Cassell, J. S., Williams, M. M. R. (2000). The infinite medium Green’s function in space and angle
for the integro-differential form of the neutron transport equation with isotropic scattering.
Ann. Nucl. Enegy 27:11–19.
Dede, K. M. (1964). An explicit solution of the one velocity must-dimensional Boltzmannequation in PN approximation. Nukleonik 6:267–271.
Ganapol, B. D. (2000). A Consistent Theory of Neutral Particle Transport in an Infinite Medium.
Transp. Theory Stat. Phys. 29:43–68.
Ganapol, B. D. (2015). The Infinite Medium Green’s Function of Monoenergetic Neutron Transport Theory via Fourier Transform. Eccl. Sci. Eng. 180:224–246.
Ganapol, B. D., Kornreich, D. E. (1995). A 3-D neutron transport benchmark solution. Transp.
Theor. Stat. Phys. 24:89–111.
Garrettson, G., Leonard, A. (1970). Green’s functions for multidimensional neutron transport ina slab. J. Math. Phys. 11:725–740.
Gibbs, A. G. (1969). Analytical solutions of the neutron transport equation in arbitrary convex
geometry. J. Math. Phys. 10:875–890.
Ishimaru, A. (1978). Wave Propagation and Scattering in Random Media. New York, NY: Academic.
Kaper, H. G. (1969). Elementary solutions of the reduced three-dimensional transport equation.
J. Math. Phys. 10:286–297.
Kim, A. D. (2004). Transport theory for light propagation in biological tissue. J. Opt. Soc. Am. A
21:820–827.
Kobayashi, K. (1977). Spherical harmonics solutions of multi-dimensional neutron transport
equation by finite fourier transformation. J. Nucl. Sci. Tech. 14:489–501.
Lam, S. K., Leonard, A. (1971). Milne’s problem for two-dimensional transport in a quarter space.
J. Quant. Spec. Rad. Trans. 11:893–904.
Lam, S. K., Leonard, A. (1973). Two-dimensional transport in a quarter space: Uniform-source
problem. Trans. Theory and Stat. Phys. 31:1–28.
Leonard, A. (1971). Two dimensional quarter space problems in one speed transport theory. J.
Math. Phys. 12:754–766.
Liemert, A., Kienle, A. (2011a). Radiative transfer in two-dimensional infinitely extended scattering media. J. Phys. A: Math. Theor. 44:505206.
Liemert, A., Kienle, A. (2011b). Analytical solution of the radiative transfer equation for infinitespace fluence. Phys. Rev. A 83:015804.
Liemert, A., Kienle, A. (2011c). Analytical Green’s function of the radiative transfer radiance forthe infinite medium. Phys. Rev. E 83:036605.

608

M. MACHIDA

Liemert, A., Kienle, A. (2012a). Infinite space Green’s function of the time-dependent radiativetransfer equation. Biomed. Opt. Exp. 3:543–551.
Liemert, A., Kienle, A. (2012b). Light transport in three-dimensional semi-infinite scattering
media. J. Opt. Soc. Am. A 29:1475–1481.
Liemert, A., Kienle, A. (2012c). Green’s functions for the two-dimensional radiative transfer equation in bounded media. J. Phys. A: Math. Theor. 45:175201.
Liemert, A., Kienle, A. (2012d). Analytical approach for solving the radiative transfer equation in
two-dimensional layered media. J. Quant. Spec. Rad. Trans. 113:559–564.
Liemert, A., Kienle, A. (2012e). Green’s function of the time-dependent radiative transport equation in terms of rotated spherical harmonics. Phys. Rev. E 86:036603.
Liemert, A., Kienle, A. (2013a). The line source problem in anisotropic neutron transport withinternal reflection. Ann. Nucl. Energy 60:206–209.
Liemert, A., Kienle, A. (2013b). Exact and efficient solution of the radiative transport equation
for the semi-infinite medium. Sci. Rep. 3:2018.
Liemert, A., Kienle, A. (2013c). Two-dimensional radiative transfer due to curved Dirac delta line
sources. Waves in Random and Complex Media 23:461–474.
Liemert, A., Kienle, A. (2015). A novel analytical solution for the radiance in an anisotropicallyscattering medium. Appl. Opt. 54:1963–1969.
Machida, M. (2016). The Green’s function for the three-dimensional linear Boltzmann equation
via Fourier transform. J. Phys. Math. Theor. 49(17):175001.
Machida, M. arXiv:1507.07259.
Machida, M. (2014). Singular eigenfunctions for the three-dimensional radiative transport equation. J. Opt. Soc. Am. A 31:67–74.
Machida, M. (2015). An FN method for the radiative transport equation in three dimensions. J.
Phys. A: Math. Theor. 48:325001.
Machida, M., Panasyuk, G., Schotland, J. C., Markel, V. A. (2010). The Green’s function for theradiative transport equation in the slab geometry. J. Phys. A: Math. Theor. 43:065402.
Machida, M., Panasyuk, G. Y., Wang, Z.-M., Markel, V. A., Schotland, J. C. (2016). Radiative transport and optical tomography with large datasets. J. Opt. Soc. Am. A 33:551–558.
McCormick, N. J., Kuščer, I. (1966). Bi-Orthogonality Relations for Solving Half-Space Transport
Problems. J. Math. Phys. 7:2036–2045.
Markel, V. A. (2004). Modified spherical harmonics method for solving the radiative transport
equation. Waves Random Media 14:L13–L19.
Mika, J. R. (1961). Neutron transport with anisotropic scattering. Nucl. Sci. Eng. 11:415–
427.
Panasyuk, G., Schotland, J. C., Markel, V. A. (2006). Radiative transport equation in rotated reference frames. J. Phys. A: Math. Gen. 39:115–137.
Pomraning, G. C. (1996). (Weakly) three-dimensional caseology. Ann. Nucl. Energy 23:413–427.
Schotland, J. C., Markel, V. A. (2007). Fourier-Laplace structure of the inverse scattering problem
for the radiative transport equation. Inv. Prob. Imag. 1:181–188.
Schreiner, S. I., McCormick, N. J., Gibbs, A. G. (1969). Transport solution of a quarter-space
problem. Trans. Am. Nucl. Soc. 12:160–161.
Siewert, C. E. (1978). The FN method for solving radiative-transfer problems in plane geometry.
Astrophys. Space Sci. 58:131–137.
Siewert, C. E., Benoist, P. (1979). The FN method in neutron-transport theory. Part I: Theory andapplications. Nucl. Sci. Eng. 69:156–160.
Varshalovich, D. A., Moskalev, A. N., Khersonskii, V. K. (1988). Quantum Theory of AngularMomentum. Singapore: World Scientific.
Williams, M. M. R. (1967). Approximate solutions of the neutron transport equation in two andthree dimensional systems. Nukleonik 9:305–312.

JOURNAL OF COMPUTATIONAL AND THEORETICAL TRANSPORT

609

Williams, M. M. R. (1968). Diffusion length and criticality problems in two- and threedimensional, one-speed neutron transport theory. I. Rectangular coordinates. J. Math. Phys.
9:1873–1885.
Williams, M. M. R. (1982). The three-dimensional transport equation with applications to energy
deposition and reflection. J. Phys. A: Math. Gen. 15:965–983.
Xu, H., Patterson, M. S. (2006a). Determination of the optical properties of tissue-simulating
phantoms from interstitial frequency domain measurements of relative fluence and phasedifference. Opt. Exp. 14:6485–6501.
Xu, H., Patterson, M. S. (2006b). Application of the modified spherical harmonics method tosome problems in biomedical optics. Phys. Med. Biol. 51:N247–N251.