Machida et al.

Vol. 26, No. 5 / May 2009 / J. Opt. Soc. Am. A

1291

Diffusion approximation revisited
Manabu Machida,1,* George Yu. Panasyuk,1 John C. Schotland,1 and Vadim A. Markel2
1

Department of Bioengineering, University of Pennsylvania, Philadelphia, Pennsylvania 19104, USA
2
Departments of Radiology and Bioengineering, University of Pennsylvania, Philadelphia,
Pennsylvania 19104, USA*Corresponding author: mmachida@seas.upenn.edu
Received January 29, 2009; revised March 17, 2009; accepted March 19, 2009;
posted March 23, 2009 (Doc. ID 106888); published April 30, 2009

We study the diffusion approximation (DA) to the radiative transport equation (RTE) in infinite homogeneous
space. Different definitions of the reduced intensity Ir that satisy a simplified RTE (without accounting forscattering) and that are often used in the derivation of the DA are examined. By comparing the results of theDA with exact solutions to the RTE, we come to the conclusion that the best accuracy in the DA is achieved if
we choose the definition of the reduced intensity (from a family of possible definitions) that results in Ir = 0.
Thus, the separation of the specific intensity into reduced and diffuse components is unnecessary. We also discuss the conditions under which the DA is applicable. © 2009 Optical Society of America
OCIS codes: 030.5620, 290.1990, 170.5280.

⑀共r,ŝ兲 = A␦共r兲␦共ŝ − ŝ0兲,

1. INTRODUCTION
The theoretical analysis of problems involving multiple
scattering of waves is frequently based on the radiativetransport equation (RTE) [1–4]. In many practical applications, the diffusion approximation (DA) to the RTE canbe used. The advantage of the DA is its relative mathematical simplicity. It is of particular interest in biomedical imaging of tissues with multiply scattered light [5–7].
The diffusion equation (DE) depends on a number of parameters and functions that are inherited from the RTE.
These include the diffusion and absorption coefficients,
the extrapolation distance (which appears in the boundary conditions) and the source function. While the definition of the absorption coefficient is quite straightforward,
a significant effort has been devoted to obtaining an accurate expression that relates the diffusion coefficient to theparameters of the RTE [8–13]. The boundary conditions
for the DE and the extrapolation distance have been considered in detail [14–16]. The source function of the theDE has also been discussed in [17–20]. The optimal shapeof the source function was deduced by comparison of theanalytical solution to the DE with experimental measurements [17–19] or to the results of Monte Carlo simulations [20]. However, the question was not studied systematically from the theoretical point of view. The main
purpose of this paper is to address this gap.
In the transport theory, the primary physical quantityof interest is the specific intensity I共r , ŝ兲 where r is the
position and the unit vector ŝ specifies the direction inwhich light propagates. In the DA one is interested onlyin the angularly averaged quantities u共r兲 and J共r兲 [defined below in Eqs. (6) and (7)], which have the physicalmeaning of electromagnetic energy density and current.
Correspondingly, the mathematical form of the source
term in the RTE is markedly different from that in theDE. Namely, in the case of RTE, the source term is a function of both r and ŝ and can be written as ⑀共r , ŝ兲. For example, a point unidirectional source of the form
1084-7529/09/051291-10/$15.00

共1兲

describes a narrow beam of light of total power A injected
into a scattering medium in the direction of ŝ0. Without
loss of generality, we assume that the injection point coincides with the origin of the laboratory frame. In the caseof the DE, the source term is a function of position only
and is denoted by S共r兲. Obviously, a transition from the
RTE with a point unidirectional source of the form (1) tothe DE requires that S共r兲 cannot be spherically symmetric. This deviation from spherical symmetry is required tocapture the dependence of ⑀共r , ŝ兲 on the direction ŝ.
The above fact has been widely recognized in the literature. The frequently accepted approach is to use thesource function for the DE in the form of a “dipole.” Forexample, in [17], the source was assumed to be constantinside a half-sphere. In this geometry, the source is characterized not only by location but also by direction. In[18,19], the problem was considered in which a narrowcollimated beam is injected into a highly scattering medium from vacuum. The source function was modeled inthese references as a point that is displaced from theboundary into the medium by a distance ⌬x, which must
be determined experimentally. In [20], the DA (in the 1D
geometry) was interpreted as the P1 approximation inwhich the higher angular moments of the RTE source
function (1) appear naturally. Yet another commonly usedalternative is to assume that S共r兲 is an exponentially decaying function along the ray that points from the origin
in the direction ŝ0 and zero everywhere else [3]. This formof S共r兲 can be obtained by decomposing the specific intensity into its reduced and diffuse components [3,21] Ir共r , ŝ兲and Id共r , ŝ兲, respectively. The DA is then made for the diffuse component. However, this approach contains a stepthat is quite arbitrary. We will see below that a family ofdifferent DAs can be obtained by using different definitions of Ir. The correct DA must be chosen by comparing
the solutions to the DE to those of the RTE at large dis© 2009 Optical Society of America

1292

J. Opt. Soc. Am. A / Vol. 26, No. 5 / May 2009

Machida et al.

tances from the source. We will show that the most logicalchoice is Ir = 0; correspondingly, the DA must be applied tothe total specific intensity. Thus, we find that introductionof the reduced intensity in the context of the DA is unnecessary.
This paper is organized as follows. In Section 2 webriefly review the RTE and define the density and currentof electromagnetic energy. In Section 3, we use the RTE toderive asymptotic expressions for these two quantities. In
Section 4, we propose a systematic method for deriving
the DA. The method accounts for the ambiguity in thedefinition of Ir by introducing a new adjustable parameterthat determines the rate of exponential decay of thesource function on the ray defined above. In Section 5 wecompare the density of electromagnetic radiation computed from the RTE and from the DE. We come to the conclusion that the introduction of the reduced intensity is
not justified and that the most logical choice is Ir = 0. The
DA is then made for the total specific intensity. In thissection, we also discuss the conditions under which the
DA is applicable. Section 6 contains a summary of obtained results. Some mathematical properties of the diffuse mode of the RTE are given in Appendix A.

J共r兲 =

3. DENSITY AND CURRENT FROM THE
RTE
An explicit solution to Eq. (2) with the source (1) can bederived by the method of rotated reference frames [23,24].
The solution can be written as冕d s⬘A共ŝ,ŝ⬘兲 = 1;2

