*N*-layered turbid media. The hybrid diffusion-P3 equation is solved for an

*N*-layered finite or infinite turbid medium in the steady-state domain for one point source using the extrapolated boundary condition. The Fourier transform formalism is applied to derive the analytical solutions of the fluence rate in Fourier space. Two inverse Fourier transform methods are developed to calculate the fluence rate in real space. In addition, the solutions of the hybrid diffusion-P3 equation are compared to the solutions of the diffusion equation and the Monte Carlo simulation. For the case of small absorption coefficients, the solutions of the

*N*-layered diffusion equation and hybrid diffusion-P3 equation are almost equivalent and are in agreement with the Monte Carlo simulation. For the case of large absorption coefficients, the model of the hybrid diffusion-P3 equation is more precise than that of the diffusion equation. In conclusion, the model of the hybrid diffusion-P3 equation can replace the diffusion equation for modeling light propagation in the

*N*-layered turbid media for a wide range of absorption coefficients.

## 1.

## Introduction

The *in vivo* determination of tissue optical properties and the study of light propagation in biological tissues are very important in a variety of biomedical fields that can be used to obtain knowledge of the physiological state of tissues.^{1, 2, 3, 4} Diffuse reflectance spectroscopy has emerged as one of the noninvasive optical techniques and it has been researched to obtain quantitative tissue characterization and optical properties.^{4, 5, 6, 7, 8, 9} In recent years, the study of light propagation in layered turbid media has gathered much attention because many parts of the body such as the skin, stomach, and head have layered tissue structures that require corresponding forward solutions.

Usually, the radiative transfer equation (RTE) is considered to be the most precise equation for describing light propagation in biological tissues and several solutions for the RTE have been reported in previous literature. However, the direct analytical solutions of the RTE cannot be easily obtained. Instead, the discrete ordinates (*S*
_{N}) and the spherical harmonics equations (*P*
_{N}) have been established to solve the RTE. The diffusion equation, which is the *P*
_{1} approximation of RTE, has been widely used to obtain the solution for light propagation in layered turbid media and several solutions have been reported in literature. Dayan used the Fourier and Laplace transforms to acquire the solutions for the two-layered diffusion equation,^{10} while Kienle acquired the solutions for the two-layered diffusion equation in the steady-state, frequency, and time domains using the Fourier transforms.^{11, 12} Martelli presented the solutions for the case of two and three layers using the Eigenfunction method and the perturbation model.^{13, 14, 15, 16} For the case of the *N*-layered diffusion equation, several groups have also reported their solutions.^{17, 18, 19, 20, 21, 22} Liemert presented the solutions for *N*-layered infinite or finite turbid media in the steady-state, frequency, and time domains using the two-dimensional (2-D) Fourier transform formalism^{19, 20} and also presented the solutions for a cylindrical geometry.^{21, 22} Liemert's group gave us a complete set of solutions for the *N*-layered diffusion equation, and they also studied the correctly and efficiently numerical 2-D inverse Fourier transform and provided the executable program on the Internet using the Delphi (Pascal) language.

However, the diffusion equation has several limitations. First, the diffusion equation requires the optical absorption coefficient μ_{a} to be much smaller than the scattering coefficient μ_{s}. A typical criterion for the applicable region of the diffusion equation is that the reduced albedo *a*′ must be larger than 0.98.^{23} Second, the diffusion equation is valid only when the radial distance is larger than 10 transport mean free paths.^{24} Therefore, the diffusion equation can be satisfied in the therapeutic window (λ = 650 − 950 *nm*) of most tissues with small absorption coefficients. For example, the absorption coefficient of gray and white matter in the adult head at a wavelength of 800 nm is 0.025 and 0.005 mm^{−1} and the reduced scattering coefficient is 2.5 and 6 mm^{−1}, respectively.^{25} The absorption coefficient of human skin *in vitro* in the therapeutic window is less than 0.05 mm^{−1} and the reduced scattering coefficient is almost 1.5 mm^{−1}.^{26} However, in the near-IR (λ > 1000 *nm*) region of the spectrum, the absorption coefficients of tissues become large and the reduced albedo can be 0.5 or smaller. For example, the absorption coefficient of the dermis of the skin when λ > 1400 *nm*, is larger than 0.5 mm^{−1} and the reduced scattering coefficient is less than 1.5 mm^{−1}.^{27} Evidently, there is a significant need for methods that accurately describe light propagation in media with large absorption coefficients, especially in the application of noninvasive blood glucose monitoring using NIR light in the region of 1000 to 2000 nm because of the high absorption of water and skin tissues.^{28} Hull and Foster^{29} derived Green's function of the steady-state RTE in the *P*
_{3} approximation for the infinite turbid media and demonstrated that the *P*
_{3} approximation models the radiance in highly absorbing media or in the region close to the source with an accuracy that is superior to that of the diffusion equation. In order to develop the diffusion equation and to simplify the *P*
_{3} approximation to be used for highly absorbing media, Hull and Foster investigated a modified diffusion approximation model called the hybrid diffusion-P3 equation.^{29} Klose also developed a method of the simplified spherical harmonics equations (*SP*
_{N}) to approximate the more complicated equation of radiative transfer for modeling light propagation in biological tissues and found that the *SP*
_{N} significantly improves the diffusion solution in transport-like domains with high absorption and small geometries.^{30}

In this study, we present the solution for the hybrid diffusion-P3 equation for *N*-layered turbid media by combining the derivation of the *N*-layered diffusion equation in Liemert and Kienle's paper^{20} and the modified diffusion approximation model investigated by Hull and Foster^{29} in order to find a solution for light propagation in *N*-layered turbid media for the case of large absorption coefficients. The hybrid diffusion-P3 equation is solved for an *N*-layered finite or infinite turbid medium in the steady-state domain for one point source using the extrapolated boundary condition. The Fourier transform formalism is applied to derive the analytical solutions of the fluence rate in the Fourier space. The 2-D inverse Fourier transform is numerically calculated. In addition, the solutions of hybrid diffusion-P3 equation are compared to those of the diffusion equation and Monte Carlo simulations.

## 2.

## Theory

In this section, we use the Fourier transform formalism to derive the steady-state hybrid diffusion-P3 equation for an *N*-layered turbid medium having finite or infinite extensions.

## 2.1.

### N-layered Diffusion Equation

First, the *N*-layered diffusion equation can be described as follows,^{20} when the turbid media in every layer is homogeneous:

## Eq. 1

[TeX:] \documentclass[12pt]{minimal}\begin{document}\begin{eqnarray} D_1 \Delta \Phi _1 (x,y,z) - \mu _{a1} \Phi _1 (x,y,z) = - q(x,y,z),\;0 \le z < l_1,\hspace*{-6pt}\nonumber\\ \end{eqnarray}\end{document} $$\begin{array}{c}\hfill {D}_{1}\Delta {\Phi}_{1}(x,y,z)-{\mu}_{a1}{\Phi}_{1}(x,y,z)=-q(x,y,z),\phantom{\rule{0.28em}{0ex}}0\le z<{l}_{1},\end{array}$$## Eq. 2

