Journal of Physics Communications

PAPER • OPEN ACCESS

Experimental tests of indicators for the degree of validness of the
diffusion approximation
To cite this article: Antoine Capart et al 2021 J. Phys. Commun. 5 025012

View the article online for updates and enhancements.

This content was downloaded from IP address 202.253.44.186 on 02/03/2021 at 06:10

J. Phys. Commun. 5 (2021) 025012

https://doi.org/10.1088/2399-6528/abe4e1

PAPER

OPEN ACCESS

Experimental tests of indicators for the degree of validness of the
diffusion approximation

RECEIVED

4 December 2020
REVISED

5 February 2021
ACCEPTED FOR PUBLICATION

10 February 2021
PUBLISHED

26 February 2021

Antoine Capart1,2,4, Shunsuke Ikegaya1,4, Eiji Okada3, Manabu Machida1

and Yoko Hoshi1

1Hamamatsu University School of Medicine, Institute for Medical Photonics Research, Hamamatsu, 431-3192, Japan
Aix Marseille University, CNRS, Centrale Marseille, Institut Fresnel, Marseille, 13013, France
3
Keio University, Department of Electronics and Electrical Engineering, Yokohama, 223-8522, Japan
4
These authors contributed equally.
2

E-mail: yhoshi@hama-med.ac.jp
Keywords: near-infrared spectroscopy, diffusion approximation, reduced scattering coefﬁcient, diffusive regime

Original content from thiswork may be used under
the terms of the Creative
Commons Attribution 4.0
licence.

Abstract
Any further distribution of The diffusion approximation has been one of the central topics in near-infrared spectroscopy (NIRS).
this work must maintain
When NIRS measurements are analyzed by the diffusion theory, the measurements must be
attribution to theauthor(s) and the title ofperformed in the diffusive regime. However, since most of past researches have focused on theoreticalthe work, journal citationor qualitative nature of the diffusion approximation, it is not easy to know if each measurement isand DOI.
designed in the diffusive regime. In this paper, we consider the diffusion approximation quantitatively
and propose indicators that quantify the degree of validness of the diffusion approximation. The
difference between the measurement and diffusion theory can be evaluated with the χ2 value, ℓ1 andℓ2 norms, and Kullback-Leibler divergence. We conduct a liquid phantom experiment to test the
proposed χ2 value. Moreover, the χ2 value is further investigated by Monte Carlo simulations. We
ﬁnd the χ2 value becomes signiﬁcantly large when measurements are performed in the nondiffusive ortransport regime. The proposed indicators similarly work. In particular, the χ2 value is shown to work
as an indicator which evaluates the degree of validness of the diffusion approximation. These
indicators are general and can be used for different numerical, experimental, and clinicalmeasurements in NIRS.

1. Introduction
Since near-infrared light is absorbed and scattered in biological tissue, the light detected on the surface has theinformation on the medium in which light propagates. Hence we can access optical properties of the medium
such as the absorption and scattering coefﬁcients by illuminating the medium by near-infrared light. Nearinfrared spectroscopy (NIRS) including optical tomography often relies on the diffusion approximation to theradiative transport equation. When optical properties of biological tissue are retrieved using the diffusion
theory, measurements must be conducted in the diffusive regime, in which the source and detector are well
separated and the absorption coefﬁcient μa is sufﬁciently smaller than the reduced scattering coefﬁcient m¢s .
Although it has been one of the central topics in NIRS to identify the diffusive regime, most researches havequalitatively studied the diffusion approximation. Hence, in this paper, we revisit the diffusion approximation
quantitatively and provide an indicator which helps design measurements in the diffusive regime.
The transition from the transport regime to the diffusive regime takes place asymptotically when the
propagation distance of light is large and ma  m¢s . That is, there is no clear boundary between the two regimes.
However, there is a practical need of drawing the border. In NIRS, optical properties of biological tissue aredetermined by the comparison between measurements and simulational results. Since the computational
cost of the radiative transport equation is quite often too expensive, the diffusion equation should be used if
measurements are performed in the diffusive regime. Past researches on the transition from the transport regime
to diffusive regime include the following studies.

© 2021 The Author(s). Published by IOP Publishing Ltd

J. Phys. Commun. 5 (2021) 025012

A Capart et al

In [1], the ratio of the ﬂuence rate and its value from the diffusion approximation was plotted against
the distance between the source and detector for different μa values. The plot shows that the diffusion
approximation works well when ma m¢s = 1 100 and the distance is larger than 10 m¢s and the approximation
becomes worse when the absorption exceeds ma m¢s = 1 25. Although such tendencies can be read from suchqualitative study, it is not easy to capture how satisfactory the diffusion approximation holds for given μa andsource-detector distance. Here we will present an indicator which shows the degree of validness of the diffusion
approximation.
In [2], time-resolved measurements were performed in the transmission geometry and measured temporal
proﬁles were compared with the prediction by the diffusion theory. They concluded that the diffusion
approximation breaks when the width of the slab becomes smaller than 10ℓ* (ℓ* is the transport mean freepath). In [3], solutions to the time-independent diffusion equation were compared with Monte Carlo solutionsin the reﬂection geometry for a semi-inﬁnite medium. From the comparison, the criterion ma < m¢s 4 can beread as a condition that the diffusion approximation holds. Moreover they found that the diffusion results are
more accurate for small g. By investigating the Brownian motion using a low coherence interferometry, it was
experimentally shown that the transition to the diffusive regime depends on the anisotropic factor g even whenm¢s is unchanged [4]. The breakdown threshold 8ℓ* was suggested in [5]. In their study, this breakdown threshold
was obtained through time-resolved experiments in the transmission geometry with a sample of various
thicknesses. By observing temporal proﬁles of transmitted light, slabs of width greater than 8 m¢s were deﬁned to
be in the diffusive regime [6, 7]. In [8], the condition ma < m¢s 10 was proposed. The limits of the diffusion
approximation was considered with a two-dimensional numerical phantom of the human brain and a threedimensional cubic numerical phantom by the comparison between the diffusion equation and radiativetransport equation calculated by ﬁnite-difference scheme [1, 9]. Their calculation shows that the ﬂuence rate by
the diffusion approximation becomes half of its true value when ma = m¢s 2 and the distance is 8 in the unit of
1 m¢s . It was reported that the deviation of the ﬂuence rate increases when the anisotropic factor g (introduced
below) is close to 1, which is typical in biological tissue [10]. The deviation of the dynamic intensity-intensity
correlation function from the prediction of the diffusion theory was experimentally observed [11].
The extrapolated boundary condition has often been used to further approximate the diffusion equation.
It is known that those approximate solutions are quite precise in the diffusive regime. The reﬂectance andtransmittance were derived as functions of time from the diffusion equation solved using the extrapolated
boundary condition [12]. The prediction from those concise formulae agreed well with Monte Carlosimulations. In [13], the validity of the diffusion model was conﬁrmed by experiments and Monte Carlosimulations. In [14], the validity of the diffusion model was studied in frequency domain using the spherical andcylindrical geometries in addition to the slab geometry. The effect of different approximations to the diffusion
equation was studied [15, 16]. The diffusion approximation with the extrapolated boundary condition wastested using Monte Carlo in the slab geometry [17].
It was found that the quality of the diffusion approximation depends on the refractive index [13, 18]. Resultsfor different index mismatches were reported [19]. Moreover, the validity of the diffusion approximation on thesurface was studied [20].
The nondiffusive regime was studied by using other approximations to the radiative transport equation.
When ma > m¢s , a clear difference was observed between solutions to the diffusion equation and telegrapher
equation, which takes the ballistic component of photon propagation at short times into account [21]. The
diffusion approximation was compared to the P1 and P3 approximations [22].
The light propagation is said to be in the diffusive regime if the solution to the radiative transport equationand the solution to the diffusion equation match. However, there is no unique way of setting the border of the
diffusive regime because the change from the transport regime to diffusive regime takes place gradually. To draw
a line between two regimes, a criterion must be set. In this paper, a constant a (see below) is used as the threshold.
The solution to the radiative transport equation can be obtained from measurements since infrared light inbiological tissue is governed by the radiative transport equation. In this case, the threshold a must be set in such a
way that the determined boundary of the diffusive regime is not affected by measurement noise or error.
In this paper, we propose to use an indicator such as the χ2 value to identify the diffusive regime. Thus, thedegree of validness of diffusion approximation is given as a χ2 value. Although the change from the nondiffusive
or transport regime to the diffusive regime takes place asymptotically, we ﬁnd that χ2 values rather rapidlychange when measurements shift from one regime to the other. The behavior of the χ2 value is tested both byliquid phantom experiments and Monte Carlo simulations. By Monte Carlo simulation, it is conﬁrmed that
other indicators ℓ1 and ℓ2 norms, and the Kullback-Leibler (KL) divergence work as well.