共3兲and ⑀共r , ŝ兲 is the source term of the form (1). The isotropyof space implies that the phase function can be expandedasA共ŝ,ŝ⬘兲 =兺 A Y 共ŝ兲Y 共ŝ⬘兲,
l

lm

lm

*
lm

共4兲

where Ylm共ŝ兲 are the spherical functions (viewed here asfunctions of the unit vector ŝ) and A0 = 1. The first moment of the phase function is the so-called scattering
asymmetry parameter g,
A1 = g = ŝ ·冕d2s⬘ŝ⬘A共ŝ,ŝ⬘兲 艋 1.

共5兲

Several physical models for the higher-order coefficients
Al have been proposed. In the commonly used Henyey–
Greenstein model [22], Al = gl. The density of the electromagnetic energy u共r兲 and the current of energy J共r兲 areexpressed as the zeroth and the first angular moments ofthe specific intensity
1
u共r兲 =c冕

d2sI共r,ŝ兲,共6兲兺 兺 ␹ ⬘共r兲Y 共ŝ;r̂兲Y ⬘ 共ŝ ;r̂兲.

I共r,ŝ兲 = Am=−⬁ l,l⬘=兩m兩

m
ll

lm

l m

0

共8兲

In this expansion, Ylm共ŝ ; r̂兲 are the spherical functions defined in a reference frame whose z-axis coincides with thedirection of the unit vector r̂. A detailed definition of thesefunctions is given in [23,24], but is not needed here. Indeed, to compute the density and the current according toEqs. (6) and (7), we require only the two integrals冕冕

A共ŝ,ŝ⬘兲I共r,ŝ⬘兲d2s⬘ + ⑀共r,ŝ兲, 共2兲

where ␮t = ␮a + ␮s is the total attenuation coefficient; ␮aand ␮s are the absorption and scattering coefficients, respectively; A共ŝ , ŝ⬘兲 is the phase function normalized by
the condition⬁⬁

We consider the RTE in an infinite, spatially uniform, isotropic medium, which is of the form

冕共7兲

d2sŝI共r,ŝ兲,

where c is the average speed of light in the medium.

2. RADIATIVE TRANSPORT EQUATION

共ŝ · ⵜ + ␮t兲I共r,ŝ兲 = ␮s冕

d2sYlm共ŝ;r̂兲 = 冑4␲␦l0␦m0 ,d2sŝYlm共ŝ;r̂兲 = r̂

冑

4␲
3

共9兲

␦l1␦m0 .

共10兲m

The expansion coefficients ␹ll 共r兲 are given by the ex⬘pressionm␹ll⬘共r兲 =共− 1兲m¯l兺 共− 1兲2␲冑␴l␴l⬘ M=−l¯
⫻

兺⬘

M

具l兩yn共M兲典具yn共M兲兩l⬘典
␭n3 共M兲n¯l⫻兺 C ⬘⬘j=0

兩l−l 兩+2j,0 兩l−l⬘兩+2j,0l,M,l ,−M Cl,m,l⬘,−m k兩l−l⬘兩+2j

冉 冊
r

␭n共M兲

, 共11兲

where l̄ = min共l , l⬘兲 and various quantities appearing in
Eq. (11) are explained below in Eqs. (12)–(15). First, thecoefficients ␴l 共l = 0 , 1 , 2 , . . . 兲 are given by␴l = ␮a + ␮s共1 − Al兲,

共12兲

where Al are the expansion coefficients in Eq. (4). Second,␭n共M兲 and 兩yn共M兲典 are the eigenvalues and the eigenvectors of a set of infinite, tridiagonal, symmetric matrices
B共M兲 (parameterized by the integer M) whose elements
are of the form
具l兩B共M兲兩l⬘典 =where l , l⬘ 艌 兩M兩 andbl共M兲␦l⬘,l−1 + bl+1共M兲␦l⬘,l+1

冑␴ l ␴ l ⬘

,

共13兲

Machida et al.

Vol. 26, No. 5 / May 2009 / J. Opt. Soc. Am. A

bl共M兲 =
m3
Third, Cjj3m
j m
1

1 2

冑

l2 − M2
4l2 − 1

1293

共14兲

.

are the Clebsch–Gordan coefficients.

2

Fourth, kl共x兲 = −ilhl共1兲共ix兲 are the modified spherical Bessel

functions of the second kind (defined here without the ␲ / 2
factor). Finally, the spectrum of eigenvalues of all matrices B共M兲 is symmetric. That is, for each eigenvalue ␭,
there is also an eigenvalue −␭, and the summation in Eq.
(11) is carried out only over such indices n that ␭n ⬎ 0.
Thus, the notation 兺⬘n should be understood as兺 ⬘f = 兺 f .
n

共15兲

n

␭n⬎0

n

We now use the integrals (9) and (10) to compute thedensity and the current as
u共r兲 =

J共r兲 =A ⬁兺 冑2l + 1␹ 共r兲P 共ŝ · r̂兲,
0
0l

c l=0

Ar̂ ⬁

l

共16兲

0

冑2l + 1␹1l0 共r兲Pl共ŝ0 · r̂兲,

冑3 兺
l=0

共17兲

where Pl共x兲 are the Legendre polynomials. The functions
00
␹0l
共r兲 and ␹1l
共r兲 are significantly simplified compared tom0
the more general functions ␹ll 共r兲. Thus, for ␹0l
共r兲, we
⬘
have
0
␹0l
共r兲 =

1

兺⬘
2 ␲ 冑␴ ␴

具0兩yn共0兲典具yn共0兲兩l典
␭n3 共0兲

0 l n

冉 冊
r

kl

␭n共0兲

. 共18兲

0
共r兲 is more involved. It is possible
The expression for ␹1l
to show that

0
␹1l
共r兲 =

1
2 ␲ 冑␴ 1 ␴ l−

冋兺−冑2l共l + 1兲
r

⬘

n

兺
n

⬘

具1兩yn共0兲典具yn共0兲兩l典
␭n3 共0兲具1兩yn共1兲典具yn共1兲兩l典
␭n2 共1兲

冉 冊
冉 冊册

kl⬘

term in the summations over n in both Eqs. (18) and (19),namely, the term corresponding to the diffusion mode. We
then obtain
0
␹0l
共r兲 ⬇

r

␭n共0兲

0
␹1l
共r兲 ⬇

r

kl