[TeX:] \documentclass[12pt]{minimal}\begin{document}\begin{eqnarray} D_k \Delta \Phi _k (x,y,z) - \mu _{ak} \Phi _k (x,y,z) &=& 0,\;\sum\limits_{j = 1}^{k - 1} {l_j } < z \le \sum\limits_{j = 1}^k {l_j },\nonumber\\ k &=& 2,3,\ldots,N, \end{eqnarray}\end{document} $$\begin{array}{ccc}\hfill {D}_{k}\Delta {\Phi}_{k}(x,y,z)-{\mu}_{ak}{\Phi}_{k}(x,y,z)& =& 0,\phantom{\rule{0.28em}{0ex}}\sum _{j=1}^{k-1}{l}_{j}<z\le \sum _{j=1}^{k}{l}_{j},\hfill \\ \hfill k& =& 2,3,...,N,\hfill \end{array}$$_{k}is the fluence rate of the layer

*k*,

*l*

_{k}is the thickness, and [TeX:] $\mu _{sk}^\prime $ ${\mu}_{sk}^{\prime}$ and μ

_{ak}are the reduced scattering and absorption coefficients, respectively. [TeX:] $D_k = 1/[3(\mu _{sk}^\prime + \mu _{ak})]$ ${D}_{k}=1/\left[3({\mu}_{sk}^{\prime}+{\mu}_{ak})\right]$ is the diffusion coefficient.

*q*is the source term.

For a semi-infinite situation, the derivation of the *N*-layered diffusion equation involves characterizing the source term and satisfying the appropriate boundary condition. For the light propagation in biological tissues, a pencil beam is usually modeled as an infinite line of isotropic sources on a semi-infinite scattering medium. According to the formalism developed by Farrell,^{31} we can obtain the distribution equation of the source

## Eq. 3

[TeX:] \documentclass[12pt]{minimal}\begin{document}\begin{equation} q(z) = \frac{a^{\prime}\mu_t^{\prime}}{4\pi }\exp (- \mu_t^{\prime}z), \end{equation}\end{document} $$q\left(z\right)=\frac{{a}^{\prime}{\mu}_{t}^{\prime}}{4\pi}\mathrm{exp}(-{\mu}_{t}^{\prime}z),$$In order to represent the pencil beam in terms of simpler source distributions, isotropic point sources are adopted to describe the distributions of the sources. The distributions of isotropic point sources have a dipole moment with respect to an origin at the air–tissue interface as the distribution in Eq. 3.^{29} To satisfy the dipole moment, a single point source is needed

## Eq. 4

[TeX:] \documentclass[12pt]{minimal}\begin{document}\begin{equation} \int_0^\infty {za'\mu _t^{\prime} } \exp (- \mu _t^{\prime} z)dz = \int_0^\infty {za'} \delta (z - z_0)dz. \end{equation}\end{document} $${\int}_{0}^{\infty}z{a}^{\prime}{\mu}_{t}^{\prime}\mathrm{exp}(-{\mu}_{t}^{\prime}z)dz={\int}_{0}^{\infty}z{a}^{\prime}\delta (z-{z}_{0})dz.$$*q*(

*x*,

*y*,

*z*)in Eq. 1 can be expressed as

## Eq. 5

[TeX:] \documentclass[12pt]{minimal}\begin{document}\begin{equation} q(x,y,z) = a_1 ^{\prime} \delta (x,y,z - z_0), \end{equation}\end{document} $$q(x,y,z)={a}_{1}^{\prime}\delta (x,y,z-{z}_{0}),$$To obtain the solutions of the radiance emitted from a semi-infinite medium, appropriate boundary conditions must be prescribed at the interface between the surrounding media and the biological tissue. The extrapolated boundary condition is one of the boundary conditions used for a semi-infinite scattering medium,^{32, 33} shown in Fig. 1, where Φ(*x*, *y*, *z* = −*z*
_{b}) = 0 and *z*
_{b} is the position of the extrapolated boundary.

We then used the Fourier transform approach to solve the *N*-layered diffusion equation based on Liemert and Kienle's paper.^{20} We apply the 2-D Fourier transform

## Eq. 6

[TeX:] \documentclass[12pt]{minimal}\begin{document}\begin{equation} \Phi _k (z,s_1,s_2) = \int_{ - \infty }^\infty {\int_{ - \infty }^\infty {\Phi _k (z,x,y)} } e^{i(s_1 x + s_2 y)} dxdy \end{equation}\end{document} $${\Phi}_{k}(z,{s}_{1},{s}_{2})={\int}_{-\infty}^{\infty}{\int}_{-\infty}^{\infty}{\Phi}_{k}(z,x,y){e}^{i({s}_{1}x+{s}_{2}y)}dxdy$$*N*-layered diffusion equation in Fourier space

## Eq. 7

[TeX:] \documentclass[12pt]{minimal}\begin{document}\begin{equation} \frac{{\partial ^2 }}{{\partial z^2 }}\Phi _1 (z,s) - a_1^2 \Phi _1 (z,s) = - \frac{1}{{D_1 }}a_1 ^{\prime} \delta (z - z_0),\quad 0 \le z < l_1, \end{equation}\end{document} $$\frac{{\partial}^{2}}{\partial {z}^{2}}{\Phi}_{1}(z,s)-{a}_{1}^{2}{\Phi}_{1}(z,s)=-\frac{1}{{D}_{1}}{a}_{1}^{\prime}\delta (z-{z}_{0}),\phantom{\rule{1em}{0ex}}0\le z<{l}_{1},$$## Eq. 8

[TeX:] \documentclass[12pt]{minimal}\begin{document}\begin{eqnarray} \frac{{\partial ^2 }}{{\partial z^2 }}\Phi _k (z,s) - a_k^2 \Phi _k (z,s) &=& 0,\quad \sum\limits_{j = 1}^{k - 1} {l_j } < z \le \sum\limits_{j = 1}^k {l_j },\nonumber\\ k &=& 2,3,\ldots,N, \end{eqnarray}\end{document} $$\begin{array}{ccc}\hfill \frac{{\partial}^{2}}{\partial {z}^{2}}{\Phi}_{k}(z,s)-{a}_{k}^{2}{\Phi}_{k}(z,s)& =& 0,\phantom{\rule{1em}{0ex}}\sum _{j=1}^{k-1}{l}_{j}<z\le \sum _{j=1}^{k}{l}_{j},\hfill \\ \hfill k& =& 2,3,...,N,\hfill \end{array}$$## 2.2.

### N-layered Hybrid Diffusion-P3 Equation

The *N*-layered hybrid diffusion-P3 equation is the method based on *N*-layered diffusion equation, combined with the coefficients of the *P*
_{3} approximation, which is the high-order approximation of the solution of RTE. It has been demonstrated that the *P*
_{3} approximation for modeling the radiance in highly absorbing media or in the region close to the source is more accurate than the diffusion equation.^{29} *v*
^{−} is the high-order coefficient of the P3-Green function, which is named as the asymptotic attenuation coefficient and approximately equal to the effective attenuation coefficient
[TeX:]
$\mu _{\hbox{\scriptsize\textit{eff}}} $
${\mu}_{\mathit{\text{eff}}}$
of the diffusion-Green function for large albedos and increasingly deviates from
[TeX:]
$\mu _{\hbox{\scriptsize\textit{eff}}} $
${\mu}_{\mathit{\text{eff}}}$
as the albedo decreases.^{29} In addition, with the use of
[TeX:]
$D = \mu _a /(\mu _{\hbox{\scriptsize\textit{eff}}})^2 $
$D={\mu}_{a}/{\left({\mu}_{\mathit{\text{eff}}}\right)}^{2}$
, it has also been demonstrated that the hybrid diffusion-P3 approximation for semi-infinite media is more accurate than diffusion approximation for the case of large absorption coefficients.^{34} So we can obtain that *v*
^{−} associates with
[TeX:]
$\mu _{\hbox{\scriptsize\textit{eff}}} $
${\mu}_{\mathit{\text{eff}}}$
of diffusion-Green function and
[TeX:]
${{\mu _{ak} } / {D_k }}$
${\mu}_{ak}/{D}_{k}$
and *D*
_{k} of diffusion equation associate with
[TeX:]
$(v_{^k }^ -)^2 $
${\left({v}_{k}^{-}\right)}^{2}$
and
[TeX:]
$D_{asym\_k} $
${D}_{asym\_k}$
of *P*
_{3} approximation, respectively.
[TeX:]
$(v_{^k }^ -)^2 $
${\left({v}_{k}^{-}\right)}^{2}$
and
[TeX:]
$D_{asym\_k} $
${D}_{asym\_k}$
are the high-order coefficients compared to
[TeX:]
${{\mu _{ak} } / {D_k }}$
${\mu}_{ak}/{D}_{k}$
and *D*
_{k} of the diffusion equation.