2

J. Phys. Commun. 5 (2021) 025012

A Capart et al

2. Methods
2.1. Diffusion theory
Here, the diffusion theory is developed for a three-dimensional semi-inﬁnite medium. In addition to thereﬂection on the boundary, we take into account the instrument response function, diameters of the opticalﬁbers, and numerical aperture in the diffusion approximation. The approximate form ΦDA(t) of the detected
light Φ(ρd, t) (deﬁned below) will be derived.
2.1.1. The radiative transport equation
Let us consider the three-dimensional half space W = {r Î 3; -¥ < x < ¥, -¥ < y < ¥, 0 < z < ¥}.
We consider the speciﬁc intensity –––– I (r, sˆ , t ) of near-infrared light at position r = (x, y, z) in direction sˆ = (sin J
cos j, sin J sin j, cos J) at time t. Here, the polar angle ϑ and azimuthal angle j vary in 0  ϑ  π and 0  j < 2π,
respectively. When absorption coefﬁcient μa and scattering coefﬁcient μs are positive constants, the speciﬁc intensity
I (r, sˆ , t ) obeys the following radiative transport equation.
⎛1 ¶⎞⎜+ sˆ ·  + ma + ms ⎟ I (r , sˆ , t )⎝ c ¶t⎠= ms ò

2p0p

ˆ t ),
ò0 p (ˆs , sˆ¢) I (r, sˆ¢ , t ) sin J¢ dJ¢dj¢ + q (r, s,

where c > 0 is the speed of light in the medium and p (ˆs, sˆ¢) is the scattering phase function, which is normalized
as
2p

p

ò0 ò0 p (ˆs, sˆ¢) sin J dJdj = 1(2)

for all J¢ Î [0, p ], j¢ Î [0, 2p ). We assume p (ˆs, sˆ¢) = p (- sˆ¢ , - sˆ). The anisotropic factor g is introduced as
g=

2p

p

ò0 ò0 (ˆs · sˆ¢) p (ˆs, sˆ¢) sin J dJdj.

(3)ˆ t ) on the rightThe reduced scattering coefﬁcient is given by m¢s = (1 - g ) ms . There is the source term q (r, s,
hand side of (1).
In addition to the initial condition I (r, sˆ , 0) = 0, the boundary condition is imposed asˆ t ) = R (ˆ)
I (r, s,s I (r, sˆr , t )

on

z=0

(4)

for 0  ϑ < π/2. Here, sˆr = (sin J cos j, sin J sin j, -cos J). Let n denote the refractive index of themedium. We assume air outside the medium. Assuming unpolarized light, the Fresnel reﬂection R (ˆ)
s is given by
R (ˆ)
s =

2
⎡
⎛ nm - cos J ⎞2 ⎤1 ⎢⎛ n cos J - mn ⎞⎜⎟ +⎜ n⎟ ⎥2 ⎢⎣⎝ n cos J + mn ⎠⎝ nmn + cos J ⎠ ⎥⎦(5)

for mc  cos J  1 and R (ˆ)
s = 1 when 0  cos J < mc . Here, mn = 1 - n2 (1 - cos2 J) and2
mc = n - 1 n.
Let rs = (ρs, 0+), where ρs = (xs, ys), be the position of the optical ﬁber for the incident beam on the surface.
ˆ t ) asWe set ρs = 0. We write the source term q (r, s,ˆ t ) = d (r - rs) d (ˆs - zˆ) h (t )q (r, s,= d (x ) d ( y ) d (z - 0+) d (cos J - 1) d (j) h (t ) ,(6)

where δ( · ) is Dirac’s delta function, zˆ = (0, 0, 1), and h(t) is the temporal proﬁle of the incident beam. In our
measurement in the reﬂection geometry, the detector is placed on the boundary at rd = (ρd, 0) with ρd = (xd, yd).
Let  (t ) denote the response function. Considering the diameter  (r, rd) and numerical aperture  (ˆ)
s of theoptical ﬁber for the detector at rd, the measured quantity Φ(ρd, t) is expressed as
t

2p

0

0

F (rd , t ) = ò  (t - t ¢) ò
´ò

¥

ò

¥

-¥ -¥

p

s sin Jòp2 (ˆ)ˆ t ¢) dxdydJdjdt ¢ , (r , rd) I (r , 0, s,

where ρd = |ρd|. If the detector radius ρﬁber is small, we can put  (r, rd) = pr2fiber d (r - rd ). If the numericalaperture is given as p - Jmax  J  p , 0  j < 2π for some Jmax , then we can write  (ˆ)
s = cNA (ˆ)
s , whereχNA = 1 if J Î [p - Jmax , p ] and χNA = 0 otherwise.
3(7)

J. Phys. Commun. 5 (2021) 025012

A Capart et al

The function  (t ) below corresponds to the instrument response function.
tò0 (t - s) h (s) ds.

 (t ) =(8)ˆ t ) = G RTE (r, s;ˆ rs , z;ˆ t ) becomes the Green’s function for the radiative transport