Fig. 1. Eigenvalues of the matrices B共M兲 for the following parameters: ␮a = 0.03 cm−1, ␮s = 500 cm−1, g = 0.98. These parameters are typical for biological tissues in the near-IR spectral
range. (a) All eigenvalues of B共0兲 versus the eigenvalue number
n. (b) The maximum eigenvalues of the matrices B共M兲 versus M.
In simulations, infinite matrices B共M兲 were truncated so that thesize of each matrix was N = 103.

␭n共1兲

,
共19兲

where kl⬘共z兲 = dkl共z兲 / dz.
We seek an expression for the density and current far
from the source, so that r in the above expressions is sufficiently large. A more precise mathematical formulationof this condition will become clear momentarily. The functions kl共x兲 decay exponentially as exp共−x兲. Let the largest
eigenvalue of the matrix B共0兲 be ␭d and the corresponding
eigenvector be 兩yd典. We will refer to this eigenstate as the
diffusion mode. As was shown in [24], the spectrum of eigenvalues is discrete for ␭ ⬎ 1 / ␮t, and the largest eigenvalue ␭d is in the discrete spectrum. Therefore, there is a
finite gap ⌬␭ between the largest and the second largest
eigenvalues of B共0兲. In many practical situations, this gap
is quite significant, as is illustrated in Fig. 1(a). In addition, ␭d is typically much larger than the maximum eigenvalue of all the matrices B共M兲 where 兩M兩 ⬎ 0. This is illustrated in Fig. 1(b). Thus, at sufficiently large distancesfrom the source, we can neglect the second term in thesquare brackets in Eq. (19). Further, we can keep only one

1

具0兩yd典具yd兩l典

2␲冑␴0␴l␭d3−1

具1兩yd典具yd兩l典

2␲冑␴1␴l

␭d3

冉冊
冉冊
r

kl

kl⬘

␭d

,

共20兲

.

共21兲

r

␭d

The above equalities are exponentially accurate whenexp共−␬兲 Ⰶ 1, where

␬=

r

⌬␭

␭d ␭d − ⌬␭

.

共22兲

It is now possible to state the condition under which
Fick’s law,
J = − D ⵜ u,

共23兲

is applicable. This condition is r / ␭d Ⰷ 1. Indeed, if this inequality holds, we can write approximately ⵜu = r̂⳵u / ⳵rand neglect the tangential derivative of the factors
Pl共ŝ0 · r̂兲 in Eq. (16). The accuracy of this approximation isalgebraic, and the error is of the order of O共r / ␭d兲, while
the error in Eqs. (20) and (21) is of the order ofO关exp共−␬兲兴. Note, however, that for isotropic sources ofthe form ⑀ = A␦共r兲, Fick’s law is exponentially accurate.
Indeed, if the expression (16) is integrated over ŝ0, only
the l = 0 term remains nonzero. Therefore, the equality
ⵜu = r̂⳵u / ⳵r is, in this case, exact.

1294

J. Opt. Soc. Am. A / Vol. 26, No. 5 / May 2009

Machida et al.

The diffusion coefficient D in Eq. (23) can easily befound from Eqs. (16), (17), (20), and (21). As was mentioned above, in evaluating the gradient of Eq. (16), thetangential derivative can be neglected. We then obtain
D=c

具1兩yd典具0兩yd典

冑

␴0
3␴1

␭d .

共24兲

If we also take into account that ␴0 = ␮a and that the characteristic equation for the diffusion mode implies
关b1 / 冑␴0␴1兴具1 兩 yd典 = ␭d具0 兩 yd典 with b1 = 1 / 冑3, we arrive at theresultD = c␮a␭d2 .

共25兲

Thus, the diffusion coefficient is defined by the average
speed of light in the medium, the absorption coefficient,and the largest eigenvalue of B共0兲. Some mathematicalproperties of the diffusion eigenvalue ␭d and of the corresponding eigenvector 兩yd典 are given in Appendix A.
Introducing the diffuse wave number kd = 1 / ␭d and using Eqs. (16) and (20), we can write the density u共r兲 in theform
u共r兲 = A

kd

兺S

D lm

共RTE兲
* 共ŝ 兲,kl共kdr兲Ylm共r̂兲Ylm
0
l

共26兲

where
Sl共RTE兲 =冑␴0
共2l + 1兲␴l

2具0兩yd典具yd兩l典.

共27兲

共ŝ · ⵜ + ␮t兲Ir共r,ŝ兲 = ⑀共r,ŝ兲.

共31兲In highly scattering media such that ␮t Ⰷ ␮a, the above
two definitions result in reduced intensities that decay
with very different exponential rates.
In the standard approach to the DA, however, neither
Eq. (30) nor Eq. (31) is used. Instead, the reduced intensity is defined by
共ŝ · ⵜ + ␮*兲Ir共r,ŝ兲 = ⑀共r,ŝ兲,

共32兲

where ␮* = 1 / ᐉ* = ␮a + 共1 − g兲␮s is the reciprocal of thetransport mean free path ᐉ*. We note the obvious inequality ␮a 艋 ␮* 艋 ␮t. The definition (32) utilizes a more physically relevant quantity compared with Eqs. (30) and (31),namely, the reduced scattering coefficient ␮⬘s = 共1 − g兲␮s.
Thus Eq. (32) predicts that the reduced intensity in a medium with strictly forward-peaked scattering 共g = 1兲 is nodifferent from that in a medium with no scattering at all
共␮s = 0兲, as one could expect on physical grounds. The deficiency of the definitions (30) and (31) is that the first iscompletely independent of ␮s while the second is independent of g.
In spite of the above, all three definitions (30)–(32) are
ad hoc and need justification. Such justification can be obtained only by direct comparison with the RTE. This will
be done below. However, there is no compelling reason to
restrict consideration to the three discrete values of thecoefficient that enters the definition of reduced intensity(either ␮a, ␮t, or ␮*). Instead, we will pose a more general
problem. Suppose that Ir is defined by¯ 兲Ir共r,ŝ兲 = ⑀共r,ŝ兲,
共ŝ · ⵜ + ␮

共33兲

In particular,
S0共RTE兲 = 2具0兩yd典具yd兩0典,

S1共RTE兲 = 2␮a␭d具0兩yd典具yd兩0典.
共28兲

We will see below that a similar expansion (but with dif共DA兲
ferent coefficients Slm) can be obtained in the DA.

4. DIFFUSION APPROXIMATION

共ŝ · ⵜ + ␮t兲Id = ␮sAId + ⑀r ,The diffusion approximation (DA) to the RTE (2) is usually obtained as follows. We first expand the specific intensity as
I共r,ŝ兲 = Ir共r,ŝ兲 + Id共r,ŝ兲,