Thus, the *N*-layered hybrid diffusion-P3 equation is obtained based on the *N*-layered diffusion equation by replacing
[TeX:]
${{\mu _{ak} } / {D_k }}$
${\mu}_{ak}/{D}_{k}$
and *D*
_{k} with
[TeX:]
$(v_{^k }^ -)^2 $
${\left({v}_{k}^{-}\right)}^{2}$
and
[TeX:]
$D_{asym\_k} $
${D}_{asym\_k}$
, respectively.

## Eq. 9

[TeX:] \documentclass[12pt]{minimal}\begin{document}\begin{eqnarray} &\displaystyle v_{_k }^ - = \frac{1}{{\sqrt {18} }}\left(\beta _k - \sqrt {\beta _k ^2 - \gamma _{ak} }\right)^{\frac{1}{2}} \nonumber\\ &\displaystyle \beta _k \equiv 27\mu _{ak} (\mu _{ak} + \mu _{sk} ^\prime) + 28\mu _{ak} (\mu _{ak} + \mu _{sk} ^\prime \delta)\nonumber\\ &\displaystyle \hspace*{22pt}+\, 35\mu _{ak} (\mu _{ak} + \mu _{sk} ^\prime \gamma _k)\mu _{ak} (\mu _{ak} + \mu _{sk} ^\prime \delta _k),\\ &\displaystyle \gamma _{ak} \equiv 3780(\mu _{ak} + \mu _{sk} ^\prime)(\mu _{ak} + \mu _{sk} ^\prime \gamma _k)(\mu _{ak} + \mu _{sk} ^\prime \delta _k), \nonumber\\ &\displaystyle \gamma _k = (1 - g_{2k})/(1 - g_{1k}),\nonumber\\ &\displaystyle\delta _k = (1 - g_{3k})/(1 - g_{1k}), \nonumber\\ &\displaystyle D_{asym\_k} = \mu _{ak} /(v_{_k }^ -)^2, \nonumber \end{eqnarray}\end{document} $$\begin{array}{cc}& {\displaystyle {v}_{k}^{-}=\frac{1}{\sqrt{18}}{\left({\beta}_{k}-\sqrt{{\beta}_{k}^{2}-{\gamma}_{ak}}\right)}^{\frac{1}{2}}}\\ & {\displaystyle {\beta}_{k}\equiv 27{\mu}_{ak}({\mu}_{ak}+{\mu}_{sk}^{\prime})+28{\mu}_{ak}({\mu}_{ak}+{\mu}_{sk}^{\prime}\delta )}\\ & {\displaystyle \phantom{\rule{22.0pt}{0ex}}+\phantom{\rule{0.16em}{0ex}}35{\mu}_{ak}({\mu}_{ak}+{\mu}_{sk}^{\prime}{\gamma}_{k}){\mu}_{ak}({\mu}_{ak}+{\mu}_{sk}^{\prime}{\delta}_{k}),}\\ & {\displaystyle {\gamma}_{ak}\equiv 3780({\mu}_{ak}+{\mu}_{sk}^{\prime})({\mu}_{ak}+{\mu}_{sk}^{\prime}{\gamma}_{k})({\mu}_{ak}+{\mu}_{sk}^{\prime}{\delta}_{k}),}\\ & {\displaystyle {\gamma}_{k}=(1-{g}_{2k})/(1-{g}_{1k}),}\\ & {\displaystyle {\delta}_{k}=(1-{g}_{3k})/(1-{g}_{1k}),}\\ & {\displaystyle {D}_{asym\_k}={\mu}_{ak}/{\left({v}_{k}^{-}\right)}^{2},}\end{array}$$*g*

_{1k},

*g*

_{2k}, and

*g*

_{3k}are the first-, second-, and third-moments of the phase function in the layer

*k*, respectively. For the Henyey–Greenstein phase function, γ

_{k}= 1 +

*g*

_{1k}and [TeX:] $\delta _k = 1 + g_{1k}\break + g_{1k} ^2 $ ${\delta}_{k}=1+{g}_{1k}+{g}_{1k}^{2}$ .

Thus, we get the expression of the hybrid diffusion-P3 equation for *N*-layered turbid media in Fourier space

## Eq. 10

[TeX:] \documentclass[12pt]{minimal}\begin{document}\begin{eqnarray} \frac{{\partial ^2 }}{{\partial z^2 }}\Phi _{_1 }^H (z,s) - \big(a_1^H\big)^2 \Phi _{_1 }^H (z,s) &=& - \frac{1}{{D_{asym\_1} }}a_1 ^{\prime} \delta (z - z_0),\nonumber\\ && 0 \le z \le l_1, \end{eqnarray}\end{document} $$\begin{array}{ccc}\hfill \frac{{\partial}^{2}}{\partial {z}^{2}}{\Phi}_{1}^{H}(z,s)-{\left({a}_{1}^{H}\right)}^{2}{\Phi}_{1}^{H}(z,s)& =& -\frac{1}{{D}_{asym\_1}}{a}_{1}^{\prime}\delta (z-{z}_{0}),\hfill \\ & & 0\le z\le {l}_{1},\hfill \end{array}$$## Eq. 11

[TeX:] \documentclass[12pt]{minimal}\begin{document}\begin{eqnarray} \frac{{\partial ^2 }}{{\partial z^2 }}\Phi _{_k }^H (z,s) - \big(a_k^H\big)^2 \Phi _{_k }^H (z,s) &=& 0,\quad \sum\limits_{j = 1}^{k - 1} {l_j } \le z \le \sum\limits_{j = 1}^k {l_j },\nonumber\\ k &=& 2,3,\ldots,N, \end{eqnarray}\end{document} $$\begin{array}{ccc}\hfill \frac{{\partial}^{2}}{\partial {z}^{2}}{\Phi}_{k}^{H}(z,s)-{\left({a}_{k}^{H}\right)}^{2}{\Phi}_{k}^{H}(z,s)& =& 0,\phantom{\rule{1em}{0ex}}\sum _{j=1}^{k-1}{l}_{j}\le z\le \sum _{j=1}^{k}{l}_{j},\hfill \\ \hfill k& =& 2,3,...,N,\hfill \end{array}$$*k*th layer.

The following boundary conditions are used in the Fourier space for an finite *N*-layered turbid medium:^{20}

## Eq. 12

[TeX:] \documentclass[12pt]{minimal}\begin{document}\begin{equation} \Phi _1^H (- z_{b1},s) = 0, \end{equation}\end{document} $${\Phi}_{1}^{H}(-{z}_{b1},s)=0,$$## Eq. 13