If h(t) = δ(t), then I (r, s,
equation in the half space. Let us deﬁne (rd , t ) = ò

2p

p

¥

¥

s sin J ò ò (r , rd)òp2 (ˆ)-¥ -¥0

ˆ rs, z;ˆ t ) dxdydJdj.
´ G RTE (r , 0, s;(9)

We note that
t

t¢

0
t

0

F (rd , t ) = ò  (t - t ¢) ò=ò0

 (rd , t ¢ - t ) h (t ) dt dt ¢s(ò (s - t ¢) h (t ¢) dt ¢) (r , t - s) ds.
d

0

(10)

Therefore we obtaint

F (rd , t ) = ò  (rd , t - s )  (s ) ds
0

=ò

2p

0

p

¥

¥

s sin J ò ò (r , rd)òp2 (ˆ)-¥ -¥tˆ rs, z;ˆ t - s )  (s ) dsdxdydJdj.
´ ò G RTE (r , 0, s;(11)

rd = 20 mm.

(12)

0

Throughout the paper we set
The above Φ(ρd, t) corresponds to the measurement by the TRS-20, which we denote by ΦTRS(ρd, t).
2.1.2. Diffusion approximation
Let us calculate Φ(ρd, t) using the diffusion approximation. In the spirit of the diffusion approximation, we write
the speciﬁc intensity as
1
3
u (r, t ) +sˆ · J(r, t ).
4p
4pˆ t) =I (r, s,(13)

By substituting the above expression (13) in (1) after some calculations [23], we obtain
1 ¶u (r, t ) - DDu (r, t ) + ma u (r, t ) = S (r, t ) ,
c ¶t

(14)

where the diffusion coefﬁcient is given by [24]
D=

1
,
3m¢s

(15)

where m¢s = (1 - g ) ms . The dependence of the diffusion coefﬁcient on parameters in the radiative transport
equation has been studied. Indeed, D = 1 [3 (m¢s + ma )] is derived from the straightforward calculation.
However, the form 1 (3m¢s ) was suggested by more detailed studies [25–27, 24]. In Appendix A, we develop analternative approach to the diffusion approximation. The diffusion coefﬁcient (15) is derived in Appendix A.
The source term is obtained as
S (r, t ) = Q 0 (r, t ) -1
 · Q1(r, t ).
m¢s

(16)

Here we introduced
Q 0 (r, t ) = ò2ppˆ t ) sin J dJdj ,ò q (r, s,

0
0
2p p

Q1(r, t ) = ò0ˆ t ) sin J dJdj.
ò0 sˆq (r, s,(17)

The current J(r, t) is obtained as
J(r, t ) = - D u (r, t ) +1
Q1(r, t ).
m¢s

(18)

Let us consider the boundary condition for the diffusion equation (14). With the diffuse surface reﬂection,
we obtain the following Robin boundary condition.
4

J. Phys. Commun. 5 (2021) 025012

A Capart et al¶u(r, t ) = 0 at z = 0,¶zu (r, t ) - 2DA(19)

where [28, 29]
A=

1 + rd
,1 - rd(20)

whererd = - 1.4399n-2 + 0.7099n-1 + 0.6681 + 0.063n.

(21)

Moreover the initial condition is given by u(r, 0) = 0.
Let us calculate  DA (rd , t ) by making the diffusion approximation to  (rd , t ). Then we have within the
diffusion approximation,
F (rd , t ) »tò0 DA(rd , t - t ) (t ) dt.

(22)

Let us begin by considering the diffusion approximation to the Green’s function GRTE. In the radiativetransport equation (1), we set
ˆ t ) = d (r - rs) d (ˆs - zˆ) d (t ).
q (r, s,(23)

Then we have Q0(r, t) = δ(r − rs)δ(t) and Q1(r, t ) = zˆd (r - rs) d (t ), and the source term for the diffusion
equation is obtained as⎛⎞
1 ¶S (r, t ) = d (x - x s) d ( y - ys) d (t ) ⎜⎜d (z - 0+) d (z - 0+) ⎟⎟.
m¢s ¶z⎝⎠(24)

From the fact that rd ≠ rs, we see J(rd, t) = − D∇u(rd, t). Thus we obtain DA (rd , t ) = ò

2p

0

p

¥

¥

⎛ 1

⎞

3

s sin J ò òJ(rd , t ) · sˆ⎟ dxdydJdj
 (r , rd) ⎜ u (rd , t ) +òp2 (ˆ)
-¥ -¥
⎝ 4p
⎠
4p

⎛ 1 - cos Jmax⎞
3
(1 - cos 2Jmax ) ⎟ u (rd , t ) ,= pr 2fiber ⎜
+
⎝
⎠
2
16A

(25)

where we assumed that ρﬁber is small and  (ˆ)
s = cNA (ˆ)
s was used. Using the Green’s function G (r, r¢ ; t ) for the
diffusion equation with the source term d (r - r¢) d (t ), we obtaintu (rd , t ) = ò ò G (rd , r; t - t ¢) S (r, t ¢) d rdt ¢0W= G (rd , x s, ys , 0+) +¶1 ¥d (z - 0+) G (x s, ys , z , rd ; t ) dz ,
ò
0
¢
ms
¶z

(26)

where we used integration by parts and the reciprocal property G (r, r¢ ; t ) = G (r¢ , r; t ). Hence by making useof the Robin boundary condition,⎛
3 ⎞⎟u (rd , t ) = ⎜1 +G (rd , rs; t ).
⎝
2A ⎠

(27)

Therefore we arrive at⎛ 1 - cos Jmax⎞
3
 DA (rd , t ) = pr 2fiber ⎜+(1 - cos 2Jmax ) ⎟
⎝
⎠
2
16A
⎛
3 ⎞⎟
´ ⎜1 +G (rd , rs; t ).
⎝
2A ⎠

(28)

We can obtain the Green’s function as
G (r, r¢ ; t ) =(x - x ¢ ) 2 + ( y - y ¢ ) 2
ce -ma ct
4Dct
eg (z , z ¢ ; t ).
32
(4pDct )

(29)

Here,
⎛ z + z ¢ + 2Dct ℓ ⎞z + z ¢ Dct
2
pDct e ℓ + ℓ2 erfc ⎜
⎟
⎝
⎠
ℓ
4Dct
2
2
2
ℓ
z
+
z
¢+
Dct
⎛
(z - z ¢ ) 2(z + z ¢ ) 2
(
)
z
+
z
¢
z + z ¢ + 2Dct ℓ ⎞2pDct e- 4Dct e ( 4Dct ) erfc ⎜= e - 4Dct + e- 4Dct ⎟,
⎝
⎠
ℓ
4Dct
(z - z ¢ ) 2(z + z ¢ ) 2

g (z , z ¢ ; t ) = e - 4Dct + e - 4Dct -2¥

where ℓ = 2DA and the complementary error function is given by erfc (z ) = p ò e-t dt . We note thatz
G = 0 when t < 0. The expression (29) is obtained by straightforward calculation [30].
5

2

(30)

J. Phys. Commun. 5 (2021) 025012

A Capart et al

We can write
FTRS (t ) = C1 ¢¥¥t

ò0 DA(rd , t - t ) (t ) dt,

(31)

t

where C1 ¢ = ò FTRS (t ) dt ò ò  DA (rd , t - t )  (t ) dtdt . When ΦTRS is compared with results from the000
diffusion theory, precise knowledge of rd , Jmax , A is not required because the constant factors on the right-handside of (28) can be absorbed by the constant C1 ¢. That is, the following relation holds in the diffusive regime.
FTRS (t ) = C1 FDA (t ) ,t

FDA (t ) = ò G (rd , rs; t - t )  (t ) dt
0(32)

for some positive constant C1.
2.2. Liquid phantom
2.2.1. Time-resolved spectroscopy instrument
A two-channel time-resolved spectroscopy instrument (TRS-20, Hamamatsu Photonics, Hamamatsu, Japan)
was used to measure optical properties of the liquid phantom. The TRS-20 system consists of three pulsed laser
diodes of wavelengths 759 nm, 794 nm, and 834 nm each of which has the duration of around 100 ps at therepetition rate of 5 MHz, a high-speed photomultiplier tube for single photon detection, and the TCSPC system,
which has a circuit for time-resolved measurements based on the time correlated single photon counting(TCSPC) method. The TCSPC system can detect arrival photons with approximately 10 ps time resolution.
Although temporal proﬁles of emerging photons from the sample can be simultaneously accumulated up to a
count rate of 500 kcps for all three wavelengths, temporal proﬁles are measured at less than 150 kcps to avoid
pile-up distortion. The minimum data acquisition time is 100 ms. The instrumental response of the TRS-20 ismeasured by placing the incident optical ﬁber (200 μm diameter, NA = 0.25) opposite the detecting optical ﬁber(bundle ﬁber, 3 mm diameter, NA = 0.29) with a neutral density ﬁlter between them. The FWHM (full widthhalf maximum) of the instrument response function of the TRS-20 is about 350 ps at each wavelength. In thisstudy, we used the wavelength of 794 nm.
2.2.2. Liquid phantom measurements
Figure 1 shows the experimental setup. A liquid phantom was prepared in a bucket of radius 8–10 cm and height
20 cm. The 10% Intralipos (Otsuka, Tokyo, Japan) and green brown ink (Chugai Kasei, Tokyo, Japan) were
gradually added to water (initially 3 L) to control the scattering and absorption, respectively. The incident anddetecting optical ﬁbers from the TRS-20 were attached to a holder. The separation of the ﬁbers was 2 cm, i.e.,ρd = 20 mm. We set the holder on the surface of the liquid phantom and supported optical ﬁbers by using two
arm-type clamps. The liquid phantom was stirred with a magnetic stirrer (SW-030, Nissin Science, Tokyo,
Japan) except for during measurement. We conﬁrmed that the bucket is deep enough and the liquid phantom
can be regarded as a semi-inﬁnite medium even for smallest ma , m¢s used in the experiments.
In the diffusive regime for large m¢s , optical properties of the liquid phantom were estimated by the TRS-20,
which obtains μa and m¢s using the Levenberg-Marquardt method by minimizing the difference between themeasured temporal proﬁle and the time-resolved reﬂectance calculated from the diffusion theory (see [31]for details). We conﬁrmed that estimated μa values by the TRS-20 were consistent with values fromspectrophotometric measurements (UV-3100, Shimadzu, Kyoto, Japan). The absorption coefﬁcient μa of theliquid phantom can be kept the same for different m¢s by adjusting the ink so that the ratio of the ink to the total
amount of the liquid phantom is ﬁxed. Making use of the fact that the reduced scattering coefﬁcient m¢s of theliquid phantom linearly depends on the concentration of Intralipos, we determined m¢s values in the nondiffusive
regime by extrapolating the line ﬁtted by linear regression taking several m¢s values in the diffusive regime in theplane of the Intralipos concentration and m¢s . The ma , m¢s values used for the experiments are summarized intables B1 and B2 in Appendix B.
2.3. Monte Carlo simulations
In addition to the liquid phantom experiment, we employ Monte Carlo simulation, in which optical propertiesof the medium can be precisely speciﬁed. The validity of Monte Carlo simulations for a liquid phantom was
conﬁrmed in [32]. Our Monte Carlo simulation was implemented using the variance reduction technique[33, 34]. While photons propagate in the medium, they are scattered according to the scattering phase function
at each scattering point. The Henyey-Greenstein model [35] was used. The anisotropy factor was set to g = 0.9
for μa = 0.06 mm−1 and μa = 0.1 mm−1, and g = 0.33 otherwise. The refractive index was set to n=1.33. On theboundary at z = 0, the reﬂection and refraction caused by the refractive index mismatch were considered.
6

J. Phys. Commun. 5 (2021) 025012

A Capart et al

Figure 1. Schematic ﬁgure of the experimental setup.

The diameter and numerical aperture of the incident ﬁber at the origin are 0.2 mm and 0.250, respectively.
The detecting optical ﬁber was placed at a distance of 20 mm from the source (that is, ρd = 20 mm). For thedetector, the diameter and numerical aperture were set to 3 mm and 0.260, respectively.
The weight of detected photon packets was accumulated into time-resolved arrays for absorption
coefﬁcients. The duration of an array element was 10 ps, corresponding to the temporal resolution of the TRS
instrument used in our experiment. In each Monte Carlo simulation, 108 photons were launched.
In Monte Carlo simulations,  MC(rd , t ) corresponding to  (rd , t ) is numerically computed. We have
 (r, rd) = crd (r ) and  (ˆ)
s = cNA (ˆ)
s , where crd (r ) = 1 if |ρ − ρd|  ρﬁber and crd (r ) = 0 otherwise. We set
Jmax = 0.26, ρﬁber = 1.5 mm. We introduce
FMC (t ) =tò0  MC(rd , t - t ) (t ) dt.

(33)

We haveFTRS (t ) = C 2 FMC (t ).

(34)

The constant C2 is obtained as¥

C2 =ò0 FTRS (t ) dt¥

ò0 FMC (t ) dt

.

(35)

2.4. The χ2 indicator
Let us normalize the time-resolved functions as
 EXP (t ) =F DA (t ) =F

FEXP (t )¥ò0 FEXP (t ) dt
FDA (t )¥ò0 FDA (t ) dt

,

,

(36)

(37)

where EXP means TRS or MC. The integrals in (36) and (37) are computed with the trapezoidal rule.
Measurement times are discrete values given by
tk = kDt ,Dt = 10 ps,

Here, M = 1024 for TRS and M = 800 for MC.
7

k = 0, ¼, M - 1.

(38)

J. Phys. Commun. 5 (2021) 025012

A Capart et al

There are different ways of measuring the difference between two functions. To evaluate the difference
between the measured time-resolved reﬂectance, which is the solution to the radiative transport equation, andthe temporal proﬁle by the diffusion approximation, we introduce the χ2 value askb

c2 = å

k=ka


∣F

EXP

 (tk)∣2( tk ) - F.
DA ( tk )
F
DA

(39)

 DA (tk ) for k = 0, K, ka − 1 and k = kb + 1, K, M − 1
The lower and upper limits ka, kb are chosen such that FDA
are less than one-tenth of the peak height of F (tk ). The integers ka, kb will be chosen such that values which are
more than 10% of the peak height are considered in the sum in (39).
In addition to the χ2 value, we propose the ℓ1 norm, ℓ2 norm, and KL divergence. They are given by
 EXP (tk ) andå kkb= ka ∣ fk - gk∣, å kkb= ka ∣ fk - gk∣2 , and å kkb= ka fk ln ( fk gk ), respectively. Here, fk = F (tk ). These indicators behave similarly if scaled properly (see ﬁgure 8 below). Hence we will focus on
gk = Fthe χ2 indicator in this paper.
A few remarks are necessary: (i) We note that the χ2 indicator introduced above has nothing to do with the2χ test in statistics. The χ2 value in (39) is used to measure the distance between the solution to the radiativetransport equation and the solution to the diffusion equation. In this paper, we admit that experimentallyobtained data obey the radiative transport equation and treat FEXP (t ) as the solution to the radiative transport
equation. This means that we do not argue how precisely measured experimental data are described by theradiative transport equation. (ii) The KL divergence sometimes becomes negative whereas the χ2 value, ℓ1 norm,
and ℓ2 norm are always nonnegative.
We separate obtained data for different μa and m¢s into the diffusive and nondiffusive regimes using the χ2value (39). To this end, we set a threshold a. The measurement is identiﬁed as in the diffusive regime if χ2 < a.
Otherwise the measurement is considered to be in the nondiffusive or transport regime.
DA

3. Results
Figure 2 shows χ2 values (39) obtained from the liquid phantom experiment for the absorption coefﬁcient
μa = 0.0023 mm−1 (purple), 0.012 mm−1 (green), and 0.021 mm−1 (red). It is found that the χ2 value becomessigniﬁcantly large for μa = 0.012 mm−1 and 0.021 mm−1 when the source-detector distance rd m¢s in the unit of
1 m¢s is less than 5. In the experiment ρd = 20 mm. Tables B1 and B2 in Appendix B summarize χ2 values
obtained in the liquid phantom experiment. Figure 2 illustrates that the indicator takes large values when μa islarge.
The diffusive and nondiffusive regimes are segmented in ﬁgure 3. The threshold a is set to 5. Since the χ2value rapidly grows as shown in ﬁgure 2 once performed measurements move from the diffusive regime to thetransport regime, the same result is obtained if a = 4 or a = 6. Figure 4 shows time-resolved reﬂectances for TRS (t ) from the TRS-20 for the
μa = 0.021 mm−1 (the top curve in ﬁgure 3). In ﬁgure 4, temporal proﬁles F
smallest m¢s (the leftmost point in ﬁgure 3) and largest m¢s (the rightmost point in ﬁgure 3) are compared with the DA (t ) from the diffusion equation. There is a noticeable discrepancy in the left panel
time-resolved reﬂectances Fof ﬁgure 4 for the smallest m¢s = 0.207 mm-1.
To explore the intermediate region between the transport and diffusive regimes more carefully, we employ
Monte Carlo simulations. Obtained χ2 values from Monte Carlo simulations are summarized in tables B3, B4,
and B5 in Appendix B.
In ﬁgure 5, we set a = 10. Then measurements belong to the nondiffusive regime when rd m¢s and ma m¢s are
small.
Let us try other values of the threshold a and see how the deﬁned diffusive regime is affected by a. In the left
panel of ﬁgure 6, we set a = 20. Since the criterion a = 20 is loose, the deﬁned diffusive regime is wider than theregion deﬁned in ﬁgure 5. That is, there are more blue open squares in the left panel of ﬁgure 6. In the right panelof ﬁgure 6, the a value is set to 5. This a = 5 is severe and measurements are classiﬁed in the nondiffusive regime
unless ρd and μa are very small. When measurements in the non-diffusive regime are of interest, a should betaken rather large (e.g., a = 20), so that the points marked in red are surely in the non-difﬁsive regime. On theother hand, if measurements in the diffusive regime are of interest, a should be set to small values (e.g., a = 5),
so that blue open squares certainly belong to the difﬁsive regime. Figure 7 shows results for Monte Carlosimulations with μa = 0.01 mm−1. All curves in ﬁgure 7 are plotted together in ﬁgure 8 after scaling. The ℓ1
norm, ℓ2 norm, and KL divergence are scaled as
8

J. Phys. Commun. 5 (2021) 025012

A Capart et al

Figure 2. Obtained χ2 values are plotted as functions of rd m¢s , where the source-detector distance ρd = 20 mm. Results for three μavalues 0.0023 mm−1 (purple), 0.012 mm−1 (green), and 0.021 mm−1 (red) are shown.

Figure 3. Results from the liquid phantom experiment are shown. Red solid circles are used for the nondiffusive regime andmeasurements in the diffusive regime are marked by open blue circles. Black solid lines show (ma m¢s ) = rd ma (rd m¢s ) for ρd = 20 mm
and μa, from the top, 0.021/mm, 0.012/mm, and 0.0023/mm.

Figure 4. Comparison between time-resolved curves from the TRS-20 and results calculated with the diffusion approximation for
μa = 0.021 mm−1 when (Left) m¢s = 0.207 mm-1 and (Right) m¢s = 1.46 mm-1.

y c2 (x = 4)y (x ) - y (x = 28) + y c2 (x = 28)
y (x = 4) - y (x = 28) + y c2 (x = 28)

,

(40)

where x = rd m¢s and y c2 is the χ2 value. Here, y is either y ℓ1, y ℓ2 , or yKL, which are deﬁned in the same way as
y c2 . Table 1 shows derivatives dy/dx (central differences) for the indicators χ2, ℓ1, ℓ2, and KL divergence, which
are calculated from numerical values in ﬁgure 8. It is observed that derivatives take negative large values atrd m¢s = 6. For large rd m¢s , derivatives are small positive numbers. When rd m¢s is greater than 20, the solution of
the radiative transport equation and the solution of the diffusion equation are almost indistinguishable. Since
9J. Phys. Commun. 5 (2021) 025012

A Capart et al

Figure 5. Same as ﬁgure 3 but data were taken from Monte Carlo simulation. Measurements in the nondiffusive regime are marked by
red solid squares and measurements in the diffusive regime are denoted by blue open squares. Black solid lines show (ma m¢s ) = rd ma(rd m¢s ) for ρd = 20 mm and μa, from the top, 0.1/mm, 0.06/mm, 0.02/mm, 0.01/mm, and 0.002/mm.

Figure 6. Same as ﬁgure 5 but (Left) a = 20 and (Right) a = 5 are used.

Figure 7. From the left, values of χ2, ℓ1, ℓ2, and KL divergence are shown in the case of Monte Carlo simulations with
μa = 0.01 mm−1.

the number of data kb − ka + 1 in (39) increases when rd m¢s is large and the time-resolved curve is broadened,derivatives take small positive values for large rd m¢s . Indeed in the case of the Monte Carlo simulation with
μa = 0.01 mm−1, (m¢s mm-1; k a, kb ) = (0.2; 4, 54), (0.4; 8, 76), (0.6; 11, 94), (0.8; 14, 109), (1.0; 17, 122),(1.2; 19, 134), and (1.4; 22, 145), respectively.

4. Discussion
Let us compare  MC and  DA . We write MC = C3  DA.

(41)

We found A = 2.785 and C3 = 1.3. The result that C3 is slightly greater than 1 attributes to the fact that the
diffusion approximation breaks on the boundary. In section 2, the source was placed on the surface when ΦDA(t)was calculated. In the diffusion approximation, quite often the source S(r, t) is placed inside the medium (z > 0)
10

J. Phys. Commun. 5 (2021) 025012

A Capart et al

Figure 8. Values of χ2, ℓ1, ℓ2, and KL divergence are shown in the case of Monte Carlo simulations with μa = 0.01 mm−1. To
compare, the curves are scaled according to (40).

Table 1. Derivatives (central differences) of the curves
in ﬁgure 8 are shown.

rd m¢s(c 2) ¢(ℓ1)¢(ℓ 2)¢

(KL)¢

6.0
10.0
14.0
18.0
22.0
26.0

−4.68
−0.93
−0.23
−0.05
0.00
0.02

−4.32
−1.25
−0.34
−0.12
0.07
0.05

−4.16
−1.19
−0.37
−0.13
0.07
0.06

−5.07
−0.71
−0.15
−0.05
0.02
0.02

ˆ t ) represents the optical ﬁber on theabout one transport-mean-free-path away although the source term q (r, s,boundary. This has an effect of improving the value of C3. To see the effect, let us consider the solution with theextrapolated boundary condition. In this case, the Green’s function for the diffusion equation in the half space isobtained by using gEBC (z , z ¢ ; t ) instead of g (z , z ¢ ; t ) in (29). Here,(z - z ¢ ) 2

gEBC (z , z ¢ ; t ) = e - 4Dct - e -(z + z ¢+ 2ℓ ) 2
4Dct

.

(42)

Suppose the source is placed at rs = (xs, ys, d), where d is small (i.e., 2Ad = ct and 4DA2 = ct). Then we havegEBC(0, d; t) = [1 + d/(2DA)]gEBC(0, 0; t). We ﬁnd that d = 0.871 m¢s results in C3 = 1.0 (the parameters,
ρd = 20 mm, ρﬁber = 1.5 mm, Jmax = 0.26, and n = 1.33 are used). Thus, the optimal depth is about thetransport mean free path.
In addition to values shown in ﬁgure 8, we further tested the symmetric χ2 value, which is given bykb
 EXP (tk ) - F
 DA (tk )∣2 (F EXP (tk ) + F
 DA (tk )). The obtained values were close to the χ2 values.
2å k = ka ∣F
As is explained in section 2.2, values of m¢s for the liquid phantom were determined using about several
measured points for which the source-detector distance is large enough ( > 20ℓ*). However, since m¢s obtainedby the TRS-20 has only about two signiﬁcant digits, the error in estimating m¢s in the nondiffusive regime may not
be negligible. To ﬁnd clearer boundaries between the diffusive and nondiffusive regimes, parameter estimation
based on the radiative transport equation is necessary.
In NIRS, the reduced χ2 value was used to discuss the accuracy of optical parameters calculated by TRS [36].
Instead of evaluating statistical errors, the aim of this paper is to compare the solutions of the radiative transport
equation and diffusion equation. In this paper, the χ2 function was introduced solely as an indicator which takes
large values in the nondiffusive regime.
The χ2 value (39) depends on the choice of ka, kb. In this paper we took into account the values which are
more than 10% of the peak height. In [37] and [38], the range between 80% of the peak on the leading edge and
20% on the falling edge was used. In [39], the range from 80% on the leading edge and 1% on the falling edge wasused. In addition to 10%–10% cut for ka, kb, we tried 10%–50%, 50%–10%, 50%–50%, 40%–10%, 80%–20%,
and 80%–20%, and conﬁrmed that the qualitative conclusion remains the same. In order to capture thedifference between FEXP (t ) and ΦDA(t), we chose the range 10%-10%, which is rather wide. When ΦDA(t) is used11

J. Phys. Commun. 5 (2021) 025012

A Capart et al

to ﬁt FEXP (t ) for parameter identiﬁcation, however, a large value such as 80% should preferably be taken on theleading edge because the diffusion approximation easily breaks at early times.

5. Conclusion
In this paper, we have proposed indicators such as the χ2 value (39) to judge if the measurement is in the diffusive
regime or in the nondiffusive regime. By using the χ2 value, the validness of the diffusion approximation can bequantitatively studied. Once the criterion a is given, the diffusive regime is deﬁned by the χ2 value.
The usefulness of these indicators can be read from ﬁgure 2. As is shown in Appendix A, the diffusion
equation is obtained from the radiative transport equation when the small parameter ò goes to zero. In this sense,
the change from the transport regime to diffusive regime is not abrupt but takes place asymptotically. However,ﬁgure 2 shows that the change is rather sharp when the degree of validness of the diffusion approximation ismeasured by the χ2 indicator function.
As is seen in section 3, the deﬁned diffusive regime depends on the value of the threshold a. Whenmeasurements require that the diffusion approximation severely holds, a should be set to a small value. Suppose
we want to conduct measurements in the diffusive regime. In order to draw a line between the diffusive andnondiffusive regimes, a must be deﬁned operationally. First we try a relatively large a, so that many points in therd m¢s - ma m¢s plane are classiﬁed into the diffusive regime. Then we try smaller a values. The value of a cannot be
small compared with measurement errors. Thus a suitable a can be chosen.

Acknowledgments
We are grateful to Yukari Tanikawa for fruitful discussions on preparing the liquid phantom. This research wasdone while A. C. stayed in Hamamatsu. His stay was supported by the PRAME (Programme Régional d’Aide à la
Mobilité Etudiante) by the Provence-Alpes-Côte d’Azur Region in France, which is greatly appreciated. M. M.
acknowledges support from JSPS KAKENHI Grant No. 17K05572, 18K03438. Y. H. acknowledges support from
JSPS KAKENHI Grant No. 17H02081.

Data availability statement
All data that support the ﬁndings of this study are included within the article (and any supplementary ﬁles).

Appendix A. Asymptotic expansion
In section 2.1.2, we began with the assumption (13). Here we consider the diffusion approximation by theasymptotic expansion [40, 41]. The calculation below is close to [42] but small nonzero absorption is assumed inaddition to large x, t. Although basically we follow the calculation in [43] and arrive at the same scaled radiativetransport equation, we consider large x instead of large μs, so the relation to the scaling in [42] becomes more
visible. The imposed conditions of large x, t and small μa in this section are consistent with the conditions in
section 2.1.2.
Let us introduce a small parameter ò > 0 which relates variables in (1) to the corresponding slow variables as r = r*, 2t = t *,

 -2ma = m*a .

(A.1)

According to the spirit of the diffusion approximation, we consider the speciﬁc intensity away from the sourceand boundaries, and deal with the homogeneous equation. By using the slow variables the radiative transport
equation is expressed as⎛ 2 ¶⎞ ⎛ r*t* ⎞ˆ 2⎟+  sˆ · * +  2m*a + ms ⎟ I ⎜ , s,⎜
⎝ c ¶t *
⎠ ⎝
 ⎠
⎛ r*
2p p
t* ⎞
= ms ò ò p (ˆs, sˆ¢) I ⎜ , sˆ¢ , 2 ⎟ sin J¢ dJ¢dj¢ ,
0
0
⎝
 ⎠

where ∇* = ∂/∂r*. Hereafter we drop the superscript ∗.
12

(A.2)

J. Phys. Commun. 5 (2021) 025012

A Capart et al

Let us expand the speciﬁc intensity as⎛rt ⎞ˆ 2 ⎟ = I (0) (r, s,ˆ t ) +  I (1) (r, s,ˆ t)I ⎜ , s,⎝ ⎠ˆ t ) + .
+  2I (2) (r, s,(A.3)

First we collect terms of order ò0. We obtainˆ t ) = ms ò
ms I (0) (r, s,

2ppò p (ˆs, sˆ¢)00´ I (0) (r, sˆ¢ , t ) sin J¢dJ¢dj¢.

(A.4)

The above equation implies that I(0) is independent of ŝ . That is, we can writeˆ t ) = I (0) (r, t ) = u (r, t ).
I (0) (r, s,(A.5)

The terms of order ò1 yieldˆ t)sˆ · I (0) (r, t ) + ms I (1) (r, s,= ms ò

2p

0p

ò0 p (ˆs, sˆ¢) I (1) (r, sˆ¢ , t ) sin J¢dJ¢dj¢.

(A.6)

We see that I(1) is obtained as
1sˆ · I (0) (r, t )ms (1 - g )1sˆ · u (r, t ).
=ms (1 - g )ˆ t) = I (1) (r, s,(A.7)

We haveˆ t ) = - 3Dsˆ · u (r, t ) ,I (1) (r, s,(A.8)

where D is given in (15). Finally, we collect terms of order ò2.
1 ¶ (0)ˆ t ) + ma I (0) (r, t ) + ms I (2) (r, s,ˆ t)
I (r, t ) + sˆ · I (1) (r, s,c ¶t= ms ò

2p

0p

ò0 p (ˆs, sˆ¢) I (2) (r, sˆ¢ , t ) sin J¢dJ¢dj¢.

(A.9)

The left-hand side of (A.9) can be expressed as
LHS =

1 ¶
u (r, t ) - 3D (ˆs · )(ˆs · ) u (r, t )c ¶tˆ t ).
+ ma u (r, t ) + ms I (2) (r, s,(A.10)

If we integrate (A.9) over ŝ , we obtain
1 ¶u (r, t ) - DDu (r, t ) + ma u (r, t ) = 0.
c ¶t

(A.11)

Thus (14) is derived.
We note that the above equation does not have the source term because we dropped the source term in theradiative transport equation in (A.2) by focusing on large space and time.
In our derivation, the diffusion coefﬁcient D in (15), which is independent of μa, was obtained.

13

J. Phys. Commun. 5 (2021) 025012

A Capart et al

Appendix B. Values of χ2
Table B1. TRS measurements for (Left) μa = 0.0023 mm−1 and (Right) μa = 0.012 mm−1
(ρd = 20 mm). The unit of m¢s is mm−1.

m¢s

rd m¢sma m¢s

χ2

0.250

5.00

0.00921

1.95

0.318

6.36

0.00724

3.50

0.379

7.58

0.00607

3.21

0.440

8.80

0.00523

2.64

0.508

10.2

0.00453

2.49

0.569

11.4

0.00404

1.65

0.624

12.5

0.00369

1.14

0.685

13.7

0.00336

0.790

0.746

14.9

0.00308

3.95

0.955

19.1

0.00241

0.650

0.800

16.0

0.00287

0.572

0.855

17.1

0.00269

0.235

0.943

18.9

0.00244

0.179

1.02

20.4

0.00226

0.205

1.09

21.7

0.00212

0.148

1.15

23.1

0.00199

0.134

1.22

24.4

0.00188

0.148

1.28

25.7

0.00179

0.188

1.40

28.1

0.00164

0.107

1.48
0.178

29.6
3.55

0.0016
0.0675

17.2

0.305

6.11

0.0393

0.125

0.433

8.66

0.0277

0.130

0.568

11.4

0.0211

0.109

0.624

12.5

0.0192

0.140

0.688

13.8

0.0174

0.201

0.752

15.0

0.0160

0.161

0.809

16.2

0.01480.0699

0.865

17.3

0.01390.0681

0.957

19.1

0.01250.09291.04

20.7

0.0116

0.202

1.11

22.1

0.01080.08901.18

23.5

0.01020.06061.25

25.0

0.009620.04831.31

26.2

0.00915

0.06261.38

27.5

0.00872

0.09791.44

28.8

0.00834

0.092514

J. Phys. Commun. 5 (2021) 025012

A Capart et al

Table B2. TRS measurements for μa = 0.021 mm−1 (ρd = 20 mm). The unit of m¢s is mm−1.

m¢s

rd m¢sma m¢s

χ2

0.207

4.15

0.101

9.68

0.326

6.53

0.0643

0.305

0.472

9.44

0.0445

0.472

0.571

11.4

0.0368

0.162

0.624

12.5

0.0336

0.166

0.684

13.7

0.0307

0.190

0.743

14.9

0.0283

0.311

0.796

15.9

0.0264

0.109

0.849

17.0

0.0247

0.143

0.935

18.7

0.0225

0.173

1.01

20.2

0.02080.06251.07

21.5

0.01960.07921.14

22.8

0.0184

0.119

1.21

24.1

0.0174

0.107

1.27

25.3

0.0166

0.160

1.33

26.5

0.01580.07991.38

27.7

0.01520.07301.46

29.2

0.01440.0716

Table B3. Monte Carlo simulation for (Left) μa = 0.002 mm−1 and (Right) μa = 0.01 mm−1
(ρd = 20 mm). The unit of m¢s is mm−1.

m¢s

rd m¢sma m¢s

χ2

0.2

4.0

0.01

15.3

0.4

8.0

0.005

2.78

0.6

12.0

0.00333

0.554

0.8

16.0

0.0025

0.275

1.0

20.0

0.002

0.262

1.2

24.0

0.00167

0.345

1.4

28.0

0.00143

0.430

0.2

4.0

0.05

23.8

0.4

8.0

0.025

5.09

0.6

12.0

0.0167

1.36

0.8

16.0

0.0125

0.435

1.0

20.0

0.01

0.233

1.2

24.0

0.00833

0.250

1.4

28.0

0.007140.339

15

J. Phys. Commun. 5 (2021) 025012

A Capart et al

Table B4. Monte Carlo simulation for (Left) μa = 0.02 mm−1 and (Right) μa = 0.06 mm−1
(ρd = 20 mm). The unit of m¢s is mm−1.

m¢s

rd m¢sma m¢s

χ2

0.4

8.0

0.05

8.27

0.6

12.0

0.0333

2.38

0.8

16.0

0.025

0.860

1.0

20.0

0.02

0.449

1.2

24.0

0.0167

0.433

1.4
0.25

28.0
5.0

0.0143
0.24

62.3

0.3

6.0

0.2

46.3

0.35

7.0

0.171

35.4

0.4

8.0

0.15

28.2

0.45

9.0

0.133

22.9

0.5

10.0

0.12

18.2

0.75

15.0

0.08

8.18

Table B5. Monte Carlo simulation for μa = 0.1 mm−1 (ρd = 20 mm). The unit of m¢s is mm−1.

m¢s

rd m¢sma m¢s

χ2

0.25

5.0

0.4

117

0.35

7.0

0.286

67.9

ORCID iDs
Manabu Machida https://orcid.org/0000-0002-1830-0760
Yoko Hoshi https://orcid.org/0000-0001-5385-9714

References[1] Hielscher A H, Alcouffe R E and Barbour R L 1998 Comparison of ﬁnite-difference transport and diffusion calculations for photon
migration in homogeneous and heterogeneous tissues Phys. Med. Biol. 43 1285–302[2] Yoo K M, Liu F and Alfano R R 1990 When does the diffusion approximation fail to describe photon transport in random media? Phys.
Rev. Lett. 64 2647–50[3] Flock S T, Patterson M S, Wilson B C and Wyman D R 1989 Monte carlo modeling of light propagation in highly scattering tissues–i:
Model predictions and comparison with diffusion theory IEEE Trans. Biomed. Eng 36 1162–8[4] Bizheva K K, Siegel A M and Boas D A 1998 Path-length-resolved dynamic light scattering in highly scattering random media: thetransition to diffusing wave spectroscopy Phys. Rev. E 58 7664–7[5] Kop R H J, de Vries P, Sprik R and Lagendijk A 1997 Observation of anomalous transport of strongly multiple scattered light in thin
disordered slabs Phys. Rev. Lett. 79 4369–72[6] Elalouﬁ R, Carminati R and Greffet J-J 2002 Time-dependent transport through scattering media: from radiative transfer to diffusion.
J. Opt. A: Pure Appl. Opt. 4 S103–8[7] Elalouﬁ R, Carminati R and Greffet J-J 2004 Diffusive-to-ballistic transition in dynamic light transmission through thin scatteringslabs: a radiative transfer approach J. Opt. Soc. Am. A 21 1430–7[8] Xu H, Farrell T J and Patterson M S 2006 Investigation of light propagation models to determine the optical properties of tissue frominterstitial frequency domain ﬂuence measurements J. Biomed. Opt. 11 041104[9] Hielscher A H, Alcouffe R E and Barbour R L 1997 Transport and diffusion calculations on mri-generated data Optical Tomography andSpectroscopy of Tissue: Theory, Instrumentation, Model, and Human Studies II, Proc. SPIE 2979 (San Jose, CA,) ed B Chance andR R Alfano 500–8[10] Yoon G, Prahl S A and Welch A J 1989 Accuracies of the diffusion approximation and its similarity relations for laser irradiatedbiological media Appl. Opt. 15 2250–5[11] Freund I, Kaveh M and Rosenbluh M 1988 Dynamic multiple scattering: Ballistic photons and the breakdown of the photon-diffusion
approximation Phys. Rev. Lett. 60 1130–3[12] Patterson M S, Chance B and Wilson B C 1989 Time resolved reﬂectance and transmittance for the noninvasive measurement of tissue
optical properties Appl. Opt. 28 2331–6

16

J. Phys. Commun. 5 (2021) 025012

A Capart et al[13] Farrell T J, Patterson M S and Wilson B 1992 A diffusion theory model of spatially resolved, steady-state diffuse reﬂectance for thenoninvasive determination of tissue optical properties in vivo Med. Phys. 19 879–88[14] Pogue B W and Patterson M S 1994 Frequency-domain optical absorption spectroscopy of ﬁnite tissue volumes using diffusion theory
Phys. Med. Biol. 39 1157–80[15] Kienle A and Patterson M S 1997 Improved solutions of the steady-state and the time-resolved diffusion equations for reﬂectance froma semi-inﬁnite turbid medium J. Opt. Soc. Am. A 14 246–54[16] Hielscher A H, Jacques S L, Wang L and Tittel F K 1995 The inﬂuence of boundary conditions on the accuracy of diffusion theory intime-resolved reﬂectance spectroscopy of biological tissues Phys. Med. Biol. 40 1957–75[17] Martelli F, Contini D, Taddeucci A and Zaccanti G 1997 Photon migration through a turbid slab described by a model based on
diffusion approximation. ii. comparison with monte carlo results Appl. Opt. 36 4600–12[18] Kienle A, Lilge L, Patterson M S, Hibst R, Steiner R and Wilson B C 1996 Spatially resolved absolute diffuse reﬂectance measurements
for noninvasive determination of the optical scattering and absorption coefﬁcients of biological tissue Appl. Opt. 35 2304–14[19] Zhang X and Zhang Z-Q 2002 Wave transport through thin slabs of random media with internal reﬂection: ballistic to diffusivetransition Phys. Rev. E 66 016612[20] Chen B, Stamnes K and Stamnes J J 2001 Validity of the diffusion approximation in bio-optical imaging Appl. Opt. 40 6356–66[21] Durian D J and Rudnick J 1997 Photon migration at short times and distances and in cases of strong absorption J. Opt. Soc. Am. A 14
235–45[22] Spott T and Svaasand L O 2000 Collimated light sources in the diffusion approximation Appl. Opt. 39 6453–65[23] Ishimaru A 1978 Wave Propagation and Scattering in Random Media Vol. 1 (New York, NY: Academic Press)[24] Furutsu K and Yamada Y 1994 Diffusion approximation for a dissipative random medium and the applications Phys. Rev. E 50 3634–40[25] Bassani M, Martelli F, Zaccanti G and Contini D 1997 Independence of the diffusion coefﬁcient from absorption: experimental andnumerical evidence Opt. Lett. 22 853–5[26] Durduran T, Yodh A G, Chance B and Boas D A 1997 Does the photon-diffusion coefﬁcient depend on absorption? J. Opt. Soc. Am. A
14 3358–65[27] Elalouﬁ R, Carminati R and Greffet J-J 2003 Deﬁnition of the diffusion coefﬁcient in scattering and absorbing media J. Opt. Soc. Am. A
20 678–85[28] Egan W G and Hilgeman T W 1979 Optical Properties of Inhomogeneous Materials (New York, NY: Academic Press)[29] Groenhuis R A J, Ferwerda H A and Ten Bosch J J 1983 Scattering and absorption of turbid materials determined from reﬂection
measurements. 1: Theory Appl. Opt. 22 2456–62[30] Machida M and Nakamura G 2020 Born series for the photon diffusion equation perturbing the robin boundary condition J. Math.
Phys. 61 013502[31] Hoshi Y, Shimada M, Sato C and Iguchi Y 2005 Reevaluation of near-infrared light propagation in the adult human head: implications
for functional near-infrared spectroscopy J. Biomed. Opt. 10 064032[32] Flock S T, Wilson B C and Patterson M S 1989 Monte carlo modeling of light propagation in highly scattering tissues–ii: Comparison
with measurements in phantoms IEEE Trans. Biomed. Eng. 36 1169–73[33] van der Zee P and Delpy D T 1987 Simulation of the point spread function for light in tissue by a monte carlo technique Adv. Exp. Med.
Biol. 215 179–91[34] Okada E and Delpy T 2003 Near-infrared light propagation in an adult head model. ii. effect of superﬁcial tissue thickness on thesensitivity of the near-infrared spectroscopy signal Appl. Opt. 42 2915–22[35] Henyey L G and Greenstein J L 1941 Diffuse radiation in the galaxy Astrophys. J. 93 70–83[36] Guyon L, da Silva A, Planat-Chrétien A, Rizo P and Dinten J-M 2009 χ2 analysis for estimating the accuracy of optical properties
derived from time resolved diffuse-reﬂectance Opt. Exp. 17 20521–37[37] Alerstam E, Andersson-Engels S and Svensson T 2008 Improved accuracy in time-resolved diffuse reﬂectance spectroscopy Opt. Exp 1610434[38] Diop M and Lawrence K St. 2013 Improving the depth sensitivity of time-resolved measurements by extracting the distribution oftimes-of-ﬂight Biomed. Opt. Exp. 4 447–59[39] Guggenheim J A, Bargigia I, Farina A, Pifferi A and Dehghani H 2016 Time resolved diffuse optical spectroscopy with geometricallyaccurate models for bulk parameter recovery Biomed. Opt. Exp 7 3784–94[40] Bensoussan A, Lions J L and Papanicolaou G C 1979 Boundary layers and homogenization of transport processes Publ. RIMS 15 53–157[41] Larsen E W and Keller J B 1974 Asymptotic solution of neutron transport problems for small mean free paths. J. Math. Phys 15 75–81[42] Ryzhik L, Papanicolaou G and Keller J B 1996 Transport equations for elastic and other waves in random media Wave Motion 24
327–70[43] Arridge S R and Schotland J C 2009 Optical tomography: forward and inverse problems Inverse Prob. 25 123010

17