共29兲

where Ir is the reduced intensity, which is defined to satisfy a reduced RTE, and Id is the “diffuse” component of
the specific intensity. It is usually assumed that Id obeys(approximately) a DE, while Ir is highly singular and,
therefore, must be considered separately [3].
In this paper, we point out that the reduced intensity
can be defined in more than one way. Thus, for example,
we can require that the reduced intensity satisfy the RTE(2) in which we set ␮s = 0. This leads to the equation
共ŝ · ⵜ + ␮a兲Ir共r,ŝ兲 = ⑀共r,ŝ兲.

¯ is an arbitrary positive constant (additional conwhere ␮¯ will be imposed later). We will
straints on the values of ␮¯ as a free parameter and compare the results obview ␮¯ with thetained by making the DA for a given value of ␮
predictions of radiative transport theory.
We now proceed with the derivation of the DA. First,
we substitute the decomposition (29) into Eq. (2) and obtain the equation for the diffuse component,

共30兲

We can also formally set A共ŝ , ŝ⬘兲 = 0 in Eq. (2) and obtain

共34兲

where the reduced source function ⑀r is given by¯ 兲Ir⑀r共r,ŝ兲 = 共␮sA − ␮t + ␮

共35兲

and Ir satisfies Eq. (33). In the above two equations, A denotes the linear operator defined by the integral in theright-hand side of Eq. (2).
We seek an approximate solution to Eq. (34) in the form
Id共r,ŝ兲 = 共c/4␲兲ud共r兲 + 共3/4␲兲ŝ · Jd共r兲.

共36兲

It can be verified by direct integration of Eq. (36), according to Eqs. (6) and (7), that ud and Jd are the energy density and current associated with the diffuse intensity Id.
The DA is obtained by substitution of the ansatz (36)
into Eq. (34) and by considering the first two moments(with respect to the angular variable) of the resultant
equation. The substitution yields

Machida et al.

Vol. 26, No. 5 / May 2009 / J. Opt. Soc. Am. A

共ŝ · ⵜ + ␮a兲ud + 3/c共ŝ · ⵜ + ␮*兲ŝ · Jd = 共4␲/c兲⑀r .

共37兲

We now evaluate the zeroth and the first moments of theabove equation with respect to ŝ. In the first case, we integrate Eq. (37) over d2s and, in the second case, overŝd2s. This leads to the two equations
ⵜ · Jd + c␮aud = E,

共38兲

ⵜud + 共3␮*/c兲Jd = 共3/c兲Q,

共39兲

where the scalar and the vector source terms E and Q are
given by the expressions
E共r兲 =

Q共r兲 =

冕冕d2s⑀r共r,ŝ兲,共40兲d2sŝ⑀r共r,ŝ兲.

共41兲

We can use Eq. (35) to express the above two functions interms of the reduced intensity Ir:¯ − ␮ a兲
E共r兲 = 共␮¯ − ␮ *兲
Q共r兲 = 共␮冕冕

d2sIr共r,ŝ兲,共42兲d2sŝIr共r,ŝ兲.

共43兲¯ = ␮* is special: it
It can now be seen why the choice ␮causes Q共r兲 to vanish. We will show, however, that accounting for this term in the DE does not lead to additional difficulties.
At this point, we require a specific expression for the
reduced intensity. It can be easily seen that the solution
to Eq. (33) with the source (1) is
Ir共r,ŝ兲 = A¯ r兲
exp共− ␮r2

␦共ŝ − ŝ0兲␦共r̂ − ŝ0兲.

共44兲

We then substitute Eq. (44) into Eqs. (42) and (43) and obtain¯ − ␮ a兲
E共r兲 = A共␮¯ r兲exp共− ␮¯ − ␮ *兲
Q共r兲 = Aŝ0共␮r2

␦共r̂ − ŝ0兲,¯ r兲
exp共− ␮r2

␦共r̂ − ŝ0兲.

共45兲

共46兲

1295

The diffusion coefficient (48) and the similar quantity(25) obtained from the RTE under more general conditions differ, but approach each other in the limit ␮a / ␮s
→ 0. In this limit, ␭d is sharply bounded by the inequalities (A8) of Appendix A. If only the lowest-order (in ␮a / ␮s)non-vanishing term is retained, ␭d ⬇ 1 / 冑3␮a␮*. Substitution of this expression into Eq. (25) results in D ⬇ c / 3␮* inagreement with Eq. (48). However, the solution to the DE(47) depends on D exponentially. Therefore, even a small
error in D can result in an exponentially large discrepancy between the solution to the DE and the solution to
the RTE. It has been suggested that the DA can be “corrected” by using a more accurate expression (25) for the
diffusion coefficient [11–13]. This correction does not follow in a mathematically consistent way from the DA.
Nevertheless, it has proved to be useful. We will adopt
this approach below. Specifically, we will use the expression (25) instead of Eq. (48) for the coefficient D that appears in the DE (47).
The second term in (49) contains a derivative of a delta
function, namely, ⵜ · Q. Although this has not been statedexplicitly in the literature, it seems plausible that the¯ = ␮* in the definition of Ir has been made because
choice ␮
it results in Q = 0. However, the presence of the term ⵜ · Q
does not pose a mathematical problem. To evaluate thelatter, we act with the operator ⵜ on Eq. (43) and use Eq.
(33). A straightforward calculation yields冋¯ − ␮*兲 ␦共r兲 − ␮¯
ⵜ · Q = A共␮¯ r兲exp共− ␮r2

册

␦共r̂ − ŝ0兲 . 共50兲Combining this expression with Eq. (45), we obtain

再S共r兲 = A −¯ − ␮*
␮␮*

␦共r兲 +¯ 2 − ␮a␮* exp共− ␮
¯ r兲
␮

␮*

r2

冎

␦共r̂ − ŝ0兲 .
共51兲

It can be seen that the source term contains two contributions. The first term describes a point source located at¯ = ␮*. The second term isthe origin; this term vanishes if ␮
a function that exponentially decays along the ray r̂ = ŝ0;¯ = 冑␮a␮*. The integral source isthis term vanishes if ␮冕¯ 兲.
d3rS共r兲 = A共1 − ␮a/␮

共52兲is the diffusion coefficient obtained in the approximation
specified by Eqs. (29) and (36), and the source term isgiven by