[TeX:] \documentclass[12pt]{minimal}\begin{document}\begin{equation} \Phi_{1}^{{H}^{\prime}} (z_0,s) = \Phi_1^H (z_0,s),\;\;z_0 < l_1, \end{equation}\end{document} $${\Phi}_{1}^{{H}^{\prime}}({z}_{0},s)={\Phi}_{1}^{H}({z}_{0},s),\phantom{\rule{0.28em}{0ex}}\phantom{\rule{0.28em}{0ex}}{z}_{0}<{l}_{1},$$## Eq. 14

[TeX:] \documentclass[12pt]{minimal}\begin{document}\begin{equation} \left. {\frac{{\partial \Phi _1^H (z,s)}}{{\partial z}}} \right|_{z = z_0 } - \left. {\frac{{\partial \Phi_{1}^{{H}^{\prime}} (z,s)}}{{\partial z}}} \right|_{z = z_0 } = \frac{{a_1 ^{\prime} }}{{D_{asym\_1} }}, \end{equation}\end{document} $${\left.\frac{\partial {\Phi}_{1}^{H}(z,s)}{\partial z}\right|}_{z={z}_{0}}-{\left.\frac{\partial {\Phi}_{1}^{{H}^{\prime}}(z,s)}{\partial z}\right|}_{z={z}_{0}}=\frac{{a}_{1}^{\prime}}{{D}_{asym\_1}},$$## Eq. 15

[TeX:] \documentclass[12pt]{minimal}\begin{document}\begin{equation} \frac{{\Phi _k^H (L_k,s)}}{{\Phi _{k + 1}^H (L_k,s)}} = \left(\frac{{n_k }}{{n_{k + 1} }}\right)^2,\;L_k = \sum\limits_{j = 1}^k {l_j },\;1 \le k \le N - 1, \end{equation}\end{document} $$\frac{{\Phi}_{k}^{H}({L}_{k},s)}{{\Phi}_{k+1}^{H}({L}_{k},s)}={\left(\frac{{n}_{k}}{{n}_{k+1}}\right)}^{2},\phantom{\rule{0.28em}{0ex}}{L}_{k}=\sum _{j=1}^{k}{l}_{j},\phantom{\rule{0.28em}{0ex}}1\le k\le N-1,$$## Eq. 16

[TeX:] \documentclass[12pt]{minimal}\begin{document}\begin{eqnarray} \left. {D_{asym\_k} \frac{{\partial \Phi _k^H (z,s)}}{{\partial z}}} \right|_{z = L_k } &=& \left. {D_{asym\_k + 1} \frac{{\partial \Phi _{k + 1}^H (z,s)}}{{\partial z}}} \right|_{z = L_k },\nonumber\\ L_k &=& \sum\limits_{j = 1}^k {l_j },\quad 1 \le k \le N - 1, \end{eqnarray}\end{document} $$\begin{array}{ccc}\hfill {\left.{D}_{asym\_k}\frac{\partial {\Phi}_{k}^{H}(z,s)}{\partial z}\right|}_{z={L}_{k}}& =& {\left.{D}_{asym\_k+1}\frac{\partial {\Phi}_{k+1}^{H}(z,s)}{\partial z}\right|}_{z={L}_{k}},\hfill \\ \hfill {L}_{k}& =& \sum _{j=1}^{k}{l}_{j},\phantom{\rule{1em}{0ex}}1\le k\le N-1,\hfill \end{array}$$## Eq. 17

[TeX:] \documentclass[12pt]{minimal}\begin{document}\begin{equation} \Phi _N^H (L_N + z_{b2},s) = 0, \end{equation}\end{document} $${\Phi}_{N}^{H}({L}_{N}+{z}_{b2},s)=0,$$*z*

_{b1}and

*z*

_{b2}are the positions of the extrapolated boundaries above the first layer and below the last layer, respectively.

## Eq. 18

[TeX:] \documentclass[12pt]{minimal}\begin{document}\begin{equation} z_{b1} = \frac{{1 + R_{{\hbox{\scriptsize\textit{eff}}}\_1} }}{{1 - R_{{\hbox{\scriptsize\textit{eff}}}\_1} }}2D_{asym\_1},\;z_{b2} = \frac{{1 + R_{{\hbox{\scriptsize\textit{eff}}}\_N} }}{{1 - R_{{\hbox{\scriptsize\textit{eff}}}\_N} }}2D_{asym\_N}, \end{equation}\end{document} $${z}_{b1}=\frac{1+{R}_{\mathit{\text{eff}}\_1}}{1-{R}_{\mathit{\text{eff}}\_1}}2{D}_{asym\_1},\phantom{\rule{0.28em}{0ex}}{z}_{b2}=\frac{1+{R}_{\mathit{\text{eff}}\_N}}{1-{R}_{\mathit{\text{eff}}\_N}}2{D}_{asym\_N},$$^{33}

## 2.3.

### Solution of N-layered Hybrid Diffusion-P3 Equation

The solution of the *N*-layered hybrid diffusion-P3 equation [Eqs.
10, 11, 12, 13, 14, 15, 16, 17] for the fluence rate in the Fourier space is solved based on the derivation results from Liemert and Kienle's paper and by applying Cramer's rule.^{20} The solution for the fluence rate in the first layer (above the single point source) is given by

## Eq. 19

[TeX:] \documentclass[12pt]{minimal}\begin{document}\begin{eqnarray} \Phi _1^H (z,s) &=& a_1 ^{\prime} \left\{ \frac{{a_1^H \beta _3 D_{asym\_1} \cosh \big[a_1^H (l_1 - z)\big] + a_2^H \gamma _3 D_{asym\_2} \big(n_2^2 /n_1^2\big)\sinh \big[a_1^H (l_1 - z)\big]}}{{a_1^H \beta _3 D_{asym\_1} \cosh \big[a_1^H (l_1 + z_{b1})\big] + a_2^H \gamma _3 D_{asym\_2} \big(n_2^2 /n_1^2\big)\sinh \big[a_1^H (l_1 + z_{b1})\big]}}\right. \nonumber\\ && \left.\times\, \frac{{\sinh \big[a_1^H (z_0 + z_{b1})\big]}}{{a_1^H D_{asym\_1} }} - \frac{{\sinh \big[a_1^H (z_0 - z)\big]}}{{a_1^H D_{asym\_1} }}\right\}, \end{eqnarray}\end{document} $$\begin{array}{ccc}\hfill {\Phi}_{1}^{H}(z,s)& =& {a}_{1}^{\prime}\left\{\frac{{a}_{1}^{H}{\beta}_{3}{D}_{asym\_1}\mathrm{cosh}\left[{a}_{1}^{H}({l}_{1}-z)\right]+{a}_{2}^{H}{\gamma}_{3}{D}_{asym\_2}\left({n}_{2}^{2}/{n}_{1}^{2}\right)\mathrm{sinh}\left[{a}_{1}^{H}({l}_{1}-z)\right]}{{a}_{1}^{H}{\beta}_{3}{D}_{asym\_1}\mathrm{cosh}\left[{a}_{1}^{H}({l}_{1}+{z}_{b1})\right]+{a}_{2}^{H}{\gamma}_{3}{D}_{asym\_2}\left({n}_{2}^{2}/{n}_{1}^{2}\right)\mathrm{sinh}\left[{a}_{1}^{H}({l}_{1}+{z}_{b1})\right]}\right.\hfill \\ & & \left.\times \phantom{\rule{0.16em}{0ex}}\frac{\mathrm{sinh}\left[{a}_{1}^{H}({z}_{0}+{z}_{b1})\right]}{{a}_{1}^{H}{D}_{asym\_1}}-\frac{\mathrm{sinh}\left[{a}_{1}^{H}({z}_{0}-z)\right]}{{a}_{1}^{H}{D}_{asym\_1}}\right\},\hfill \end{array}$$*N*th layer (

*L*

_{N − 1}⩽

*z*⩽

*L*

_{N}) can be given by

