*Article* **Characteristic Length and Time Scales of the Highly Forward Scattering of Photons in Random Media**

#### **Hiroyuki Fujii \*,†, Moegi Ueno, Kazumichi Kobayashi and Masao Watanabe**

Division of Mechanical and Space Engineering, Faculty of Engineering, Hokkaido University, Kita 13 Nishi 8, Kita-ku, Sapporo, Hokkaido 060-8628, Japan; pueinsuki@eis.hokudai.ac.jp (M.U.);

kobakazu@eng.hokudai.ac.jp (K.K.); masao.watanabe@eng.hokudai.ac.jp (M.W.)

**\*** Correspondence: fujii-hr@eng.hokudai.ac.jp

† He partially conducted this research while he was a visiting scholar at University of Michigan.

Received: 28 November 2019; Accepted: 17 December 2019; Published: 20 December 2019

**Abstract:** Background: Elucidation of the highly forward scattering of photons in random media such as biological tissue is crucial for further developments of optical imaging using photon transport models. We evaluated length and time scales of the photon scattering in three-dimensional media. Methods: We employed analytical solutions of the time-dependent radiative transfer, *M*-th order delta-Eddington, and photon diffusion equations (RTE, dEM, and PDE). We calculated the fluence rates at different source-detector distances and optical properties. Results: We found that the zeroth order dEM and PDE, which approximate the highly forward scattering to the isotropic scattering, are valid in longer length and time scales than approximately 10/*μ <sup>t</sup>* and 40/*μ tv*, respectively, where *μ <sup>t</sup>* is the reduced transport coefficient and *v* the speed of light in a medium. The first and second order dEM, which approximate the highly forward-peaked phase function by the first two and three Legendre moments, are valid in the longer scales than approximately 4.0/*μ <sup>t</sup>* and 6.3/*μ tv*; 2.8/*μ t* and 3.5/*μ tv*, respectively. The boundary conditions less influence the length scales, while they reduce the times scales from those for bulk at the longer length scale than approximately 4.0/*μ t*. Conclusion: Our findings are useful for constructions of accurate and efficient photon transport models. We evaluated length and time scales of the highly forward scattering of photons in various kinds of three-dimensional random media by analytical solutions of the radiative transfer, *M*-th order delta-Eddington, and photon diffusion equations.

**Keywords:** radiative transfer equation; highly forward scattering of photons; diffusion and delta-Eddington approximations; characteristic length and time scales of photon transport

#### **1. Introduction**

Elucidation of the photon scattering and transport in random media such as biological tissue volumes is crucial for biomedical optical imaging using the near-infrared light in the wavelength range from 700 to 1100 nm such as diffuse optical tomography [1–3], because the imaging technique uses light scattered by the media. There are mainly three kinds of systems of the imaging techniques: time-domain, steady-state, and frequency-domain systems. Among them, the time-domain system has several advantages over steady-state and frequency-domain systems thanks to the richest information of the time-domain measurement data [4]. Time-dependent photon transport is basically classified into the two kinds of characteristic regimes (length and time scales): the ballistic and diffusive regimes corresponding to short and long source-detector (SD) distances and times after light incidence, respectively [5]. In the ballistic regime, photon is little scattered and almost goes straight. In the diffusive regime, photon is diffusive due to multiple scattering of photons, and the scattering process can be approximated as isotropic. Between the ballistic and diffusive regimes, photon undergoes a few scattering events in the highly forward direction and travels along the zig-zag paths while changing the

direction [6,7]. In this paper, we call this few-scattering-event regime between the ballistic and diffusive regimes as the scattering regime, so that we consider the three kinds of characteristic regimes for photon transport. The radiative transfer equation (RTE), which is the linear Boltzmann equation, can describe photon transport for random media in the three characteristic regimes because the phase function in the RTE can treat the highly forward scattering of photons. However, the numerical calculations of the RTE require high computational loads (times and memories) especially for three-dimensional (3D) random media because the RTE is an integro-differential equation with respect to the light intensity as a function of position vector, unit angular direction vector, and time. To reduce the computational loads, many researches have employed the photon diffusion equation (PDE), which is obtained by the diffusion approximation to the RTE and omits the angular direction. However, the PDE is valid only in the diffusive regime. A coupling model using the RTE and PDE has been constructed [8,9]. For the coupling model, the RTE is calculated in the ballistic and scattering regimes, while the PDE is calculated in the diffusive regime. For constructions of the coupling model or appropriate uses of the PDE, many researches have devoted their work to evaluations of the crossover length and time from the scattering to diffusive regimes [9–15]. For example, a comparison study [10] between actual light measurements and the PDE-results for slab media showed that the crossover length is evaluated as approximately 10/*μ <sup>t</sup>*, where the reduced transport coefficient *μ <sup>t</sup>* is a sum of reduced scattering coefficient and absorption coefficient. A numerical study [9] using the time-dependent RTE and PDE for 2D media evaluated the crossover length the same as the above study [10], and in addition, evaluated the crossover time as approximately 20/*μ tv*, where *v* is the speed of light in the medium. Although the coupling model can provide accurate and efficient results, a further efficiency improvement of the photon transport model is necessary because the RTE-calculations in the scattering regime still require high computational loads. For modeling the highly forward scattering of photons, the phase function in the RTE is highly forward-peaked and changes exponentially with respect to the scattering angle. Basically, accurate numerical treatments of the highly forward-peaked phase function require a large number of discrete angular directions, resulting in high computational loads of the RTE-calculations. To overcome the difficulty, one needs to investigate the highly forward scattering of photons and influences of the highly forward-peaked phase function on the RTE-results in the scattering regime. For this purpose, in this paper, we employ the *M*-th order delta-Eddington equation (dEM), which is obtained by the delta-Eddington approximation (dEA) to the RTE.

The dEA decomposes the highly forward-peaked phase function into purely forward-peaked and other components [16,17]. The purely forward-peaked component is given by the delta-function and contributes to a modification of the scattering length scale. The other components are expanded by finite series of Legendre polynomials and used as the phase function in the dEM. The dEM has been widely employed in the field of biomedical optics since Klose and coworkers have introduced from the field of astrophysics [18–22]. It has been shown that the numerical solutions of the dEM agree with those of the RTE while the number of the discrete angular directions is reduced to approximately one-sixth of the required number for the RTE-calculations. Nevertheless, a validity of the dEM is still unclear; Klose and coworkers showed that the second order dEM can provide the accurate results the same as the RTE [21], while Jia and coworkers stated that the zeroth order dEM is sufficient [22], and Boulet and coworkers employed the zeroth and first order dEM [20]. These differences probably come from the fact that the validity of the dEM depends on the SD distances, times after light incidence, and the optical properties, like the validity of the PDE.

Our objective is to examine photon transport especially in the scattering regime for various kinds of random media by using the time-dependent RTE, dEM, and PDE, and to evaluate a regime (length and time scales) for the dEM to be valid as a function of *μ <sup>t</sup>* and *μ tv*. We employed the analytical solutions instead of the numerical solutions because the analytical solutions are attained much faster than the numerical solutions and allow us to investigate a wide range of the SD distances and optical properties. In addition, we can focus the modification of the scattering coefficient and the approximation of the phase function by the dEA without a discussion of numerical errors induced by the numerical discretization. To the best of our knowledge, no studies have reported to evaluate the regime for the dEM. In this paper, we investigate photon transport for infinite and semi-infinite 3D homogeneous media as a first step of investigations for realistic heterogeneous biological tissue volumes. The investigations for infinite and semi-infinite 3D homogeneous media have been widely discussed for the diffuse optics community [4]. The analytical solutions have been employed for evaluations of the optical properties of various kinds of random media such as biological tissue volumes, colloidal suspensions [23], and agricultural products [24]. The following section describes the three kinds of the photon transport models for 3D random media: the time-dependent RTE, dEM, and PDE; and analytical solutions of the three models for the temporal profiles of the fluence rate. Section 3 investigates the results of the analytical solutions and evaluates the regime for the dEM. Finally, conclusions are described.

#### **2. Materials and Methods**

#### *2.1. Time-Domain Photon Transport Models*

#### 2.1.1. Radiative Transfer Equation (RTE)

We consider photon transport in the length and time scales of the mean free path and mean time of flight. Then, photon transport can be described by the RTE [25]. For 3D random media, the time-dependent RTE is given by

$$\left[\frac{\partial}{\nu\partial t} + \boldsymbol{\Omega} \cdot \nabla + \mu\_s(\boldsymbol{r}) + \mu\_s(\boldsymbol{r})\right] l(\boldsymbol{r}, \boldsymbol{\Omega}, t) = \mu\_s(\boldsymbol{r}) \int\_{\mathbb{S}^2} d\boldsymbol{\Omega}' p(\boldsymbol{\Omega}, \boldsymbol{\Omega}') l(\boldsymbol{r}, \boldsymbol{\Omega}', t) + q(\boldsymbol{r}, \boldsymbol{\Omega}, t), \qquad (1)$$

where, *I*(*r*, **Ω**, *t*) in W cm−<sup>2</sup> sr−<sup>1</sup> represents the light intensity as a function of spatial position vector *<sup>r</sup>* = (*x*, *<sup>y</sup>*, *<sup>z</sup>*) <sup>∈</sup> <sup>R</sup><sup>3</sup> in cm; angular direction (unit direction vector) **<sup>Ω</sup>** = (Ω*x*, <sup>Ω</sup>*y*, <sup>Ω</sup>*z*) <sup>∈</sup> <sup>S</sup><sup>2</sup> in sr; and time *t* in ps. *μa*(*r*) and *μs*(*r*) in cm−<sup>1</sup> are the absorption and scattering coefficients, respectively; *v* is the speed of light in the medium; *p*(**Ω**, **Ω** ) in sr−<sup>1</sup> is the phase function with **Ω** and **Ω** denoting the scattered-out and -in directions, respectively; and *q*(*r*, **Ω**, *t*) in W cm−<sup>3</sup> sr−<sup>1</sup> is a source function.

#### 2.1.2. Henyey-Greenstein Phase Function and Anisotropy Factor

For a formulation of *p*(**Ω**, **Ω** ), the Henyey-Greenstein (HG) phase function [26] is widely employed in biomedical optics to model the forward scattering of photons:

$$p\_{HG}(\boldsymbol{\Omega} \cdot \boldsymbol{\Omega}') = \frac{1}{4\pi} \frac{1 - \boldsymbol{g}^2}{\left(1 + \boldsymbol{g}^2 - 2\boldsymbol{g}\boldsymbol{\Omega} \cdot \boldsymbol{\Omega}'\right)^{3/2}} \tag{2}$$

where *g* ∈[−1, 1] is the anisotropy factor, defined by S<sup>2</sup> *d***Ω** (**Ω** · **Ω** )*p*(**Ω**, **Ω** ). Mostly, *g*-values of biological tissue volumes are larger than 0.8 [27], meaning the highly forward scattering. For examples, the *g*-value of muscle or bone is considered as 0.9; and the *g*-values of the human liver and of blood vessels are 0.955 [28] and 0.992 [29], respectively. Meanwhile, there are void regions in the human body such as the trachea in the human neck, whose *g*-value is almost zero, meaning the isotropic scattering. Figure 1 shows the HG phase function (Equation (2)) as a function of **Ω** · **Ω** ∈ [−1, 1] in a logarithmic scale. At *g* = 0.0 (isotropic scattering), the phase function is constant over the whole region of **Ω** · **Ω** . Meanwhile, as the *g*-value approaches to unity from zero (enhancement of the highly forward scattering), the exponential change of *pHG* with respect to **Ω** · **Ω** becomes large and the peak of *pHG* around **Ω** · **Ω** = 1.0 becomes sharp.

**Figure 1.** Henyey-Greenstein phase function as a function of **Ω** · **Ω** for isotropic scattering (*g* = 0.0), moderately forward scattering (*g* = 0.6), and highly forward scattering (*g* = 0.9 and 0.95), in a logarithmic scale.

The HG phase function is expanded in the infinite series of (unassociated) Legendre polynomials *Pl* (*l* = 0, 1, ··· , ∞):

$$p\_{HG}(\boldsymbol{\Omega} \cdot \boldsymbol{\Omega}') = \sum\_{l=0}^{\infty} \frac{2l+1}{4\pi} g^l P\_l(\boldsymbol{\Omega} \cdot \boldsymbol{\Omega}'),\tag{3}$$

where *g<sup>l</sup>* is the *l*-th order expansion coefficient.

#### 2.1.3. *M*-th Order Delta-Eddington Equation (dEM)

The dEA decomposes the highly forward-peaked phase function into a purely forward-peaked component expressed by the delta function and the other components [16,17]:

$$p\_\delta^M(\boldsymbol{\Omega} \cdot \boldsymbol{\Omega}') = \frac{1}{2\pi}h\delta(1 - \boldsymbol{\Omega} \cdot \boldsymbol{\Omega}') + (1 - h)p\_{\delta2}^M(\boldsymbol{\Omega} \cdot \boldsymbol{\Omega}'),\tag{4}$$

where *h* is a coefficient of the decomposition. *p<sup>M</sup> <sup>δ</sup>*<sup>2</sup> is a phase function excluding the delta-function component and expanded in the finite series of Legendre polynomials up to the order of *M*:

$$p\_{\delta\Omega}^M(\mathbf{\Omega}\cdot\mathbf{\Omega'}) = \sum\_{l=0}^M \frac{2l+1}{4\pi} \sigma\_{dEI} P\_l(\mathbf{\Omega}\cdot\mathbf{\Omega'}).\tag{5}$$

In the case of the HG phase function, the expansion coefficient *σdEl* and the decomposition coefficient *h* are determined so as to satisfy the moment conditions up to the order of *M* + 1:

$$
\sigma\_{d\to l} := \int\_{\mathbb{S}^2} d\Omega' p\_{\delta\2}^M(\mathbf{\Omega} \cdot \Omega') \mathbb{P}\_l(\mathbf{\Omega} \cdot \Omega') = \frac{\mathbf{g}^l - h}{1 - h}, \quad h = \mathbf{g}^{M+1}.\tag{6}
$$

The case of *l* = 0 in Equation (6) corresponds to the normalization condition of *p<sup>M</sup> δ*2: S<sup>2</sup> *d***Ω** *p<sup>M</sup> <sup>δ</sup>*2(**Ω** · **Ω** ) = 1. By using *p<sup>M</sup> <sup>δ</sup>* (Equation (4)), the RTE (Equation (1)) is approximated to

$$\left[\frac{\partial}{\partial t} + \boldsymbol{\Omega} \cdot \nabla + \mu\_{\boldsymbol{\theta}}(\boldsymbol{r}) + \mu\_{\boldsymbol{s}}^{M}(\boldsymbol{r})\right] I(\boldsymbol{r}, \boldsymbol{\Omega}, t) = \mu\_{\boldsymbol{s}}^{M}(\boldsymbol{r}) \int\_{\mathbb{S}^{2}} d\boldsymbol{\Omega}' p\_{\partial \boldsymbol{\Omega}}^{\boldsymbol{M}}(\boldsymbol{\Omega} \cdot \boldsymbol{\Omega}') I(\boldsymbol{r}, \boldsymbol{\Omega}', t) + q(\boldsymbol{r}, \boldsymbol{\Omega}, t), \quad (7)$$

where *μ<sup>M</sup> <sup>s</sup>* (*r*)=(<sup>1</sup> <sup>−</sup> *<sup>h</sup>*)*μs*(*r*)=(<sup>1</sup> <sup>−</sup> *<sup>g</sup>M*+1)*μs*(*r*). In this paper, Equation (7) is called as the *<sup>M</sup>*-th order delta-Eddington equation (dEM). Comparing the RTE (Equation (1)) with the dEM, *μ<sup>s</sup>* is modified to *μ<sup>M</sup> <sup>s</sup>* and *p* is approximated to *p<sup>M</sup> <sup>δ</sup>*2. The modification of *<sup>μ</sup><sup>s</sup>* to *<sup>μ</sup><sup>M</sup> <sup>s</sup>* means a change of the scattering length scale. Figure 2a,b plot a scattering length ratio of *μ<sup>M</sup> <sup>s</sup>* /*μ<sup>s</sup>* and a relative error *E<sup>M</sup> phase* of the phase functions between the HG function (Equation (2)) and the *M*-th order dEA (Equation (5)) at different *g*-values of 0.3, 0.6, 0.9, and 0.95. Here, the error *E<sup>M</sup> phase* is defined by

$$E\_{phase}^{M} = \int\_{-1}^{1} d\mu \left| \frac{p\_{\mathcal{S}2}^{M}(\mu) - p\_{HG}(\mu)}{p\_{HG}(\mu)} \right|, \quad \mu = \boldsymbol{\Omega} \cdot \boldsymbol{\Omega}'. \tag{8}$$

As the expansion order *M* increases, the ratio *μ<sup>M</sup> <sup>s</sup>* /*μ<sup>s</sup>* approaches to unity and the error *E<sup>M</sup> phase* decreases, meaning the scattering length scale and the phase function for the dEM converge to those for the RTE. At *g* = 0.9 and 0.95 (highly forward scattering), the convergences are slower than the cases of *g* = 0.3 and 0.6 (moderate forward scattering).

**Figure 2.** (**a**) Ratios *μ<sup>M</sup> <sup>s</sup>* /*μ<sup>s</sup>* and (**b**) relative errors *E<sup>M</sup> phase* of the phase functions between the HG function (Equation (2)) and the *M*-th order dEA (Equation (5)) at different *g*-values of 0.3, 0.6, 0.9 and 0.95. The error *E<sup>M</sup> phase* is defined by Equation (8).

#### 2.1.4. Photon Diffusion Equation (PDE)

The PDE is obtained by the diffusion approximation to the RTE [30]:

$$
\left[\frac{\partial}{\partial t} - \nabla \cdot D(r) \nabla + \mu\_d(r)\right] \Phi(r, t) = q\_{DE}(r, t), \tag{9}
$$

where the fluence rate Φ(*r*, *t*) is defined by <sup>S</sup><sup>2</sup> *d***Ω***I*(*r*, **Ω**, *t*); the diffusion coefficient *D*(*r*) is by [3(1 − *g*)*μs*(*r*))]−1; and the isotropic source *qDE*(*r*, *t*) is by <sup>S</sup><sup>2</sup> *d***Ω***q*(*r*, **Ω**, *t*).

*2.2. Analytical Solutions (AS) of the Time-Domain Photon Transport Models*

#### 2.2.1. 3D Infinite Homogeneous Media

In this paper, we mainly employed the analytical solutions of the time-dependent RTE, dEM, and PDE for 3D infinite homogeneous media to calculate the temporal profiles of the fluence rate Φ(*r*, *t*) when an isotropic point source is incident at the initial time and at the origin of the *xyz*-coordinate as shown in Figure 3a. Then, the source-detector (SD) distance *rsd* is given by *r* = |*r*|. The analytical form of the fluence rate for the RTE Φ*RTE*(*r*, *t*) with an arbitrary *g*-value is derived by Liemert and Kienle [31]:

$$\Phi\_{\rm RTE}(r, t) = \Phi\_{\rm LK}(r, t; \mu\_{\rm u}, \mu\_{\rm s}, \sigma\_{\rm HG}, n), \tag{10}$$

where Φ*LK* represents the analytical form by Liemert and Kienle; *σHG* the vector of the expansion coefficients for the HG phase function up to a maximum order of *<sup>N</sup>*: [1, *<sup>g</sup>*, *<sup>g</sup>*2, *<sup>g</sup>*3, ··· , *<sup>g</sup>N*]; the *<sup>N</sup>*-value is determined so as to attain sufficient convergence of the results; and *n* is the refractive index of the medium. The verifications of Φ*RTE*(*r*, *t*) were confirmed by several independent works [23,32]. Using Φ*LK*, we can obtain the analytical form for the dEM Φ*dEM*(*r*, *t*):

$$\Phi\_{dEM}(r,t) = \Phi\_{LK}(r,t; \mu\_{a\_{\prime}} \mu\_{s}^{M}, \sigma\_{\mathcal{M}}, n), \tag{11}$$

where *σ<sup>M</sup>* represents the vector of the expansion coefficients for *p<sup>M</sup> <sup>δ</sup>*<sup>2</sup> (Equation (6)) up to *N*(> *M*): [1, *σdE*1, *σdE*2, ··· *σdEM*, *σdEM*+<sup>1</sup> = 0, ··· , *σdEN* = 0]. We confirmed the verification of Φ*dEM*(*r*, *t*) by comparing with the numerical solutions in Appendix A. The analytical solution for the dE0 (*M* = 0) by Liemert and Kienle is largely oscillated and unstable for several cases such as short SD distances, short times after light incidence, and the low scattering coefficient. In those cases, we employed

the analytical solution for the dE0 derived by Paasschens [33] instead. It is noted that Martelli and coworkers employed a heuristic approach to obtain the approximate solution of the RTE for highly forward scattering [34] from the exact solution of the RTE for isotropic scattering by Paasschens [33]. Their approach is identical to the zeroth order dEA, so that the Φ(*r*, *t*)-result using their approach coincides with the result of Φ*dE*0(*r*, *t*). The analytical form for the PDE Φ*PDE*(*r*, *t*) is given as [35]

$$\Phi\_{PDE}(r,t) = \frac{1}{\left(4\pi Dt\right)^{3/2}} \exp\left[-\mu\_a vt - \frac{r^2}{4Dvt}\right].\tag{12}$$

**Figure 3.** Source and detector positions in 3D (**a**) infinite and (**b**) semi-infinite homogeneous media.

#### 2.2.2. 3D Semi-Infinite Homogeneous Media

To investigate influences of the boundary conditions on the temporal profiles of the fluence rate, we employed the approximate solutions of the RTE, dEM, and PDE for 3D semi-infinite homogeneous media under the diffusive refractive-index mismatched boundary condition:

$$
\nabla \Phi(\mathbf{r}, t) + D\gamma(\mathbf{u}) \mathbf{e}\_z \cdot \nabla \Phi(\mathbf{r}, t) = 0,\tag{13}
$$

where *e<sup>z</sup>* represents the outward normal unit vector at the *xy*-plane in Figure 3b, the coefficient *<sup>γ</sup>*(*n*) is give by 2(<sup>1</sup> + *Rf*)/(<sup>1</sup> − *Rf*) with *Rf* = −1.44*n*−<sup>2</sup> + 0.71*n*−<sup>1</sup> + 0.69 + 0.064*<sup>n</sup>* [36]. Using the extrapolation method for the boundary condition (Equation (13)) [37], the analytical solutions of the RTE, dEM, and PDE for semi-infinite homogeneous media are obtained:

$$\Phi\_{\rm seci}(\mathbf{r}, t; z\_0, z\_b) = \Phi\_{\rm inf}(|\mathbf{r} - \mathbf{r}^+|, t) - \Phi\_{\rm inf}(|\mathbf{r} - \mathbf{r}^-|, t), \tag{14}$$

where *r* = (*ρ*, *z*), *ρ* = (*x*, *y*), a real source position of *r*<sup>+</sup> = (0, 0, *z*0) with *z*<sup>0</sup> = 1/*μ <sup>s</sup>*, an imaginary source position of *r*<sup>−</sup> = (0, 0, −2*zb* − *z*0) with *zb* = *γ*(*n*)*D*, and Φ*inf* = {Φ*RTE*, Φ*dEM*, Φ*PDE*} (Equations (10)–(12)). Then, the SD distance *rsd* is given by <sup>|</sup>*<sup>r</sup>* <sup>−</sup> *<sup>r</sup>*+|. As a preliminary study, we compared the Φ*semi*(*r*, *t*)-results for the RTE with those by another diffusion approximation method for the boundary condition [38], and confirmed that both results are almost the same to each other. The verification of the approximate solutions of the RTE has been confirmed by a comparative study with Monte-Carlo simulations [38].

#### *2.3. Optical Properties, Source-Detector (SD) Distances, and Computations of the Analytical Solutions*

We investigated photon transport in various kinds of random media at different SD distances. Table 1 lists parameter sets of the optical properties for the media and the SD distances, where the optical properties are in the wavelength range of the near-infrared light. At each parameter set, four parameters are fixed and one parameter is varied, e.g., at the parameter set A, *μa*, *μs*, *g*, and *n* are fixed and *rsd* is varied from 0.15 to 10.0 cm. The optical properties at the parameter set A are often used in the research field of biomedical optics [34,38], and the properties at the parameter set B correspond to

those of SiO2 with epoxy resin [39]. Investigations at the parameter sets C, D, E and F allow us to find the dependence of the analytical results for the fluence rate on *μa*, *μs*, and *g*, respectively.

**Table 1.** Optical properties: *<sup>μ</sup>a*[cm−1], *<sup>μ</sup>s*[cm−1], *<sup>g</sup>*[−], and *<sup>n</sup>*[−]; and SD distance *rsd*[cm]. Corresponding colors are used in Figures 5, 6 and 8–11.


For the numerical calculations of the analytical solutions of the RTE, dEM, and PDE, we employed MATLAB codes. For the calculation of Φ*LK*, we modified an open source by Liemert and Kienle [31]. The computational times for analytical solution of the RTE and dEM are less than 30 s when the optical properties and SD distance are given. It is noted that the calculations of the numerical solutions of the RTE and dEM require much higher computational loads than those of the analytical solutions. For example, the computational times to obtain the numerical solutions of the dEM discussed in Appendix A are roughly 26 h, although the parallel computing was carried out.

#### **3. Results**

#### *3.1. Classification of the Length and Time Scales of Photon Transport*

In this subsection, we investigated the temporal profiles of the fluence rate Φ(*r*, *t*) using the analytical solutions of the RTE (Equation (10)), dE0 (Equation (11) with *M* = 0), and PDE (Equation (12)) for 3D infinite homogeneous media. Also, we investigated the peak time of Φ(*r*, *t*), *tpeak*, because it is one of characteristic times of photon transport, and a difference in Φ(*r*, *t*) is correlated with a difference in *tpeak*. Based on results of *tpeak* at different SD distances *rsd* = *r*, we classify length and time scales of photon transport into following three kinds of characteristic regimes; (1) ballistic regime: *tpeak* ∼ *t ball peak*, (2) scattering regime (few-scattering-event regime): *tpeak t ball peak*, *t di f f peak*, and (3) diffusive regime: *tpeak* ∼ *t di f f peak*, where *t ball peak* and *t di f f peak* are defined by

$$t\_{peak}^{ball} = \frac{r}{\upsilon} \, \, \tau \qquad \qquad t\_{peak}^{diff} = \frac{-3 + 2\sqrt{\frac{9}{4} + \frac{r^2 \mu\_a}{D}}}{4\mu\_a \upsilon} \, \, \, \tag{15}$$

*t ball peak* represents the arriving time of photons to the detector position without scattering, and *t di f f peak* the peak time for the diffusive photons calculated from the PDE (Equation (12)), respectively. Figure 4 show examples of the temporal profiles of Φ(*r*, *t*) at the three characteristic length regimes. Here, we used the analytical solutions with the optical properties for the parameter set A as listed in Table 1 corresponding to the highly forward scattering of photons. While the RTE can appropriately treat the highly forward scattering, the dE0 and PDE approximate the highly forward scattering to the isotropic scattering. As shown in Figure 4a corresponding to the ballistic regime, the Φ-profiles using the RTE and dE0 are sharp, although both profiles are different from each other. The sharp profiles mean that photons are little scattered and travel almost straight. The peak times *tpeak* for the RTE and dE0 are almost the same as the arriving time *t ball peak*, whose value is approximately 8.4 ps. Meanwhile, the broad profile of Φ using the PDE is observed, indicating the unphysical multiple scattering of photons. In addition, the *tpeak*-value for the PDE is shorter than the *t ball peak*-value, meaning the violation of the causality in the PDE. As shown in Figure 4b corresponding to the scattering regime, the profiles using all the equations have smooth peaks. These peaks are formed by scattering of photons, so that these

peak times are longer than *t ball peak*, in other words, the paths of photons are longer than the SD distance. In Figure 4b, the peak times for the dE0 and PDE are different from that for the RTE, indicating that the forward scattering of photons is not approximated to the isotropic scattering in this regime. As shown in Figure 4c corresponding to the diffusive regime, all the profiles are almost the same to each other, meaning the photon scattering becomes diffusive and the isotropic scattering approximation holds.

**Figure 4.** Analytical solutions of the RTE (Equation (10)), dE0 (Equation (11) with *M* = 0), and PDE (Equation (12)) at the three characteristic length regimes of photon transport: (**a**) ballistic regime at the SD distance of *r* = 0.18 cm (log10 *rμ <sup>t</sup>* = 0.26), (**b**) scattering regime at *r* = 0.40 cm (log10 *rμ <sup>t</sup>* = 0.60), and (**c**) diffusive regime at *r* = 1.60 cm (log10 *rμ <sup>t</sup>* = 1.20). The optical properties at the parameter set A are used, which are listed in Table 1.

We evaluated crossover lengths and times from the ballistic to scattering regimes, and from the scattering to diffusive regimes, respectively, by investigating *tpeak* for the RTE, dE0, and PDE. Then, we normalized the SD distance *r* by *μ <sup>t</sup>* and the peak time *tpeak* by *μ tv*, respectively, with the reduced transport coefficient of *μ <sup>t</sup>* = *μs*(1 − *g*) + *μa*. The inverses of *μ <sup>t</sup>* and *μ tv* correspond to the characteristic length and time for the dE0 and PDE where photons experience single scattering and absorbing processes. The normalizations allow us to investigate general features of photon transport in various kinds of random media. Figure 5 shows the results of the normalized peak times as a function of the normalized SD distances for the parameter sets A to F listed in Table 1 in a regime of short length and time scales. Here, the colors of the figures correspond to the parameter sets, e.g., the black colored pluses, circles, and squares correspond to the results for the parameter set A. In the regime of log10 *rμ <sup>t</sup>* - 0.35 and log10 *tpeakμ tv* - 0.35, the *tpeak*-results for the RTE and dE0 are linearly related to *r*, while the results for the PDE are not. This result suggests that the regime of the length and time scales correspond to the ballistic regime; and the crossover length and time from the ballistic to scattering regimes are evaluated as approximately 100.35/*μ <sup>t</sup>* ∼ 2.2/*μ <sup>t</sup>* and 100.35/*μ tv* ∼ 2.2/*μ tv*, respectively.

Figure 6a shows the *tpeak*-results in a regime of longer length and time scales than the ballistic regime. In the case of the same parameter set (the same colored figures), the *tpeak*-results for the dE0 and PDE deviate from those for the RTE in the regime of 0.35 log10 *rμ <sup>t</sup>* - 1.0 and 0.35 - log10 *tpeakμ tv* - 1.6, while the results for all the equations are almost the same to each other in the regime of 1.0 log10 *rμ <sup>t</sup>* and 1.6 log10 *tpeakμ tv*. For the further investigation of *tpeak*, we evaluated the relative errors *Epeak* of *tpeak* against the RTE, defined by

$$E\_{\rm peak} = \left| \frac{t\_{\rm peak} - t\_{\rm peak,RTE}}{t\_{\rm peak,RTE}} \right| \times 100, \qquad t\_{\rm peak} = \{t\_{\rm peak,dE0}, t\_{\rm peak,PDE}\}. \tag{16}$$

Although we calculated *Epeak* for all the parameter sets, we showed the results only for the parameter set A in Figure 6b because the results for the other parameter sets behave similarly to the results for the parameter set A. As shown in Figure 6b, the *Epeak*-values for the dE0 and PDE are larger than 2% in the regime of 0.35 log10 *rμ <sup>t</sup>* - 1.0. Here the *Epeak*-value of 2% is considered as a thresh to determine whether the dE0 and PDE are valid or not. The large *Epeak*-values indicate the strong influence of modeling the forward scattering on the Φ-results in this regime. Meanwhile, the *Epeak*-values for the dE0 and PDE are less than 2% in the regime of 1.0 log10 *rμ <sup>t</sup>*, meaning the isotropic scattering approximation holds. From these results, the crossover length and time from the scattering to diffusive regimes are evaluated as approximately 101/*μ <sup>t</sup>* = 10/*μ <sup>t</sup>* and 101.6/*μ tv* ∼ 40/*μ tv*, respectively. The crossover length from the scattering to diffusive regimes has been extensively discussed so far, and our evaluation is consistent with the previous studies such as a comparative study of the measurement data and the PDE-results for 3D slab media [10], and a comparative study of the numerical results using the RTE and PDE for 2D square media [9]. This consistency implies that boundary conditions less influence the crossover length because the previous studies consider the finite media, while this study the infinite media. Meanwhile, the crossover time from the scattering to diffusive regimes is longer than that evaluated by the previous study [9], which is approximately 20/*μ tv*. The difference probably comes from the dimension of the medium and boundary effects; the previous study considers 2D finite media, while this study the 3D infinite media.


**Figure 5.** Normalized peak times *tpeak* by *μ <sup>t</sup>v* as a function of the normalized SD distances *r* by *μ t* for the RTE, dE0, and PDE in the ballistic and scattering regimes. *μ <sup>t</sup>* = *μs*(1 − *g*) + *μ<sup>a</sup>* represents the reduced transport coefficient and *v* the speed of light in a medium, respectively. Colors of the figures correspond to the results for the parameter sets A to F listed in Table 1 or in the legend at the right-side of the figure. The dashed line represents the linear relation between *tpeak* and *r*. The yellow colored area corresponds to the scattering regime.

**Figure 6.** (**a**) Normalized peak times and (**b**) relative errors *Epeak* (Equation (16)) for the dE0 and PDE as a function of the normalized SD distances in the scattering and diffusive regimes. The other details are the same as Figure 5.

The dependence of the *tpeak*-results on the parameter sets is observed in the diffusive regime as shown in Figure 6a. This is because that in the diffusive regime, the normalizations of *tpeak* by *μav* and of *r*<sup>2</sup> by *μa*/*D* are appropriate, which are suggested from the form of *t di f f peak* (Equation (15)). As a preliminary study, we confirmed that when the above normalizations are employed, the *tpeak*-results for all the parameter sets are collapsed on a single curve in the diffusive regime. However, for the evaluation of the crossover length and time from the scattering to diffusive regimes, the above normalizations are not appropriate because the length regime, where the *Epeak*-values are larger than 2%, depends on the parameter sets, so that we can not evaluate the crossover length uniquely for all the parameter sets.