We expect on physical grounds the above expression to be¯ must satisfy ␮¯ ⬎ ␮a. It also
positive and conclude that ␮
can be expected that the integral (52) should be equal toA, which is the overall power of the source. This already¯ is ␮¯ = ⬁. We will
suggests that the proper choice for ␮
come to the same conclusion using more rigorous arguments below. At this point, we note that under the conditions when the DA is typically used, ␮a Ⰶ ␮* and the choice¯ = ␮* may result in a relatively small error. Nevertheless,␮
this error is both non-negligible and easily avoidable. We¯ = ⬁ results in a more natuwill also show that the choice ␮
ral and transparent theory.
Note that the current is given in the DA by the formula

S共r兲 = E共r兲 − ᐉ* ⵜ · Q共r兲.

Jd共r兲 = − D ⵜ ud共r兲 + ᐉ*Q共r兲.

At the next step, we exclude the current from Eqs. (38)and (39) and obtain a DE with respect to the density
ud共r兲, namely,
共− Dⵜ2 + c␮a兲ud共r兲 = S共r兲,

共47兲

D = c/3␮*

共48兲

where

共49兲

共53兲

1296

J. Opt. Soc. Am. A / Vol. 26, No. 5 / May 2009

Machida et al.

We now use the above results to compute the density
ud共r兲 due to the source function S共r兲 (51). We have
ud共r兲 =

冕d3r⬘G共r,r⬘兲S共r⬘兲,

共54兲

where
G共r,r⬘兲 =exp共− kd兩r − r⬘兩兲
4␲D兩r − r⬘兩

共55兲

is the free-space Green’s function of the DE (47) andkd = 冑c␮a/D = 1/␭d

G共r,r⬘兲 =兺

D lm* 共r̂⬘兲,
il共kdr⬍兲kl共kdr⬎兲Ylm共r̂兲Ylm

ud共r兲 =

兺 Y 共r̂兲

D lm
⫻

冕

lm

冕

kd

兺S

D lm

共DE兲
* 共ŝ 兲,kl共kdr兲Ylm共r̂兲Ylm
0
l

共59兲

where
Sl共DE兲 = −dr⬘共r⬘兲 il共kdr⬍兲kl共kdr⬎兲2␮*Ql共p兲 =

␦l0 +¯ 2 − ␮ a␮ *␮␮ *k d¯ /kd兲.
Q l共 ␮

共60兲冕

⬁il共x兲exp共− px兲dx,

共61兲

0

is the Legendre function of the second kind, and we haveused the specific form of the source function (51). Notethat the Legendre function must be evaluated for arguments that are larger than unity. In particular,
Q0共p兲 =⬁¯ − ␮*
␮

Here

共57兲

where iᐉ共x兲 and kᐉ共x兲 are the modified spherical Bessel
functions of the first and second kind, respectively, and r⬍
and r⬎ are the lesser and greater of r and r⬘. Note thatkᐉ共x兲 is defined here without the ␲ / 2 factor, so that k0共x兲
= exp共−x兲 / x, etc. From this, we obtainkd

ud共r兲 = A

共56兲is the diffuse wave number. Note that, as discussed above,
we use the expression (25) for the diffusion coefficient.
It is convenient to expand the Green’s function as
kd

slower than the source function S共r兲 when r → ⬁. This¯ ⬎ kd which is in addition to theleads to the requirement ␮¯ ⬎ ␮a. If we assume
previously imposed requirement ␮that the point of observation r is sufficiently far from theorigin, the radial integral in Eq. (58) converges while r⬍= r⬘ and r⬎ = r. The density can then be written with exponentially high precision in the form (26), namely, as

p+1

1
2

ln

p−1

,

Q1共p兲 = − 1 +

p+1

p
2

ln

p−1

.

共62兲

0

* 共r̂⬘兲,d2r̂⬘S共r⬘兲Ylm

共58兲

We now note the following. The Green’s function (55) decays exponentially as exp共−kd 兩 r − r⬘ 兩 兲, while the the¯ r⬘兲. On physicalsource function (51) decays as exp共−␮
grounds, we expect that the density ud共r兲 should decay

5. RESULTS AND DISCUSSION
We now compare the results derived in Sections 3 and 4.
In Fig. 2, we plot the absolute difference between the expansion coefficients Sl共RTE兲 that are defined in Eq. (27) and
Sl共DE兲 that are defined in Eq. (60). Comparison is madeonly for l = 0 and l = 1. The higher-order coefficients Sl共RTE兲¯ / ␮* for l = 0 , 1 and for the following sets of parameters: (a) ␮a = 0.03 cm−1, ␮sFig. 2. Absolute difference Sl共DE兲 − Sl共RTE兲 as a function of ␮
= 500 cm−1, g = 0.98 (these are the “physiological” parameters typical of biological tissues in the near-IR spectral range; same parameters
have been used in Fig. 1). (b) ␮a = 1 cm−1, ␮s = 5 cm−1, g = 0 (isotropic scattering). (c) ␮a = 1 cm−1, ␮s = 10 cm−1, g = 0.5 [anisotropic scattering
but same ␮* as in (b)]. (d) ␮a = ␮s = 1 cm−1, g = 0.5.

Machida et al.

Vol. 26, No. 5 / May 2009 / J. Opt. Soc. Am. A

depend on the higher-order coefficients Al (in addition tog = A1). The DA is independent of these higher-order coefficients and, therefore, direct comparison for l ⬎ 1 is inappropriate.
¯ = ␮* has no
It can be seen that in all cases, the point ␮
共DE兲
共RTE兲
special significance. The difference S0 − S0
is ratherflat, especially in the three cases (a)–(c) for which the DA
is typically considered to be applicable. If only the zerothorder expansion coefficients are examined, one can con¯ does not matter at all, and one
clude that the choice of ␮
¯ , procan choose without loss of accuracy any value of ␮¯ ⲏ ␮*. However, the first-order difference
vided that ␮¯ . This sugS1共DE兲 − S1共RTE兲 monotonically decreases with ␮¯ = ⬁. It may seem thatgests that the optimal choice is ␮
this choice is paradoxical and mathematically ill defined.
However, this is not so. The quantity ud共r兲 is well defined¯ → ⬁. Indeed, we havein the limit ␮
lim S0共DE兲 = 1,¯ →⬁␮lim S1共DE兲 = kd/3␮*,¯ →⬁␮

Sl共DE兲 = 0 for l ⬎ 1.
共63兲

From this, we findud共r兲 = 共1 + ᐉ*kdr̂ · ŝ0兲exp共− kdr兲/4␲Dr= 共1 + ᐉ*kdr̂ · ŝ0兲G共r,0兲.