## Eq. 20

[TeX:] \documentclass[12pt]{minimal}\begin{document}\begin{equation} \Phi _N^H (z,s) = a_1 ^{\prime} \left\{ \frac{{\displaystyle\prod\nolimits_{k = 2}^{N - 1} {a_k^H D_{asym\_k} (n_N /n_1)^2 \sinh \big[a_1^H (z_0 + z_{b1})\big]\sinh \big[a_N^H (L_N + z_{b2} - z)\big]} }}{{a_1^H \beta _3 D_{asym\_1} \cosh \big[a_1^H (l_1 + z_{b1})\big] + a_2^H \gamma _3 D_{asym\_2} \big(n_2^2 /n_1^2\big)\sinh \big[a_1^H (l_1 + z_{b1})\big]}}\right\}. \end{equation}\end{document} $${\Phi}_{N}^{H}(z,s)={a}_{1}^{\prime}\left\{\frac{{\displaystyle {\prod}_{k=2}^{N-1}{a}_{k}^{H}{D}_{asym\_k}{({n}_{N}/{n}_{1})}^{2}\mathrm{sinh}\left[{a}_{1}^{H}({z}_{0}+{z}_{b1})\right]\mathrm{sinh}\left[{a}_{N}^{H}({L}_{N}+{z}_{b2}-z)\right]}}{{a}_{1}^{H}{\beta}_{3}{D}_{asym\_1}\mathrm{cosh}\left[{a}_{1}^{H}({l}_{1}+{z}_{b1})\right]+{a}_{2}^{H}{\gamma}_{3}{D}_{asym\_2}\left({n}_{2}^{2}/{n}_{1}^{2}\right)\mathrm{sinh}\left[{a}_{1}^{H}({l}_{1}+{z}_{b1})\right]}\right\}.$$*N*⩾ 3 and finite turbid media, the quantities β

_{3}and γ

_{3}are obtained by using recursion formulas. The values of β

_{3}and γ

_{3}for finite

*N*-layered turbid media are initiated by

## Eq. 21

[TeX:] \documentclass[12pt]{minimal}\begin{document}\begin{eqnarray} \begin{array}{l} \left({\begin{array}{c} {\beta _N } \\ {\gamma _N } \\ \end{array}} \right) = a_{N - 1}^H D_{asym\_N - 1}^H \sinh \big[a_N^H (l_N + z_{b2})\big]\left[ {\begin{array}{c} {\cosh \big(a_{N - 1}^H l_{N - 1}\big)} \\[3pt] {\sinh \big(a_{N - 1}^H l_{N - 1}\big)} \\[3pt] \end{array}} \right] \\[15pt] \qquad\qquad\; +\, a_N^H D_{asym\_N}^H \left(\displaystyle\frac{{n_N^2 }}{{n_{N - 1}^2 }}\right)\cosh \big[a_N^H (l_N + z_{b2})\big]\left[ {\begin{array}{c} {\cosh \big(a_{N - 1}^H l_{N - 1}\big)} \\[3pt] {\sinh \big(a_{N - 1}^H l_{N - 1}\big)} \\[3pt] \end{array}} \right], \\ \end{array} \end{eqnarray}\end{document} $$\begin{array}{c}\hfill \begin{array}{c}\left(\begin{array}{c}{\beta}_{N}\\ {\gamma}_{N}\end{array}\right)={a}_{N-1}^{H}{D}_{asym\_N-1}^{H}\mathrm{sinh}\left[{a}_{N}^{H}({l}_{N}+{z}_{b2})\right]\left[\begin{array}{c}\mathrm{cosh}\left({a}_{N-1}^{H}{l}_{N-1}\right)\\ \mathrm{sinh}\left({a}_{N-1}^{H}{l}_{N-1}\right)\end{array}\right]\hfill \\ \phantom{\rule{2.em}{0ex}}\phantom{\rule{2.em}{0ex}}\phantom{\rule{0.28em}{0ex}}+\phantom{\rule{0.16em}{0ex}}{a}_{N}^{H}{D}_{asym\_N}^{H}\left({\displaystyle \frac{{n}_{N}^{2}}{{n}_{N-1}^{2}}}\right)\mathrm{cosh}\left[{a}_{N}^{H}({l}_{N}+{z}_{b2})\right]\left[\begin{array}{c}\mathrm{cosh}\left({a}_{N-1}^{H}{l}_{N-1}\right)\\ \mathrm{sinh}\left({a}_{N-1}^{H}{l}_{N-1}\right)\end{array}\right],\hfill \end{array}\end{array}$$## Eq. 22

[TeX:] \documentclass[12pt]{minimal}\begin{document}\begin{equation} \left({\begin{array}{*{20}c} {\beta _{k - 1} } \\ {\gamma _{k - 1} } \\ \end{array}} \right) = \left[ {\begin{array}{*{20}c} {\begin{array}{*{20}c} {a_{k - 2}^H D_{asym\_k - 2}^H \cosh \big(a_{k - 2}^H l_{k - 2}\big)} \\[5pt] {a_{k - 2}^H D_{asym\_k - 2}^H \sinh \big(a_{k - 2}^H l_{k - 2}\big)} \\[5pt] \end{array}} & {\begin{array}{*{20}c} {a_{k - 1}^H D_{asym\_k - 1}^H \big(n_{k - 1}^2 /n_{k - 2}^2\big)\sinh \big(a_{k - 2}^H l_{k - 2}\big)} \\[3pt] {a_{k - 1}^H D_{asym\_k - 1}^H \big(n_{k - 1}^2 /n_{k - 2}^2\big)\cosh \big(a_{k - 2}^H l_{k - 2}\big)} \\[3pt] \end{array}} \\ \end{array}} \right]\left({\begin{array}{*{20}c} {\beta _k } \\ {\gamma _k } \\ \end{array}} \right). \end{equation}\end{document} $$\left(\begin{array}{c}{\beta}_{k-1}\\ {\gamma}_{k-1}\end{array}\right)=\left[\begin{array}{cc}\begin{array}{c}{a}_{k-2}^{H}{D}_{asym\_k-2}^{H}\mathrm{cosh}\left({a}_{k-2}^{H}{l}_{k-2}\right)\\ {a}_{k-2}^{H}{D}_{asym\_k-2}^{H}\mathrm{sinh}\left({a}_{k-2}^{H}{l}_{k-2}\right)\end{array}& \begin{array}{c}{a}_{k-1}^{H}{D}_{asym\_k-1}^{H}\left({n}_{k-1}^{2}/{n}_{k-2}^{2}\right)\mathrm{sinh}\left({a}_{k-2}^{H}{l}_{k-2}\right)\\ {a}_{k-1}^{H}{D}_{asym\_k-1}^{H}\left({n}_{k-1}^{2}/{n}_{k-2}^{2}\right)\mathrm{cosh}\left({a}_{k-2}^{H}{l}_{k-2}\right)\end{array}\end{array}\right]\left(\begin{array}{c}{\beta}_{k}\\ {\gamma}_{k}\end{array}\right).$$*N*⩾ 3 and the semi-infinite