#### *3.2. Length and Time Scales for the dEM to Be Valid in the Scattering Regime*

In this subsection, we investigated the temporal profiles of the fluence rate Φ(*r*, *t*) and their peak times *tpeak* for the RTE and dEM in the scattering regime for the purpose of evaluations of the length and time scales for the dEM to be valid. As mentioned in Section 2.1.3, for the dEM, the scattering length scale is modified to the inverse of *μ<sup>M</sup> <sup>s</sup>* from the inverse of *μs*, and the highly forward scattering of photons described by *pHG* is approximated to *p<sup>M</sup> <sup>δ</sup>*2. These modifications influence the results of Φ(*r*, *t*).

Firstly, we calculated Φ(*r*, *t*) and *tpeak* for the dEM at different *g*-values of 0.3, 0.6, 0.9, and 0.95 in the scattering length regime, and investigated the dependence of the results of Φ(*r*, *t*) and *tpeak* on the expansion order *M* and on the *g*-value. The optical properties and SD distance are set as *μ<sup>a</sup>* = 0.10 cm<sup>−</sup>1, *μ <sup>s</sup>* = *<sup>μ</sup>s*(<sup>1</sup> − *<sup>g</sup>*) = 10.0 cm<sup>−</sup>1, *<sup>n</sup>* = 1.40, and *<sup>r</sup>* = 0.32 cm, where the normalized SD distance is constant as ∼ 100.5/*μ <sup>t</sup>*. Figure 7a shows the temporal profiles of Φ(*r*, *t*) for the RTE and dEM with *M* = 0, 1, and 2 (dE0, dE1, and dE2) at *g* = 0.9, corresponding to the highly forward scattering. While the Φ-profile for the dE0 largely deviates from that for the RTE, the profiles for the dE1 and dE2 are similar to that for the RTE. Moreover, we investigated the relative error *E<sup>M</sup> peak* of the peak time for the dEM:

$$E\_{p\text{calc}}^{M} = \left| \frac{t\_{p\text{calc},dEM} - t\_{p\text{ack},RTE}}{t\_{p\text{ack},RTE}} \right| \times 100,\tag{17}$$

where *tpeak*,*dEM* represents the peak time of Φ(*r*, *t*) for the dEM. As shown in Figure 7b, the *E<sup>M</sup> peak*-values are quite small when the expansion order *M* is larger than 2 at all the *g*-values. This result means that the dEM with *M* > 2 can provide the same accuracy as the RTE, almost independently of the *g*-values. Meanwhile, as shown in Figure 2, the *g*-dependences are observed in the ratios of *μ<sup>M</sup> <sup>s</sup>* /*μ<sup>s</sup>* and the errors of the phase function *E<sup>M</sup> phase*. These results suggest that for the accurate calculations of Φ(*r*, *t*) using the RTE and dEM, it is sufficient to satisfy a few orders moment conditions of the phase function and the exact agreement of the profiles of the phase functions is not required. It is noted that the analytical results of the small *E<sup>M</sup> peak*-values with *M* > 2 are quite different from the numerical calculations. The numerical results for the dEM are usually accurate within a finite range of *M*, and when the *M*-value is larger than the maximum value of the finite range, the numerical results for the dEM largely deviate from those for the RTE because of the numerical errors induced by angular discretization.

On the above discussion, we investigated *E<sup>M</sup> peak* for the dEM at the constant normalized SD distance. Next, we investigated *E<sup>M</sup> peak* for the dE0, dE1, and dE2 at the normalized SD distances varied from the scattering to diffusive regimes. Then, we evaluated the crossover lengths for the dE0, dE1, and dE2; in a regime of a smaller length scale than the crossover length, the *E<sup>M</sup> peak*-values are larger than 2%, corresponding to the invalidity of the dEM, while in a regime of a longer length scale, the errors are smaller than 2%, corresponding to the validity of the dEM. Figure 8a–c show the *E<sup>M</sup> peak*-values as a function of the normalized SD distances at the parameter sets A to F listed in Table 1. While the max value of *E<sup>M</sup> peak* for the dE0 is over 20%, the max values for the dE1 and dE2 are around 10% and 4%, respectively. The length regime for the dE2 where the *E<sup>M</sup> peak*-values are larger than 2% is shorter than those for the dE0 and dE1. From the results in the figures (a), (b), and (c), we evaluated the crossover lengths for the dE0, dE1, and dE2 as approximately 101/*μ <sup>t</sup>* = 10/*μ t*, 100.6/*μ <sup>t</sup>* ∼ 4.0/*μ t*, and 100.45/*μ <sup>t</sup>* ∼ 2.8/*μ <sup>t</sup>*, respectively. These results mean that the scattering regime is classified into the three kinds of regimes by the two kinds of the crossover lengths for the dE1 and dE2. Here, we

evaluated the crossover lengths for the dE0, dE1, and dE2 as a function of *μ <sup>t</sup>* for the purpose of the classifications of the scattering regime. As an alternative, we can consider the normalization of the SD distances by *μ<sup>M</sup> <sup>t</sup>* = *μ<sup>M</sup> <sup>s</sup>* + *μ<sup>a</sup>* because for the dEM, the inverse of *μ<sup>M</sup> <sup>t</sup>* is the characteristic length of single scattering and absorption processes. Figure 8d shows *E<sup>M</sup> peak* for the dE0, dE1, and dE2 as a function of the SD distances by *μ<sup>M</sup> <sup>t</sup>* . The crossover lengths for the dEM with *M* = 0, 1, and 2 are within the regime of 0.9 log10 *<sup>r</sup>μ<sup>M</sup> <sup>t</sup>* - 1.0, less dependently on the expansion order *M*. It is noted that the crossover length from the scattering and diffusive regimes depends on the expansion order *M* in the normalization by *μ<sup>M</sup> <sup>t</sup>* , so that we cannot evaluate the regime for the dEM to be valid uniquely in the normalization. The rough evaluation of the crossover length for the dEM as 10/*μ<sup>M</sup> <sup>t</sup>* means that after ten times scattering and absorption processes for the dEM, the *M*-th order dEA holds.

Figure 9 shows the normalized peak times for the RTE, dE0, dE1, and dE2 as a function of the normalized SD distances by *μ <sup>t</sup>*. From the figure and the results for the crossover lengths based on *μ t*, we evaluated the crossover times for the dE0, dE1, and dE2 as approximately 101.6/*μ tv* ∼ 40/*μ tv*, 100.8/*μ tv* ∼ 6.3/*μ tv*, and 100.55/*μ tv* ∼ 3.5/*μ tv*, respectively. Here, the evaluations of the crossover times are based on the *tpeak*-results for the RTE (dashed line in Figure 9 at the parameter set A) because the RTE is the most accurate photon transport model.

#### *3.3. Influence of the Boundary Conditions on the Crossover Lengths and Times*

In this subsection, we investigated the boundary effects on the fluence rate Φ(*r*, *t*) and on the crossover lengths and times for the dEM with *M* = 0, 1, and 2 by comparing the analytical solutions for the semi-infinite media with those for the infinite media. We considered the parameter sets A and B listed in Table 1, and varied the SD distances, denoted by *rsi* <sup>=</sup> <sup>|</sup>*<sup>r</sup>* <sup>−</sup> *<sup>r</sup>*+|, in the two directions: the *z*-direction with *r* = (0, 0, *z*) and *ρ*-direction with *r* = (*ρ*, 0) as shown in Figure 3b. For the case of the detector position in the *z*-direction, as the *z*-value increases, the detector position becomes far from the boundary, so that the boundary effects probably become weak. Meanwhile, for the case of the detector position in the *ρ*-direction, a distance between the detector and the boundary is unchanged, so that the boundary effects are also unchanged. Figure 10a–c show the relative errors *E<sup>M</sup> peak* for the dE0, dE1, and dE2 in the semi-infinite media as a function of the normalized SD distance by *μ <sup>t</sup>*. The crossover lengths for the dEM in the semi-infinite media are almost the same as those in the infinite media, meaning less influences of the boundary conditions on the crossover lengths.

**Figure 7.** (**a**) Temporal profiles of Φ(*r*, *t*) using the RTE, dE0, dE1, and dE2 at *g* = 0.9. (**b**) Relative errors *E<sup>M</sup> peak* of the peak time of Φ(*r*, *t*) (Equation (17) as a function of the expansion order *M* at different *g*-values of 0.3, 0.6, 0.9, and 0.95. The optical properties and SD distance are set as *μ<sup>a</sup>* = 0.10 cm<sup>−</sup>1, *μ <sup>s</sup>* <sup>=</sup> *<sup>μ</sup>s*(<sup>1</sup> <sup>−</sup> *<sup>g</sup>*) = 10.0 cm<sup>−</sup>1, *<sup>n</sup>* <sup>=</sup> 1.40, and *<sup>r</sup>* <sup>=</sup> 0.32 cm.

**Figure 8.** Relative errors *E<sup>M</sup> peak* (Equation (17)) for the (**a**) dE0, (**b**) dE1, and (**c**) dE2 as a function of the normalized SD distances by *μ <sup>t</sup>* at the parameter sets A to F listed in Table 1. The yellow, orange, and red colored areas correspond to the length regimes where the dE0, dE1, and dE2 are invalid, respectively. (**d**) Normalization of the SD distances by *μ<sup>M</sup> <sup>t</sup>* . The gray colored areas correspond to the crossover length for the dEM. The other details are the same as Figure 6.

**Figure 9.** Normalized peak times of Φ(*r*, *t*) for the RTE, dE0, dE1, and dE2 as a function of the normalized SD distances. The dashed line represents the connection line for the RTE-results at the parameter set A. The other details are the same as Figure 8a–c.

Figure 11 shows the normalized peak times for the RTE in the infinite and semi-infinite media as a function of the SD distances *rsd* for the parameter sets A and B. Here, the SD distances *rsd* are given by *r* in infinite media, and by *rsi* in semi-infinite media, respectively. The *tpeak*-results in the infinite media and semi-infinite media at the *z*-direction behave similarly to each other on the whole regime. Meanwhile, the *tpeak*-results in the semi-infinite media at the *ρ*-direction (dotted line in Figure 11 for the results at the parameter set A) are smaller than those in the other two cases (dashed line) especially on the regime of the longer length scale than approximately 100.6/*μ <sup>t</sup>* ∼ 4.0/*μ <sup>t</sup>*, which is almost the same as the crossover length for the dE1. These results suggest that the boundary conditions influence the crossover time for the dE0, while they have less influence the times for the dE1 and dE2.

**Figure 10.** Relative errors *E<sup>M</sup> peak* for the (**a**) dE0, (**b**) dE1, and (**c**) dE2 in the semi-infinite media (Figure 3b) as a function of the normalized SD distances *rsi* <sup>=</sup> <sup>|</sup>*<sup>r</sup>* <sup>−</sup> *<sup>r</sup>*+<sup>|</sup> at the parameter sets A and B listed in Table 1. The SD distances are varied in the *z*-direction with *r* = (0, 0, *z*) and in the *ρ*-direction with *r* = (*ρ*, 0). The other details are the same as Figure 8a–c.

**Figure 11.** Normalized peak times *tpeak* for the RTE in infinite media, and semi-infinite media at the *z*and *ρ*-directions. The SD distances *rsd* are given by *r* in infinite media and by *rsi* in semi-infinite media, respectively. The dashed and dotted lines represent the connection lines for the results in infinite and semi-infinite media at the *ρ*-direction, respectively, for the parameter set A. The other details are the same as Figure 10.

#### **4. Conclusions**

We investigated the temporal profiles of the fluence rate Φ(*r*, *t*) and their peak times *tpeak* for various kinds of random media at different SD distances by analytical solutions of the RTE, dE0, and PDE for 3D infinite homogeneous media. We evaluated the crossover lengths and times from the ballistic to scattering (few-scattering-event) regimes, and from the scattering to diffusive regimes as approximately 2.2/*μ <sup>t</sup>* and 2.2/*μ tv*; 10/*μ <sup>t</sup>* and 40/*μ tv*, respectively.

Next, we investigated Φ(*r*, *t*) and *tpeak* for the dEM mainly in the scattering regime. We found that the results for the dE0 are quite similar to the results for the PDE because the dE0 and PDE approximate the forward scattering of photons as the isotropic scattering. We also found that at the expansion order of *M* larger than 2, the results for the dEM are almost the same as those for the RTE in the scattering and diffusive regimes, meaning the validity of the dEM with *M* > 2 in the regimes. We evaluated the crossover lengths and times for the dE1 and dE2 as approximately 4.0/*μ <sup>t</sup>* and 6.3/*μ tv*; 2.8/*μ <sup>t</sup>* and 3.5/*μ tv*, respectively.

Finally, we investigated the boundary effects on the characteristic length and time scales of photon transport by comparing the results for the infinite media with those for the semi-infinite media. We found that the boundary conditions less influence on the crossover lengths for the dE0, dE1, and dE2, while the boundary conditions reduce the *tpeak*-values from those for bulk in the regime of the longer length scale than 4.0/*μ <sup>t</sup>*, which is the same as the crossover length for the dE1. Hence, the crossover

times for the dE1 and dE2 are less influenced by the boundary conditions. Our findings are useful for constructions of accurate and efficient photon transport models especially in the scattering regime.

**Author Contributions:** Conceptualization, H.F.; investigation and analysis, H.F. and M.U.; writing—original draft preparation, H.F.; writing—review and editing, H.F., K.K. and M.W.; supervision, K.K. and M.W.; All authors discussed the results and contributed to the final manuscript. All authors have read and agreed to the published version of the manuscript.

**Funding:** The first author (H.F.) acknowledges support from Grant-in-Aid for Scientific Research (18K13694) of the Japan Society for the Promotion of Science, KAKENHI.

**Acknowledgments:** The authors would like to thank Go Chiba, Yukio Yamada, and Yoko Hoshi for fruitful discussion. The first author (H.F.) appreciates Leung Tsang, his laboratory members, and staffs at University of Michigan for their kind supports on H.F.'s research work.

**Conflicts of Interest:** The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, or in the decision to publish the results.

#### **Appendix A. Verification of the Analytical Solutions of the dEM for 3D Infinite Homogeneous Media**

We verified the analytical solutions of the dEM with *M* = 0, 1, and 2 for 3D infinite homogeneous media (Equation (11)) by comparing with the numerical solutions for 3D homogeneous cubic media under the refractive-index mismatched boundary condition. We numerically solved the dEM and calculated the temporal profiles of the fluence rate Φ(*r*, *t*) based on the finite-difference and discrete ordinates methods. For spatial and temporal discretization, we employed the 3rd order weighted essentially non oscillatory and the 3rd order total variation diminishing-Runge-Kutta methods, respectively. For accurate treatments of the HG phase function and dE phase functions in a case of the highly forward scattering, we employed the Galerkin quadrature method. Source and detector positions are set inside the medium to suppress the boundary effects because we compare the analytical solutions for infinite media with the numerical solutions for finite media. For the details for the numerical calculations, refer to [32,40]. The optical properties of the medium are set as *μ<sup>a</sup>* = 0.20 cm<sup>−</sup>1, *μ<sup>s</sup>* = 100 cm<sup>−</sup>1, *g* = 0.90, and *n* = 1.40; and the SD distance is set as *r* = 0.40 cm, corresponding to the scattering regime.

Figure A1 compare the analytical and numerical solutions of the dEM with *M* = 0, 1, and 2 for Φ(*r*, *t*). We found good agreements between the analytical and numerical solutions for all the cases, meaning the verification of the analytical solutions of the dEM.

**Figure A1.** Comparisons of the analytical and numerical solutions of the dEM with (**a**) *M* = 0, (**b**) *M* = 1, and (**c**) *M* = 2 for the fluence rate Φ(*r*, *t*). The optical properties of the medium are set as *μ<sup>a</sup>* = 0.20 cm<sup>−</sup>1, *μ<sup>s</sup>* = 100 cm<sup>−</sup>1, *g* = 0.90, and *n* = 1.40; and the SD distance is set as *r* = 0.40 cm, corresponding to the scattering regime.

#### **References**


c 2019 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

MDPI St. Alban-Anlage 66 4052 Basel Switzerland Tel. +41 61 683 77 34 Fax +41 61 302 89 18 www.mdpi.com

*Applied Sciences* Editorial Office E-mail: applsci@mdpi.com www.mdpi.com/journal/applsci

MDPI St. Alban-Anlage 66 4052 Basel Switzerland

Tel: +41 61 683 77 34 Fax: +41 61 302 89 18