共64兲¯ → ⬁ limit of Eq. (54)Alternatively, we can take the ␮with S共r兲 given by Eq. (51). We obtainlim S共r兲 = A共1 − ᐉ*ŝ0 · ⵜ兲␦共r兲,共65兲

lim ud共r兲 = A共1 + ᐉ*ŝ0 · ⵜr⬘兲G兩共r,r⬘兲兩r⬘=0 .

共66兲

¯ →⬁
␮

¯ →⬁
␮

The same result can be obtained even more directly by¯ → ⬁ corresponds to Ir = 0. We thus
noting that the limit ␮
can apply the DA to Eq. (34) in which the reduced source⑀r must be substituted by the original source of the RTE,
⑀, given in Eq. (1). This yields E共r兲 = A␦共r兲, Q共r兲
= Aŝ0␦共r兲 and, according to Eq. (49), S共r兲 = A共1− ᐉ*ŝ0 · ⵜ兲␦共r兲. This is in agreement with Eq. (65). Thus, it
can be seen that the use of the reduced intensity Ir is unnecessary. Note that its introduction actually decreases
the accuracy of the DA. Also, the reduced intensity cannot
be interpreted as the correct specific intensity in the vicinity of the source. This is clear already from the ambiguity in the definition of Ir, as discussed in Section 4.
We now compare the density u共r兲 computed using theRTE [according to Eq. (16)] without any approximations
and the density ud共r兲 obtained in the DA. The results areshown in Fig. 3. The relative error ⌬, defined as⌬=兩u − ud兩
u

,

共67兲

¯ for various distances to theis plotted as a function of ␮
source (r = 10ᐉ* and r = 20ᐉ*) and various orientations of
the point of observation relative to the direction of inci¯ = ␮* is in no way
dence ŝ0. It can be seen that the point ␮
special or optimal. Depending on the angle ␪= arccos共r̂ · ŝ0兲, the discrepancy has minima at certain¯ . None of these values canseemingly random values of ␮1297

be universally used because each of them minimizes theerror only for a specific distance r, a specific angle ␪, anda specific set of parameters ␮a, ␮s, and g. At the same¯ = ⬁ is almost always
time, it is obvious that the choice ␮¯ = ␮ *.
preferable to ␮¯ = ␮* inTo conclude this section, we discuss the choice ␮
more detail. We will denote the source function that corresponds to this case by S*共r兲 and the corresponding moments Sl共DE兲 by Sl*. As can be seen from Eqs. (51) and (60),
S*共r兲 = A␮⬘s

exp共− ␮*r兲
r2

␦共r̂ − ŝ0兲,Sl* = 共␮⬘s /kd兲Ql共␮*/kd兲,

共68兲

共69兲

where ␮⬘s = 共1 − g兲␮s is the reduced scattering coefficient. In
biomedical optics, it is typical to work in the regime inwhich ␮* Ⰷ kd, although this inequality is never verystrong. Thus, for the physiological parameters used in
Figs. 1, 2(a), and 3(a), ␮* / kd ⬃ 10. Nevertheless, for thisratio of ␮* / kd, the asymptotic formulas for Q0共␮* / kd兲 andQ1共␮* / kd兲 are already very accurate and we obtainS0* ⬇ ␮⬘s /␮*,S1* ⬇ 共␮⬘s /␮*兲共kd/3␮*兲.

共70兲

These expressions differ from the first two equations in
Eqs. (63) by the constant factor ␮⬘s / ␮*.
We can make an additional approximation and assume
that all higher-order coefficients Sl* 共l ⬎ 1兲 are zero. Thenit turns out that the density ud obtained in the DA with¯ = ␮* differs from the density obtained for ␮¯ = ⬁ (which we␮deem here to be more accurate) by the overall factor of␮⬘s / ␮*. In the case of physiological parameters, this factor
is close to unity. Essentially, this approach was adopted
by us earlier in [25] where the factor ␮⬘s / ␮* = ␮⬘s ᐉ* appears,
for example, in Eq. (16). Note that in [25], we have implicitly assumed that Sl* = 0 for l ⬎ 1 by truncating the Taylor
expansion of the Green’s function G共r , r⬘兲 at the firstorder, namely, by writing G共r , r⬘兲 ⬇ G共r , r0兲 + 共r⬘− r0兲兩 · ⵜr⬘G共r , r⬘兲兩r⬘=0, where r0 is the location of thesource. This approximation was then used in the integralof the form (54).
¯ = ␮* and obThus, in [25], we have utilized the choice ␮tained the result for ud the differs from Eq. (65) or (66) by
the overall factor ␮⬘s / ␮*, which may seem to be an insignificant correction. However, an additional approximation
was used in [25] (neglecting the higher-order moments Sl*
as discussed above) that was not well justified. In a sense,
there were two inaccuracies in [25] that canceled eachother. The first inaccuracy was neglecting the higherorder moments, and the second inaccuracy was the incor¯ . Further, in less ideal cases, the ratio
rect choice of ␮␮⬘s / ␮* may be significantly different from unity, yet the
DA could still be applicable sufficiently far from the¯ can lead tosource. In such cases, an incorrect choice of ␮
significant errors. Finally, the derivation of the DA presented in this paper is mathematically consistent andbased on direct comparison with the RTE rather than on¯.
an ad hoc choice of ␮

1298

J. Opt. Soc. Am. A / Vol. 26, No. 5 / May 2009

Machida et al.

¯ / ␮* for different values of the distance from the source r and the angle ␪, where cos ␪
Fig. 3. Relative error ⌬ (67) as a function of ␮
= ŝ0 · r̂. The distance to the source is r = 共a兲 10ᐉ*, (b) 20ᐉ*. All parameters are the same as in Fig. 1 and Fig. 2(a).

So far, we have not discussed the conditions underwhich the DA is applicable. This question is both simple
and difficult. It is often remarked that the DA is valid
when ␮a Ⰶ ␮s. However, it is not always clear what exactly
is meant by applicability of the DA.
If it is deemed that the DA is applicable if the exact
density of electromagnetic energy u defined by Eq. (6) satisfies (approximately) the diffusion equation (47), thenthe inequality ␮a Ⰶ ␮s is not really required (for the reminder of this section, the subscript d in the referenced
formulas should be omitted). Indeed, the diffusion equation for u follows from the continuity equation (38) andFick’s law (23). But the continuity equation that couples
the density and the current is exact irrespective of parameters. As for Fick’s law, the condition for its applicability
was discussed at the end of Section 3 and reads r Ⰷ ␭d.
Thus, at sufficiently large distances from the source,
Fick’s law is accurate irrespective of the ratio ␮a / ␮s. It
should be kept in mind though that when the ratio ␮a / ␮sis not small, expressions (25) and (48) can yield very different values of the diffusion coefficient, and the diffusion
constant appearing in Fick’s law is given by the former
expression.