*N*-layered turbid medium (

*L*

_{N}→ ∞), the start term for calculating β

_{3}and γ

_{3}becomes

## Eq. 23

[TeX:] \documentclass[12pt]{minimal}\begin{document}\begin{eqnarray} \left({\begin{array}{*{20}c} {\beta _N } \\ {\gamma _N } \\ \end{array}} \right) &=& a_{N - 1}^H D_{asym\_N - 1}^H \left[ {\begin{array}{*{20}c} {\cosh \big(a_{N - 1}^H l_{N - 1}\big)} \\[3pt] {\sinh \big(a_{N - 1}^H l_{N - 1}\big)} \\[3pt] \end{array}} \right]\nonumber\\ && +\, a_N^H D_{asym\_N}^H \left(\frac{{n_N^2 }}{{n_{N - 1}^2 }}\right)\left[ {\begin{array}{*{20}c} {\cosh \big(a_{N - 1}^H l_{N - 1}\big)} \\[3pt] {\sinh \big(a_{N - 1}^H l_{N - 1}\big)} \\[3pt] \end{array}} \right],\quad \end{eqnarray}\end{document} $$\begin{array}{ccc}\hfill \left(\begin{array}{c}{\beta}_{N}\\ {\gamma}_{N}\end{array}\right)& =& {a}_{N-1}^{H}{D}_{asym\_N-1}^{H}\left[\begin{array}{c}\mathrm{cosh}\left({a}_{N-1}^{H}{l}_{N-1}\right)\\ \mathrm{sinh}\left({a}_{N-1}^{H}{l}_{N-1}\right)\end{array}\right]\hfill \\ & & +\phantom{\rule{0.16em}{0ex}}{a}_{N}^{H}{D}_{asym\_N}^{H}\left(\frac{{n}_{N}^{2}}{{n}_{N-1}^{2}}\right)\left[\begin{array}{c}\mathrm{cosh}\left({a}_{N-1}^{H}{l}_{N-1}\right)\\ \mathrm{sinh}\left({a}_{N-1}^{H}{l}_{N-1}\right)\end{array}\right],\phantom{\rule{1em}{0ex}}\hfill \end{array}$$In the case of the finite two-layered turbid medium (*N* = 2), the values of β_{3} and γ_{3} are

## Eq. 24

[TeX:] \documentclass[12pt]{minimal}\begin{document}\begin{equation} \left({\begin{array}{*{20}c} {\beta _3 } \\ {\gamma _3 } \\ \end{array}} \right) = \left({\begin{array}{*{20}c} {\sinh \big[a_2^H (l_2 + z_{b2})\big]} \\[3pt] {\cosh \big[a_2^H (l_2 + z_{b2})\big]} \\[3pt] \end{array}} \right). \end{equation}\end{document} $$\left(\begin{array}{c}{\beta}_{3}\\ {\gamma}_{3}\end{array}\right)=\left(\begin{array}{c}\mathrm{sinh}\left[{a}_{2}^{H}({l}_{2}+{z}_{b2})\right]\\ \mathrm{cosh}\left[{a}_{2}^{H}({l}_{2}+{z}_{b2})\right]\end{array}\right).$$*N*-layered turbid medium (

*L*

_{N}→ ∞), β

_{3}= γ

_{3}= 1.

Next, we used the 2-D inverse Fourier transform to obtain the fluence rate Φ(*x*, *y*, *z*) in the real space. In this step, the numerical calculation is performed to obtain the solutions for the fluence rate in real space due to the difficulties of getting the analytical solutions. The expression for the 2-D inverse Fourier transform is given by

## Eq. 25

[TeX:] \documentclass[12pt]{minimal}\begin{document}\begin{eqnarray} \Phi _k^H (x,y,z) &=& \frac{1}{{(2\pi)^2 }}\int_{ - \infty }^\infty \int_{ - \infty }^\infty \Phi _k^H (z,s)\nonumber\\ &&\times\exp [ - i(s_1 x + s_2 y)]{\mathop{\rm d}\nolimits} s_1 {\mathop{\rm d}\nolimits} s_2. \end{eqnarray}\end{document} $$\begin{array}{ccc}\hfill {\Phi}_{k}^{H}(x,y,z)& =& \frac{1}{{\left(2\pi \right)}^{2}}{\int}_{-\infty}^{\infty}{\int}_{-\infty}^{\infty}{\Phi}_{k}^{H}(z,s)\hfill \\ & & \times \mathrm{exp}[-i({s}_{1}x+{s}_{2}y)]\mathrm{d}{s}_{1}\mathrm{d}{s}_{2}.\hfill \end{array}$$## Eq. 26

[TeX:] \documentclass[12pt]{minimal}\begin{document}\begin{equation} \begin{array}{l} \Phi _k^H (x,y,z) = \displaystyle\frac{{\Delta s_1 \Delta s_2 }}{{(2\pi)^2 }}\sum\limits_{m = - \infty }^\infty {\sum\limits_{n = - \infty }^\infty {\Phi _k^H (z,m\Delta s_1,n\Delta s_2)} } \\ \quad \quad \quad \quad \quad \; \times \exp [ - i(m\Delta s_1 x + n\Delta s_2 y)]. \\ \end{array} \end{equation}\end{document} $$\begin{array}{c}{\displaystyle {\Phi}_{k}^{H}(x,y,z)=\frac{\Delta {s}_{1}\Delta {s}_{2}}{{\left(2\pi \right)}^{2}}\sum _{m=-\infty}^{\infty}\sum _{n=-\infty}^{\infty}{\Phi}_{k}^{H}(z,m\Delta {s}_{1},n\Delta {s}_{2})}\hfill \\ \phantom{\rule{1em}{0ex}}\phantom{\rule{1em}{0ex}}\phantom{\rule{1em}{0ex}}\phantom{\rule{1em}{0ex}}\phantom{\rule{1em}{0ex}}\phantom{\rule{0.28em}{0ex}}\times \mathrm{exp}[-i(m\Delta {s}_{1}x+n\Delta {s}_{2}y)].\hfill \end{array}$$^{20}

## Eq. 27

[TeX:] \documentclass[12pt]{minimal}\begin{document}\begin{eqnarray} \Phi _k^H (x,y,z) &=& \frac{{\Delta s_1 \Delta s_2 }}{{(2\pi)^2 }}\sum\limits_{m = - \infty }^\infty {\sum\limits_{n = - \infty }^\infty {\Phi _k^H (z,m\Delta s_1,n\Delta s_2)} } \nonumber\\ && \times\, \cos (m\Delta s_1 x)\cos (n\Delta s_2 y). \end{eqnarray}\end{document} $$\begin{array}{ccc}\hfill {\Phi}_{k}^{H}(x,y,z)& =& \frac{\Delta {s}_{1}\Delta {s}_{2}}{{\left(2\pi \right)}^{2}}\sum _{m=-\infty}^{\infty}\sum _{n=-\infty}^{\infty}{\Phi}_{k}^{H}(z,m\Delta {s}_{1},n\Delta {s}_{2})\hfill \\ & & \times \phantom{\rule{0.16em}{0ex}}\mathrm{cos}\left(m\Delta {s}_{1}x\right)\mathrm{cos}\left(n\Delta {s}_{2}y\right).\hfill \end{array}$$*y*= 0), and assuming the same sampling rate in both directions (Δ

*s*

_{1}= Δ

*s*

_{2}= Δ

*s*), it is sufficient to get the simpler expression of Eq. 27