If, however, we deem the DA to be accurate when thespecific intensity I共r , ŝ兲 is well approximated by the expansion (36), then the applicability of the DA is more difficult to show. Indeed, for the expansion (36) to be accum
rate, it is required that all the coefficients ␹ll 共r兲 in the⬘
expansion (8) be negligibly small for l , l⬘ ⬎ 1. Now conm
sider the expression (11) for ␹ll 共r兲. If we make an approxi⬘
mation in which only one term in the summation over n isretained, namely, the term that corresponds to the diffusion mode (see Appendix A), the above condition can not
be proved. Indeed, the Bessel functions kl共x兲 increase factorially with l (for large values of l and a fixed x), while
the components of the diffuse eigenvector 具l 兩 yd典 decrease
with l only exponentially, as is shown in Appendix A.
Thus, the factorial growth of the first factor always overpowers the exponential decay of the second. In fact, if an
accurate pointwise approximation to I共r , ŝ兲 is beingsought, it is not correct to retain only the diffuse mode in
Eq. (11). Instead, summation over all modes must be performed, including the modes of the continuous spectrum.
The terms in this summation have different sign. As a result, if the summation is performed with sufficient nu-

Machida et al.

Vol. 26, No. 5 / May 2009 / J. Opt. Soc. Am. A

1299

␤2 = 2/冑15␮*关␮a + ␮s共1 − A2兲兴,

共A3兲

m

merical accuracy, the resultant coefficients ␹ll 共r兲 decrease⬘
with l and l⬘ so that the expansion (8) is, in fact, convergent. (When the summation is performed numerically, one
should be conscious of the round-off errors, which can become large and even dominant for very large values of l
and l⬘.) Therefore, the proof of pointwise convergence ofthe specific intensity to the “diffuse” value given by formula (36) is difficult to obtain, and we are not certain
whether it is known. Such a proof could be, however,purely of academic interest. The reason is that all physical detectors always measure the specific intensity integrated over a finite area and a finite solid angle. This integration regularizes the formula (11) and can greatly
improve numerical convergence. It therefore can be statedthat the condition of validity of the DA depends on thetype of detector used.

6. SUMMARY
We have systematically developed the diffusion approximation (DA) to the radiative transport equation (RTE) ininfinite homogeneous space. We have examined different
admissible definitions of the reduced intensity Ir that iscommonly used in the derivation of the DA. The existing
ambiguity in the definition of Ir affects the form of thesource function for the DE. By comparing the results of
the DA with more rigorous solutions to the RTE, we havefound that the best accuracy is achieved if we set Ir = 0. In
this case the source function for the DE is given by Eq.
(65) and the density of electromagnetic energy is given by
Eq. (66) or, equivalently, by Eq. (64). We conclude that the
separation of the total specific intensity into the reduced
and the diffuse components is not justified and, in a typical derivation of the DA, leads to additional errors. Thus,
the DA should be made for the total specific intensity.
The theory developed in this paper is not intended tosignificantly improve the accuracy of the DA, althoughsome improvement can be anticipated. It is rather aimed
at removing the existing ambiguities and formulating the
DA in a mathematically consistent way. We also study thegeneral limits of applicability of Fick’s law and derive certain mathematical properties of the so-called diffuse modeof the RTE.

APPENDIX A: DIFFUSION MODE
In this appendix, we derive some mathematical propertiesof the diffusion mode, more specifically, of the eigenvalue
␭d and of the corresponding eigenvector 兩yd典. By definition, ␭d is the largest eigenvalue of the infinite, tridiagonal, real and symmetric matrix B共0兲 that has been defined in Section 3, Eq. (13). We denote the elements of thefirst superdiagonal of this matrix by ␤ᐉ, where␤l = 具l − 1兩B共0兲兩l典 =bl冑␴l−1␴l,l = 1,2, . . . .

共A1兲

The coefficients bl are given in Eq. (14) in which we mustset M = 0, so that bl = l / 冑4l2 − 1. Note that

␤1 = 1/冑3␮a␮* ,共A2兲

lim ␤l = 1/2␮t .

共A4兲

l→⬁

In obtaining the limit (A4), we have assumed thatliml→⬁ Al = 0 and liml→⬁ ␴l = ␮t.
We start by deriving the lower and upper bounds of ␭d.
The upper bound is obtained from the Gershgorin theorem, which states that兩␭d兩 艋 maxl兺 共1 − ␦ ⬘兲兩具l兩B共0兲兩l⬘典兩 = ␤ + ␤ .
ll

1

共A5兲

2

l⬘

Here we have assumed that ␴0 ⬍ ␴1 ⬍ ␴2¯, as is typicallythe case, and therefore ␤1 ⬎ ␤2 ⬎ ␤3. . ..
The lower bound can be derived by considering the sequence of truncated matrices Bj共0兲 共j = 1 , 2 , . . . 兲, which are
obtained by keeping only the first j rows and columns of
the infinite matrix B共0兲. As a direct consequence of the
Cauchy interlace theorem, we can state that the eigenvalues of Bj共0兲 and Bj+1共0兲 interlace. Consequently, the maximum eigenvalue of each matrix Bj共0兲, ␭dj, increases monotonically with j. We thus have ␭d1 ⬍ ␭d2 ⬍ ¯ ⬍ ␭d where␭d = limj→⬁ ␭dj. For j = 3,
␭d3 = 冑␤12 + ␤22 .

共A6兲

Thus, we have proved the inequality冑␤12 + ␤22 艋 ␭d 艋 ␤1 + ␤2 .

共A7兲Substituting ␤1 and ␤2 from Eqs. (A2) and (A3), we canrewrite this as冑1 + ␩
1 + 冑␩艋 ␭d 艋
,冑3 ␮ a ␮ *冑3 ␮ a ␮ *

共A8兲

where

␩=

␮a

4

5 ␮a + ␮s共1 − A2兲

共A9兲

.

In the limit when ␮a Ⰶ ␮s, ␩ can be viewed as a small parameter.
A sharper estimate can be obtained if we use ␭d4 as thelower bound for ␭d. The former is given by the formula
␭d4 =

冑

␤12 + ␤22 + ␤32
2

冋 冑
1+

1−

4␤12␤32
共␤12 + ␤22 + ␤32兲2

册

.

共A10兲
Further, the gap between the largest and the second
largest eigenvalues of B共0兲 can be estimated as follows.
Let the second largest eigenvalue be ␭s. From the Gershgorin theorem we have ␭s 艋 max共␤1 , ␤2 + ␤3兲. Combining
this with Eq. (A7), we find that⌬␭ ⬅ ␭d − ␭s 艌 max关冑␤12 + ␤22 − max共␤1, ␤2 + ␤3兲,0兴.
共A11兲
Note that in the limit of strong scattering we have ␤1
Ⰷ ␤2 , ␤3, and the first number in the square brackets is
positive, so that the “max” can be omitted and we have

1300

J. Opt. Soc. Am. A / Vol. 26, No. 5 / May 2009⌬␭ 艌 冑␤12 + ␤22 − ␤1 .

Machida et al.

共A12兲
9.

Finally, we consider the components of the eigenvector
兩yd典. The first two components (l = 0 and l = 1) can be expanded in powers of ␩ as

10.

具0兩yd典 = 1/冑2关1 − 1/2␩ + O共␩2兲兴,

共A13兲

11.

具1兩yd典 = 1/冑2关1 + O共␩2兲兴.

共A14兲

12.

At large values of the index l, the components 具l 兩 yd典 decay
exponentially [24] as

13.

具l兩yd典 ⬀ exp共− ␶l兲,

共l → ⬁兲.

共A15兲

By considering the three-term recurrence relation for具l 兩 yd典 in the limit l → ⬁ when ␤l → 1 / 2␮t, we find

␶ = ln关␮t␭d + 冑共␮t␭d兲2 − 1兴.

共A16兲

14.
15.

16.

ACKNOWLEGMENT
This research was supported by the National Science
Foundation (NSF) under grants EEC-0615857 and DMS0554100.

REFERENCES
1.
2.
3.
4.

5.
6.

7.
8.

S. Chandrasekhar, Radiative Transfer (Dover, 1960).
V. V. Sobolev, Light Scattering in Planetary Atmospheres,Vol. 76 of International Series of Monographs in Natural
Philosophy (Elsevier, 1975).
A. Ishimaru, Wave Propagation and Scattering in RandomMedia (IEEE Press, 1997).
M. C. W. van Rossum and Th. M. Nieuwenhuizen,“Multiple scattering of classical waves: microscopy,
mesoscopy, and diffusion,” Rev. Mod. Phys. 71, 313–371(1999).
S. R. Arridge, “Optical tomography in medical imaging,”
Inverse Probl. 15, R41–R93 (1999).
D. A. Boas, D. H. Brooks, E. L. Miller, C. A. DiMarzio, M.
Kilmer, R. J. Gaudette, and Q. Zhang, “Imaging the bodywith diffuse optical tomography,” IEEE Signal Process.
Mag. 18, 57–75 (2001).
A. P. Gibson, J. C. Hebden, and S. R. Arridge, “Recent
advances in diffuse optical imaging,” Phys. Med. Biol. 50,
R1–R43 (2005).
G. Holte, “On a method of calculating the density of

17.
18.
19.

20.
21.

22.
23.
24.
25.

neutrons emitted from a point source in an infinite
medium,” Ark. Mat., Astron. Fys. 35A, 1–9 (1948).
K. M. Case and P. F. Zweifel, Linear Transport Theory(Addison-Wesley, 1967).
I. Kuscer and N. J. McCormick, “Some analytical results
for the radiative transfer in thick atmospheres,” Transp.
Theory Stat. Phys. 20, 351–381 (1991).
R. Aronson and N. Corngold, “Photon diffusion coefficientin an absorbing medium,” J. Opt. Soc. Am. A 16, 1066–1071(1999).
R. Elaloufi, R. Carminanti, and J.-J. Greffet, “Definition of
the diffusion coefficient in scattering and absorbing media,”
J. Opt. Soc. Am. A 20, 678–685 (2003).
R. Pierrat, J.-J. Greffet, and R. Carminati, “Photon
diffusion coefficient in scattering and absorbing media,” J.
Opt. Soc. Am. A 23, 1106–1110 (2006).
I. Freund, “Surface reflections and boundary conditions for
diffusive photon transport,” Phys. Rev. A 45, 8854–8858(1992).
R. C. Haskell, L. O. Svaasand, T. T. Tsay, T. C. Feng, and
M. S. McAdams, “Boundary-conditions for the diffusion
equation in radiative transfer,” J. Opt. Soc. Am. A 11,
2727–2741 (1994).
R. Aronson, “Boundary conditions for diffuse light,” J. Opt.
Soc. Am. A 12, 2532–2539 (1995).
S. Fantini, M. A. Franceschini, and E. Gratton, “Effective
source term in the diffusion equation for photon transport
in turbid media,” Appl. Opt. 36, 156–163 (1997).
A. F. Kostko and V. A. Pavlov, “Location of the diffusingphoton source in a strongly scattering medium,” Appl. Opt.
36, 7577–7582 (1997).
X. Intes, B. Le Jeune, F. Pellen, Y. Guern, J. Cariou, and J.
Lotrian, “Localization of the virtual point source used in
the diffusion approximation to model a collimated beamsource,” Waves Random Complex Media 9, 489–499 (1999).
T. Spott and L. O. Svaasand, “Collimated light sources inthe diffusion approximation,” Appl. Opt. 39, 6453–6465(2000).
A. D. Kim and A. Ishimaru, “Optical diffusion of continuouwave, pulsed, and density waves in scattering media andcomparisons with radiative transfer,” Appl. Opt. 37,
5313–5319 (1998).
L. G. Henyey and J. L. Greenstein, “Diffuse radiation in thegalaxy,” Astrophys. J. 93, 70–83 (1941).
V. A. Markel, “Modified spherical harmonics method forsolving the radiative transport equation,” Waves Random
Complex Media 14, L13–L19 (2004).
G. Panasyuk, J. C. Schotland, and V. A. Markel, “Radiativetransport equation in rotated reference frames,” J. Phys. A
39, 115–137 (2006).
V. A. Markel and J. C. Schotland, “Symmetries, inversion
formulas,andimage
reconstructionfor
opticaltomography,” Phys. Rev. E 70, 056616 (2004).