## Eq. 28

[TeX:] \documentclass[12pt]{minimal}\begin{document}\begin{eqnarray} \Phi _k^H (x,0,z) &=& \left(\frac{{\Delta s}}{{2\pi }}\right)^2 \left\{ \Phi _k^H (z,0,0)\right.\nonumber\\ &&+\, 2\sum\limits_{m = 1}^\infty {\Phi _k^H (z,m\Delta s,0)} [1 + \cos (m\Delta sx)]\nonumber\\ && +\left. 4\sum\limits_{m = 1}^\infty \sum\limits_{n = 1}^\infty \Phi _k^H (z,m\Delta s,n\Delta s) \cos (m\Delta sx)}.\nonumber\\ \end{eqnarray}\end{document} $$\begin{array}{ccc}\hfill {\Phi}_{k}^{H}(x,0,z)& =& {\left(\frac{\Delta s}{2\pi}\right)}^{2}\left\{{\Phi}_{k}^{H}(z,0,0)\right.\hfill \\ & & +\phantom{\rule{0.16em}{0ex}}2\sum _{m=1}^{\infty}{\Phi}_{k}^{H}(z,m\Delta s,0)[1+\mathrm{cos}\left(m\Delta sx\right)]\hfill \\ & & +\left.4\sum _{m=1}^{\infty}\sum _{n=1}^{\infty}{\Phi}_{k}^{H}(z,m\Delta s,n\Delta s)\mathrm{cos}\left(m\Delta sx\right)\right\}.\hfill \end{array}$$In the second method (IFT2), the one-dimensional inverse Hankel transform is obtained from Eq. 25

## Eq. 29

[TeX:] \documentclass[12pt]{minimal}\begin{document}\begin{equation} \Phi _k^H (\rho,z) = \frac{1}{{2\pi }}\int_0^\infty {\Phi _k^H (z,s)sJ_0 (s\rho)ds}, \end{equation}\end{document} $${\Phi}_{k}^{H}(\rho ,z)=\frac{1}{2\pi}{\int}_{0}^{\infty}{\Phi}_{k}^{H}(z,s)s{J}_{0}\left(s\rho \right)ds,$$*J*

_{0}is the Bessel function of the first kind and zero order. The discrete integral method is used to calculate the fluence rate. We typically use 2 × 10

^{6}terms in the infinite sums to obtain high accuracy. The upper limit of the integral is set to [TeX:] $30\mu _s^\prime $ $30{\mu}_{s}^{\prime}$ . Additionally, the adaptive Gauss-Kronrod quadrature in the MATLAB Mathematics Toolbox

^{35}can be used to do numerical integration.

Using one of the two methods, the fluence rate in real space is obtained. Thereafter, the spatially resolved reflectance *R*
^{H}(ρ) of the first layer is given by^{32}

## Eq. 30

[TeX:] \documentclass[12pt]{minimal}\begin{document}\begin{eqnarray} R^H (\rho) &=& \int_{2\pi } {d\Omega [1 - R_{fres} (\theta)]} \frac{1}{{4\pi }}\big[\Phi _1^H (\rho,z = 0) \nonumber\\ && +\, 3D_{asym\_1}^H \frac{{\partial \Phi _1^H (\rho,z)}}{{\partial z}}|_{z = 0} \cos \theta\big]\cos \theta, \end{eqnarray}\end{document} $$\begin{array}{ccc}\hfill {R}^{H}\left(\rho \right)& =& {\int}_{2\pi}d\Omega [1-{R}_{fres}\left(\theta \right)]\frac{1}{4\pi}[{\Phi}_{1}^{H}(\rho ,z=0)\hfill \\ & & +\phantom{\rule{0.16em}{0ex}}3{D}_{asym\_1}^{H}\frac{\partial {\Phi}_{1}^{H}(\rho ,z)}{\partial z}{|}_{z=0}\mathrm{cos}\theta ]\mathrm{cos}\theta ,\hfill \end{array}$$*R*

_{fres}(θ) is the Fresnel reflection coefficient for a photon with an incident angle θ relative to the boundary.

## 3.

## Results

In this section, we first compare the different methods for computing the inverse Fourier transform. Then, we compare the solutions of the steady-state hybrid diffusion-P3 equation and diffusion equation for a semi-infinite *N*-layered turbid medium with the Monte Carlo simulation. The principles of the Monte Carlo simulation of photon transport have been thoroughly described.^{36} We use the MCML program to execute the Monte Carlo simulation developed by Lihong Wang.^{36} A pencil photon beam is normally incident upon the semi-infinite turbid medium. The Henyey–Greenstein phase function is assumed for the calculation of the scattering angle. The spatial resolution of the steady-state Monte Carlo simulations is 0.1 mm and the number of initial photons is 10^{8}.

## 3.1.

### Comparison of the Different Methods for Computing the Inverse Fourier Transform

The spatially resolved reflectance from *N*-layered turbid media was calculated using the two methods (IFT1 and IFT2). Figure 2 compares the reflectance from a semi-infinite three-layered turbid medium calculated with IFT1 and IFT2. It can be seen that the differences between the two methods are less than 10^{−3} using MATLAB language (15 to 16 significant digits; 8 bytes). The results are ascribed to the fact that the Fourier and Henkel transforms are mathematically equivalent in circumstances of circular symmetry.

Although these two methods are the mathematically simple methods to compute the inverse Fourier transform, they require lengthy calculation times. If a short calculation time is preferred, one can use the C or the Pascal language and the existing algorithms for a fast Fourier transform, which results in calculation times less than 10 ms.^{20}

## 3.2.

### Results for the Case of Small Absorption Coefficients

Figure 3 shows a comparison of the semi-infinite three-layered diffusion equation and hybrid diffusion-P3 equation with the Monte Carlo simulations in the case of low absorption coefficients. The absorption coefficients of the three layers are μ_{a1} = 0.01 mm^{−1}, μ_{a2} = 0.1 mm^{−1}, and μ_{a3} = 0.001 mm^{−1}, respectively. The reduced scattering coefficients of the three layers are μ_{s1}
^{′} = 1.2 mm^{−1}, μ_{s2}
^{′} = 1.1 mm^{−1}, and μ_{s3}
^{′} = 1.3 mm^{−1}, respectively. The thicknesses of the three layers are *l*
_{1} = 5 mm, *l*
_{2} = 5 mm, and *l*
_{3} = ∞, respectively. A good agreement among the diffusion equation, the hybrid diffusion-P3 equation, and the Monte Carlo simulations can be observed. Figure 4 gives the relative differences of the diffusion equation and the hybrid diffusion-P3 equation to the Monte Carlo simulations. At the radial distances of 0.5 mm < ρ < 12 mm, the differences are smaller than 5% with the noises of the Monte Carlo simulation and the numerical calculation accounted for. For large radial distances (ρ > 12 mm), the differences are smaller than 10% with the larger noises.

We can observe from Figs. 3 and 4 that the solutions between the diffusion equation and the hybrid diffusion-P3 equation are nearly the same when the absorption coefficient is small. As can be seen in Sec. 2.2, the hybrid diffusion-P3 equation is obtained by replacing *D*
_{k} with
[TeX:]
$D_{asym\_k} $
${D}_{asym\_k}$
. The main difference between the diffusion equation and the hybrid diffusion-P3 equation is the diffusion coefficient *D* of the diffusion equation and the asymptotic diffusion coefficient *D*
_{asym} of the hybrid diffusion-P3 equation. Figure 5 gives the comparison between the diffusion coefficient *D* and the asymptotic diffusion coefficient *D*
_{asym} in the range of μ_{a}. It can be observed from the results that as the absorption coefficient decreases, the values of the two coefficients (*D*
_{k} and
[TeX:]
$D_{asym\_k} $
${D}_{asym\_k}$
) become identical. Thus, it can be concluded that as the absorption coefficient decreases, the solutions for the *N*-layered diffusion equation and hybrid diffusion-P3 equation are almost equivalent. The conclusion of this section accords with that from Liemert and Kienle's paper.^{20}

## 3.3.

### Results for the Case of Large Absorption Coefficients

Figure 6 shows a comparison of the semi-infinite three-layered diffusion equation and hybrid diffusion-P3 equation with the Monte Carlo simulations in the case of large absorption coefficients. The absorption coefficients of the three layers are μ_{a1} = 1 mm^{−1}, μ_{a2} = 0.5 mm^{−1}, and μ_{a3} = 0.001 mm^{−1}, respectively. The reduced scattering coefficients of the three layers are μ_{s1}
^{′} = 1.2 mm^{−1}, μ_{s2}
^{′} = 1.1 mm^{−1}, and μ_{s3}
^{′} = 1.3 mm^{−1}, respectively. The thicknesses of the three layers are *l*
_{1} = 1 mm, *l*
_{2} = 5 mm, and *l*
_{3} = ∞. The results in Fig. 6 show that in the case of large absorption coefficients, the hybrid diffusion-P3 equation agrees well with the Monte Carlo simulation results, while *R*(ρ) calculated with the *N*-layered diffusion equation is not close to that of the Monte Carlo simulations. Figure 7 gives the relative differences between the *N*-layered diffusion equation and the *N*-layered hybrid diffusion equation to the Monte Carlo simulations. At small distances (ρ < 2 mm), both the diffusion equation and the hybrid diffusion-P3 equation greatly deviate from the Monte Carlo simulations. At the radial distances of 2 mm < ρ < 4.7 mm, the differences between the hybrid diffusion-P3 equation and the Monte Carlo simulations are less than 20%. For the large distances of 4.7 mm < ρ < 6 mm, the differences are less than 8%. However, the differences between the diffusion equation and the Monte Carlo simulations are larger than 20%. The results suggest that the diffusion equation is not valid for the case of large absorption coefficients. From the results in Figs. 6 and 7, it can be concluded that the model of the hybrid diffusion-P3 equation is more precise than the diffusion equation for the case of a medium with a large absorption coefficient.

Figure 8 shows another example of the semi-infinite five-layered turbid media for large absorption coefficients. The absorption coefficients of the five layers are μ_{a1} = 0.5 mm^{−1}, μ_{a2} = 0.2 mm^{−1}, μ_{a3} = 0.1 mm^{−1}, μ_{a4} = 0.3 mm^{−1}, and μ_{a5} = 0.001 mm^{−1}, respectively. The thicknesses of the five layers are *l*
_{1} = 1 mm, *l*
_{2} = 1 mm, *l*
_{3} = 1 mm, *l*
_{4} = 2 mm, and *l*
_{5} = ∞, respectively. The results show that the hybrid diffusion-P3 equation is more similar to the Monte Carlo simulations than the diffusion equation. Figure 9 gives the relative differences between the *N*-layered diffusion equation and the *N*-layered hybrid diffusion-P3 equation to the Monte Carlo simulations. At the radial distances of 1.6 mm < ρ < 2.5 mm, the differences between the hybrid diffusion-P3 equation and the Monte Carlo simulations is from 0 to 10%, and for the large distances of 2.5 mm < ρ < 6 mm, the differences between the hybrid diffusion-P3 equation and the Monte Carlo simulations are less than 20%. At small distances (ρ < 1.6 mm), the differences are largely affected by the light source. However, the differences between the diffusion equation and the Monte Carlo simulations are 0 to 45% for radial distances greater than 1.6 mm.

From Figs. 6 and 8, it can be seen that the hybrid diffusion-P3 equation is more precise than the diffusion equation for modeling light propagation in the *N*-layered semi-infinite media with large absorption coefficients. In order to investigate the versatility of this conclusion, the medium with different optical properties and geometries have been investigated. Concerning the article length, only a few results are illustrated in Figs.10, 11, 12. The results also demonstrate the accuracy of the hybrid diffusion-P3 equation for modeling light propagation in the *N*-layered semi-infinite media with large absorption coefficients.

## 4.

## Discussion

This paper discusses light propagation in *N*-layered turbid media. The solution of the hybrid diffusion-P3 equation is derived for *N*-layered finite or infinite turbid media. The solution is calculated in the steady-state domain for one point source using the extrapolated boundary condition. The Fourier transform formalism is applied to derive the analytical solutions of the fluence rate in the Fourier space based on the derivation of the *N*-layered diffusion equation from Liemert and Kienle's paper. The numerical calculation is performed to obtain the solutions for the fluence rate in real space because of the difficulty of obtaining analytical solutions. Two inverse Fourier transform methods are developed to calculate the fluence rate in real space.

In addition, the solution of the hybrid diffusion-P3 equation having an infinite or finite thick *N*th layer is compared to that of the diffusion equation and Monte Carlo simulations. The main difference between the diffusion equation and the hybrid diffusion-P3 equation is the diffusion coefficient *D* of the diffusion equation and the asymptotic diffusion coefficient *D*
_{asym} of the hybrid diffusion-P3 equation. The values of two coefficients (*D*
_{k} and
[TeX:]
$D_{asym\_k} $
${D}_{asym\_k}$
) become identical as the absorption coefficient decreases, thus it can be concluded that as the absorption coefficient decreases, the solutions of the *N*-layered diffusion equation and hybrid diffusion-P3 equation are almost equivalent. Simulation results show that, in the case of small absorption coefficients, the solutions of the *N*-layered diffusion equation and hybrid diffusion-P3 equation are almost equivalent and are in agreement with the Monte Carlo simulations. Additionally, we discussed the situation of the semi-infinite three-layered, five-layered, and seven-layered turbid media for the case of large absorption coefficients. It can be observed that the model of the hybrid diffusion-P3 equation is closer to the Monte Carlo simulation than the diffusion equation. Finally, we can make the conclusion that the model of the hybrid diffusion-P3 equation can replace the diffusion equation for light propagation in the turbid media for a wide range of absorption coefficients. The hybrid diffusion-P3 equation also has greater potential applications in the field of biomedical photonics than the diffusion equation.

Note that the derived solutions must satisfy the condition of one point source in the first layer *l*
_{1} > *z*
_{0}. As a result, the model of the *N*-layered hybrid diffusion-P3 equation cannot be applied to situations in which the first layer is very thin. In the future, we will be committed to solve this problem and study other boundary conditions and source terms to obtain more accurate solutions at small radial distances (ρ < 2 mm).

## Acknowledgments

The authors acknowledge the fruitful discussion with Dr. Liemert from Institut für Lasertechnologien in der Medizin und Meßtechnik at Universitat Ulm regarding the solution of the diffusion equation for *N*-layered turbid media. The authors would like to thank Anna H. Xue from the University of Washington for her help with revisions. The authors also acknowledge the funding supports from the National Natural Science Foundation of China (30870657, 60938002), and the Tianjin Municipal Government of China (09JCZDJC18200).

## References

*N*-layered turbid media: steady-state domain," Journal of Biomedical Optics 16(10), 105002 (1 October 2011). https://doi.org/10.1117/1.3640810