*N*-layered turbid media having a finite or an infinitely thick

*N*th layer. We focus on time-resolved light propagation in both the frequency and time domains. Based on our results for the steady-state domain, solutions of the

*N*-layered diffusion equations in the frequency and time domains are obtained by applying the Fourier transform technique. Different methods for calculation of the inverse Fourier transform are studied to validate the solutions, showing relative differences typically smaller than 10

^{-6}. The solutions are compared to Monte Carlo simulations, revealing good agreement. Finally, by applying the Laplace and Fourier transforms we derive a fast (≈ 1 ms) and accurate analytical solution for the time domain reflectance from a two-layered turbid medium having equal reduced scattering coefficients and refractive indices in both layers.

## 1.

## Introduction

In the first of our two papers dealing with light diffusion in
$N$
-layered turbid media we studied the light propagation in the steady-state domain.^{1} In this second part, we concentrate on solutions of the time-resolved diffusion equation in both the time and frequency domains.

In general, measurements in the steady-state domain, for which the considered turbid medium is illuminated by a continuously emitting source and the remitted or transmitted light is detected steadily, are relatively inexpensive and simple to handle.^{2} Contrarily, time-resolved measurements are more complex and expensive, but they contain more information on the investigated tissues.^{3} Time-resolved measurements are performed in the time and frequency domains. In the time domain, a short laser pulse of the order of picoseconds is incident onto the turbid medium and the remitted or transmitted light is detected in a time-resolved manner.^{4, 5} In the frequency domain, the incident source is intensity modulated in the range of
$100\phantom{\rule{0.3em}{0ex}}\mathrm{MHz}$
and the demodulation and phase shift of the remitted or transmitted light is detected.^{6}

In the literature, experiments in all three measurement regimes were reported assuming that the investigated turbid media can be considered homogeneous or heterogeneous.^{7} In this paper, we consider layered geometries, a special case of heterogeneous turbid media. Many parts of the human body can be approximated as layered turbid media, e.g., the head, the arm, and the leg. For the mathematical description of light propagation in such large tissue volumes the diffusion equation is most often applied. The literature on the solution of the layered diffusion equation is reviewed in the accompanying paper.^{1}

In this paper, we present the derived solutions of the time-resolved reflectance and transmittance from an
$N$
-layered turbid medium having a finite thick
$N$
’th layer, and the time-resolved reflectance from an
$N$
-layered turbid medium having an infinitely thick
$N$
’th layer. The solutions in the frequency domain are obtained by using the Fourier transform technique, as reported in the accompanying paper.^{1} Different methods are applied for calculation of the inverse Fourier transform to compare their performance and accuracy. The solutions in the time domain are obtained by applying the FFT (fast Fourier transform) to the solutions in the frequency domain. Additionally, we derive an analytical formula for the reflectance from a two-layered turbid medium having an infinitely thick second layer for the case of equal reduced scattering coefficients and refractive indices. Finally, the obtained solutions are compared to Monte Carlo simulations.

## 2.

## Theory

In this section, the derivation of the solutions of the $N$ -layered diffusion equation in the frequency and time domains is presented. We start with the diffusion equation in the time domain:

## Eq. 1

$$(\frac{1}{c}\frac{\partial}{\partial t}+{\mu}_{a}-D\Delta )\Phi (\mathbf{r},t)=S(\mathbf{r},t),$$## Eq. 2

$$(\frac{1}{c}\frac{\partial}{\partial t}+{\mu}_{a}-D\Delta )G(\mathbf{r},t)=\delta \left(\mathbf{r}\right)\delta \left(t\right).$$For an arbitrary source in time and space $S(\mathbf{r},t)$ , the fluence rate can be calculated with the knowledge of Green’s function by applying

## Eq. 3

$$\Phi (\mathbf{r},t)={\int}_{-\infty}^{\infty}{\int}_{V}G({\mathbf{r}}^{\prime},{t}^{\prime})S(\mathbf{r}-{\mathbf{r}}^{\prime},t-{t}^{\prime}){\mathrm{d}}^{3}{r}^{\prime}\phantom{\rule{0.2em}{0ex}}\mathrm{d}{t}^{\prime}.$$We propose two different methods to calculate $G(\mathbf{r},t)$ . In the first method, Eq. 2 is solved using the Laplace and Fourier transform techniques. We first apply the Laplace transform with respect to the $t$ coordinate and the 2-D Fourier transform relative to the $x$ and $y$ coordinates to obtain an ordinary differential equation depending on the $z$ coordinate in the inverse Laplace and Fourier space. Then, the boundary conditions are applied relative to the $z$ coordinate. The next step is to solve the resulting equation system. Finally, the inverse 2-D Fourier and Laplace transforms are employed to obtain the time-resolved fluence rate in real space. Note that the Laplace transform with respect to the $t$ coordinate can, in principle, be substituted by an additional Fourier transform. Thus in this case also, the inverse Laplace transform must be replaced by an inverse Fourier transform.

The challenge of this method is that it is difficult to find analytical expressions for the inverse Laplace and Fourier transforms.^{8} However, we succeed in deriving a fast analytical solution for a two-layered turbid medium in the time domain having an infinitely thick bottom layer for the same reduced scattering coefficient and refractive index of both layers. The derivation of this solution together with that for an infinite geometry is presented in the appendix.

In this section, we concentrate on the second method. First, solutions of the fluence rate in the frequency domain are derived numerically using the Fourier formalism that was explained in the accompanying paper.^{1} From these results the solutions in the time domain are derived using an FFT. For measurements in the frequency domain, the intensity of the incident source is modulated sinusoidally in time:

## Eq. 5

$$\Phi (\mathbf{r},t)={\int}_{-\infty}^{\infty}{\int}_{V}G({\mathbf{r}}^{\prime},{t}^{\prime})S(\mathbf{r}-{\mathbf{r}}^{\prime})\mathrm{exp}\left[i\omega (t-{t}^{\prime})\right]{\mathrm{d}}^{3}{r}^{\prime}\phantom{\rule{0.2em}{0ex}}\mathrm{d}{t}^{\prime}.$$Subsequently, the 3-D integration in the spatial domain is performed. By assuming a $\delta $ source in space, $S(\mathbf{r}-{\mathbf{r}}^{\prime})=\delta (\mathbf{r}-{\mathbf{r}}^{\prime})$ , we obtain

## Eq. 6

$$\Phi (\mathbf{r},t)={e}^{i\omega t}{\int}_{-\infty}^{\infty}G(\mathbf{r},{t}^{\prime})\mathrm{exp}(-i\omega {t}^{\prime})\mathrm{d}{t}^{\prime}=G(\mathbf{r},\omega ){e}^{i\omega t},$$## Eq. 7

$$(\Delta -{\alpha}^{2})G(\mathbf{r},\omega )=-\frac{1}{D}\delta \left(\mathbf{r}\right),\phantom{\rule{1em}{0ex}}{\alpha}^{2}=\frac{{\mu}_{a}}{D}+i\frac{\omega}{Dc},$$## Eq. 8

$$G(\mathbf{r},t)=\frac{1}{2\pi}{\int}_{-\infty}^{\infty}G(\mathbf{r},\omega ){e}^{i\omega t}\phantom{\rule{0.2em}{0ex}}\mathrm{d}\omega .$$For a vanishing intensity oscillation
$(\omega =0\phantom{\rule{0.3em}{0ex}}\mathrm{Hz})$
we get from Eq. 7 the steady-state diffusion equation, which was solved for
$N$
-layered turbid media in the accompanying paper.^{1} Thus, to obtain the solution of the
$N$
-layered diffusion equation in the frequency domain assuming a
$\delta $
-source at
$x=y=0$
and
$z={z}_{0}$
:

## Eq. 9

$$\Delta {\Phi}_{1}(x,y,z,\omega )-(\frac{{\mu}_{a1}}{{D}_{1}}+i\frac{\omega}{{c}_{1}{D}_{1}}){\Phi}_{1}(x,y,z,\omega )=-\frac{1}{{D}_{1}}\delta (x,y,z-{z}_{0}),\phantom{\rule{1em}{0ex}}0\u2a7dz<{l}_{1},$$## Eq. 10

$$\Delta {\Phi}_{k}(x,y,z,\omega )-(\frac{{\mu}_{ak}}{{D}_{k}}+i\frac{\omega}{{c}_{k}{D}_{k}}){\Phi}_{k}(x,y,z,\omega )=0,\phantom{\rule{1em}{0ex}}\sum _{j=1}^{k-1}{l}_{j}<z\u2a7d\sum _{j=1}^{k}{l}_{j},\phantom{\rule{1em}{0ex}}k=2,3,\dots ,N,$$^{1}can be used simply by replacing ${\mu}_{ak}\u2215{D}_{k}$ with ${\mu}_{ak}\u2215{D}_{k}+i\omega \u2215\left({c}_{k}{D}_{k}\right)$ , where ${c}_{k}$ is the velocity of light in layer $k$ . In addition, the four methods

^{1}presented to perform the numerical 2-D inverse Fourier transform as well as the approximate solutions in the steady-state domain can also be applied in the frequency domain.

The reflectance from the $N$ -layered turbid medium in the frequency domain is computed from the fluence rate by

## Eq. 11

$$R(\rho ,\omega )={\int}_{2\pi}\mathrm{d}\Omega [1-{R}_{f}\left(\theta \right)]\frac{1}{4\pi}[{\Phi}_{1}(\rho ,z=0,\omega )+3{D}_{1}{\phantom{|}\frac{\partial {\Phi}_{1}(\rho ,z,\omega )}{\partial z}|}_{z=0}\phantom{\rule{0.2em}{0ex}}\mathrm{cos}\left(\theta \right)]\mathrm{cos}\left(\theta \right),$$^{1}the Fresnel reflection coefficient and $\rho ={({x}^{2}+{y}^{2})}^{1\u22152}$ . For large distances from the source $(\rho \u2aa21\u2215{\mu}_{s}^{\prime})$ it is sufficient to use the Fick’s law for calculation of the reflectance

^{9}

## Eq. 12

$$R(\rho ,\omega )={D}_{1}{\phantom{|}\frac{\partial {\Phi}_{1}(\rho ,z,\omega )}{\partial z}|}_{z=0}.$$In the frequency domain, the intensity-modulated source causes that the reflectance measured at the detector position $\rho $ is also modulated by the same frequency $f$ but the oscillation is phase shifted and the modulation is reduced. The quantities of interest in the frequency domain are the phase angle $\theta $ between the source and the detected signal and the modulation $M$ :

## Eq. 13

$$\theta (\rho ,\omega )=\mathrm{arctan}\left\{\frac{\mathrm{Im}\left[R(\rho ,\omega )\right]}{\mathrm{Re}\left[R(\rho ,\omega )\right]}\right\},$$## Eq. 14

$$M(\rho ,\omega )={\left\{\frac{\mathrm{Re}{\left[R(\rho ,\omega )\right]}^{2}+\mathrm{Im}{\left[R(\rho ,\omega )\right]}^{2}}{\mathrm{Re}{\left[R(\rho ,\omega =0)\right]}^{2}}\right\}}^{1\u22152}.$$In summary, to obtain the solutions for the fluence rate in the frequency domain, the solutions in the steady-state domain presented in the accompanying paper^{1} are used by replacing
${\mu}_{ak}\u2215{D}_{k}$
with
${\mu}_{ak}+i\omega \u2215\left({c}_{k}{D}_{k}\right)$
. The solutions in the time domain are calculated by computing the solution in the frequency domain for many frequencies (e.g. 512) and by using^{10} an FFT. Note that because
$\Phi (\mathbf{r},t)$
is a causal function in time, it can be shown using the Hilbert transform that it is sufficient to calculate only the real part or the imaginary part of
$\Phi (\mathbf{r},\omega )$
. Then,
$\Phi (\mathbf{r},t)$
is obtained by an FFT of two times the real part or two times the imaginary part of
$\Phi (\mathbf{r},\omega )$
, accelerating the calculations by a factor of 2.

The solutions of the
$N$
-layered diffusion equation are compared to Monte Carlo simulations. The Monte Carlo method is within its stochastic nature an exact solution of radiative transfer theory, whereas diffusion theory is an approximation of radiative transfer theory. Thus, the validity of the diffusion equation for the considered geometrical and optical parameters can be tested by comparison with Monte Carlo simulations. In addition to the explanations given in the accompanying paper,^{1} for the calculation of the Monte Carlo simulations in the time domain we use a variance reduction scheme. The absorption coefficient is set to zero during the simulations for all layers and the length of the photon paths in each layer is scored. Thereby, the time-resolved reflectance can be calculated using the Lambert-Beer law.^{11} The results in the frequency domain are obtained by Fourier transforming (via an FFT) the Monte Carlo data calculated in the time domain.

## 3.

## Results

We, first, show a comparison of solutions of the diffusion equation obtained with different methods for computing the inverse Fourier transform and calculated with the formulas presented in the appendix. Second, the solutions of the diffusion equation are compared to Monte Carlo simulations.

## 3.1.

### Comparison of Different Methods for Solving the Diffusion Equation

In the accompanying paper,^{1} we presented four different methods for calculation of the inverse Fourier transform to obtain the steady-state fluence rate in real space. The accuracy of the computed solutions can be tested with these methods. All methods can be applied also in the frequency domain and, subsequently, for the calculation in the time domain. Figure 1
compares the results of two methods for calculating the 2-D inverse Fourier transform, the numerical integration of the Hankel transform (solid curve, method C in Ref. 1) and the use of the correlation for solving the Hankel transform (circles, method D in Ref. 1). Figure 1 shows the time-resolved reflectance from a two-layered turbid medium having an infinitely thick second layer at a distance
$\rho =22\phantom{\rule{0.3em}{0ex}}\mathrm{mm}$
.

Both curves are indistinguishable in the figure. Thus, Fig. 2 shows the relative difference between both curves. The differences are smaller than ${10}^{-6}$ for all time values over an intensity range of more than six orders of magnitude. A high accuracy of the obtained time-resolved fluence rate is important; e.g., for solving the inverse problem, the robust reconstruction of the optical properties of the different layers.

In the appendix, an analytical solution of the two-layered semi-infinite diffusion equation having the same reduced scattering coefficient and refractive index in both layers is derived, see Eq. 41. Figure 3 compares the reflectance calculated with this formula (circles) to that computed with the numerical solution presented in the theory section using method D in Ref. 1 for the inverse Fourier calculations (solid curve). The distance between the incident beam and detection is $\rho =30\phantom{\rule{0.3em}{0ex}}\mathrm{mm}$ .

The disagreement between the two curves are small, so that we also plot the relative differences (see Fig. 4 ). The relative differences are again smaller than ${10}^{-6}$ for the time values shown. Clearly, they can even be reduced by laying more stress on the numerical calculations, which on the other hand increases the computation time. The advantage of the analytical formula is its velocity. When the FFT is used for evaluation [see Eq. 43], the calculation time is in the range of a millisecond.

## 3.2.

### Comparison of the Solutions of the Diffusion Equation with Monte Carlo Simulations

First, we calculated the reflectance from a four-layered turbid medium having an infinitely thick fourth layer at a distance of $20.5\phantom{\rule{0.3em}{0ex}}\mathrm{mm}$ from the incident light source in the time domain. Figure 5 shows the comparison of the results obtained by the Monte Carlo method (circles) and by the diffusion solution (solid curve).

The two solutions agree well except for reflectance values at short times where it is known that diffusion theory is not valid. To show the small disagreement more distinctly the relative difference is depicted in Fig. 6 . For times longer than about $0.5\phantom{\rule{0.3em}{0ex}}\mathrm{ns}$ , the relative differences are mainly due to the statistics of the Monte Carlo simulations. The systematic relative differences between the solution of the diffusion equation and the Monte Carlo simulation are smaller than 1%.

In the frequency domain, we compare the results obtained with both models for the same four-layered turbid medium that was considered in the time domain. The geometrical and optical properties of the turbid medium are given in the caption of Fig. 5. Figures 7 and 8 show the phase and the modulation for the reflectance at a distance of $20.5\phantom{\rule{0.3em}{0ex}}\mathrm{mm}$ calculated by the Monte Carlo simulations (circles) and the solution of the diffusion equation (solid curve). The phase and modulation are drawn versus the modulation frequency $f$ . These data were directly calculated in the frequency domain in the case of the solution of the diffusion equation, whereas in the case of the Monte Carlo method simulations were performed in the time domain and the FFT was applied to obtain the phase and modulation for different modulation frequencies.

In general, good agreement between the two methods is obtained. The differences are larger for the phase, mainly caused by errors due to the diffusion approximations at small times (compare Fig. 5). The advantage in the time domain is that the photons that are remitted at early times and that do not fulfill the diffusion approximation can be separated from the photons at longer times that are well described by the diffusion approximation. Contrarily, in the frequency domain this separation is not possible.

## 4.

## Discussion

Solutions of the
$N$
-layered diffusion equation were derived for turbid media having a finite and infinite thick
$N$
’th layer and mismatched refractive indices. The final results in the frequency domain were obtained by a numerical Fourier transform of the analytical solutions of the fluence rate in the Fourier space. As explained in detail in our accompanying paper^{1} we employed four different methods to be able to validate the accuracy of the calculations. The results in the time domain were calculated by an FFT using the solutions in the frequency domain at many frequencies.

We also derived an analytical solution for the special case of a two-layered turbid medium having the same reduced scattering and refractive index in both layers, of which the bottom one is infinitely thick. By applying the FFT fast solutions in the range of $1\phantom{\rule{0.3em}{0ex}}\mathrm{ms}$ could be achieved with a state-of-the-art computer using one processor.

The accuracy of the numerical and analytical solutions was found to be smaller than ${10}^{-6}$ for the evaluated times and geometrical and optical parameters. This number is a compromise between the accuracy and velocity of the calculations.

The derived solutions were compared to Monte Carlo simulations in the time and frequency domains and similar agreement as for homogeneous geometries was found. Solutions of the diffusion equation in the frequency domain showed larger differences compared to those in the time domain for longer times due to the diffusion approximation as discussed in the results section.

In the future we plan to apply the derived solutions for investigating the inverse problem, the reconstruction of the optical properties for some or all involved layers. Finally, note that executable files of a computer code for the solutions of the
$N$
-layered diffusion equation in the time and steady-state domains are available on the Internet.^{12}

## Acknowledgment

We acknowledge the support by the European Union (nEUROPt, grant agreement no. 201076).

## References

**,” J. Biomed. Opt., 15 (2), 025003 (2010). 1083-3668 Google Scholar**

*Light diffusion in $N$-layered turbid media: steady state domain***,” Lasers Med. Sci., 6 379 –389 (1991). https://doi.org/10.1007/BF02042460 0268-8921 Google Scholar**

*The propagation of optical radiation in tissue. II: Optical properties of tissues and resulting fluence distributions***,” Inverse Probl., 15 R41 –R93 (1999). https://doi.org/10.1088/0266-5611/15/2/022 0266-5611 Google Scholar**

*Optical tomography in medical imaging inverse problems***,” Phys. Med. Biol., 50 2559 –2571 (2005). https://doi.org/10.1088/0031-9155/50/11/008 0031-9155 Google Scholar**

*Characterization of normal breast tissue heterogeneity using time-resolved near-infrared spectroscopy***,” Phys. Rev. Lett., 100 138101 (2008). https://doi.org/10.1103/PhysRevLett.100.138101 0031-9007 Google Scholar**

*Time-resolved diffuse reflectance using small source-detector separation and fast single-photon gating***,” J. Biomed. Opt., 11 044005 (2006). https://doi.org/10.1117/1.2337546 1083-3668 Google Scholar**

*In vivo absorption, scattering, and physiologic properties of 58 malignant breast tumors determined by broadband diffuse optical spectroscopy***,” J. Biomed. Opt., 13 041302 (2008). https://doi.org/10.1117/1.2967535 1083-3668 Google Scholar**

*Tutorial on diffuse light transport***,” J. Mod. Opt., 39 1567 –1582 (1992). https://doi.org/10.1080/09500349214551581 0950-0340 Google Scholar**

*Photon migration in a two-layer turbid medium. A diffusion analysis***,” J. Opt. Soc. Am. A, 14 246 –254 (1997). https://doi.org/10.1364/JOSAA.14.000246 0740-3232 Google Scholar**

*Improved solutions of the steady-state and the time-resolved diffusion equations for reflectance from a semi-infinite turbid medium***,” Phys. Med. Biol., 41 2221 –2227 (1996). https://doi.org/10.1088/0031-9155/41/10/026 0031-9155 Google Scholar**

*Determination of the optical properties of turbid media from a single Monte Carlo simulation***,” Appl. Opt., 28 2331 –2336 (1989). https://doi.org/10.1364/AO.28.002331 0003-6935 Google Scholar**

*Time-resolved reflectance and transmittance for the noninvasive measurement of tissue optical properties*## Appendices

### Appendix

In this section, we present the derivation of an analytical solution of the diffusion equation for the reflectance from a two-layered turbid medium having the same refractive index and the same reduced scattering coefficient in both layers. To this end, we apply the Laplace and Fourier formalism, which was explained in Sec. 2. As an introduction, we first present the same formalism for deriving the well-known solution in the infinite medium.

#### Appendix A: Infinite Turbid Medium

The starting point is the time-dependent diffusion equation for a delta source in space and time [see Eq. 2]. By applying the Laplace transform,

## Eq. 15

$$L\left\{G(x,y,z,t)\right\}=G(x,y,z,p)={\int}_{0}^{\infty}G(x,y,z,t){e}^{-pt}\phantom{\rule{0.2em}{0ex}}\mathrm{d}t,$$## Eq. 17

$$G({s}_{1},{s}_{2},z,p)={\int}_{-\infty}^{\infty}{\int}_{-\infty}^{\infty}G(x,y,z,p)\mathrm{exp}\left[i({s}_{1}x+{s}_{2}y)\right]\mathrm{d}x\phantom{\rule{0.2em}{0ex}}\mathrm{d}y,$$## Eq. 18

$$\frac{{\partial}^{2}G({s}_{1},{s}_{2},z,p)}{\partial {z}^{2}}-{\alpha}^{2}G({s}_{1},{s}_{2},z,p)=-\frac{1}{D}\delta \left(z\right),$$^{1}

## Eq. 21

$$\frac{\partial G({s}_{1},{s}_{2},z={0}^{-},p)}{\partial z}-\frac{\partial G({s}_{1},{s}_{2},z={0}^{+},p)}{\partial z}=\frac{1}{D}.$$The ordinary differential Eq. 18 can be solved using

## Eq. 22

$$G({s}_{1},{s}_{2},z,p)={c}_{1}{e}^{\alpha z}+{c}_{2}{e}^{-\alpha z}\phantom{\rule{1em}{0ex}}z<0,$$## Eq. 23

$$G({s}_{1},{s}_{2},z,p)={c}_{3}{e}^{\alpha z}+{c}_{4}{e}^{-\alpha z}\phantom{\rule{1em}{0ex}}z\u2a7e0.$$With the preceding four boundary conditions, the unknown constants are determined resulting in the following solutions:

## Eq. 25

$$G({s}_{1},{s}_{2},z,p)=\frac{1}{2\alpha D}{e}^{-\alpha z}\phantom{\rule{1em}{0ex}}z\u2a7e0.$$Inserting $\alpha $ in Eq. 24 yields

## Eq. 26

$$G({s}_{1},{s}_{2},z,p)={\left(\frac{c}{4D}\right)}^{1\u22152}\frac{{e}^{(-z\u2215\sqrt{Dc}){[p+c{\mu}_{a}+Dc({s}_{1}^{2}+{s}_{2}^{2})]}^{1\u22152}}}{{[p+c{\mu}_{a}+Dc({s}_{1}^{2}+{s}_{2}^{2})]}^{1\u22152}}\phantom{\rule{1em}{0ex}}z<0,$$## Eq. 27

$$G({s}_{1},{s}_{2},z,p)={\left(\frac{c}{4D}\right)}^{1\u22152}\frac{{e}^{(-z\u2215\sqrt{Dc}){[p+c{\mu}_{a}+Dc({s}_{1}^{2}+{s}_{2}^{2})]}^{1\u22152}}}{{[p+c{\mu}_{a}+Dc({s}_{1}^{2}+{s}_{2}^{2})]}^{1\u22152}}\phantom{\rule{1em}{0ex}}z\u2a7e0.$$Subsequently, these equations are two-dimensionally inverse Fourier and Laplace transformed to obtain $G(x,y,z,t)$ . The order of the three integrations can be chosen arbitrarily. We first apply the inverse Laplace transform and then two times the inverse Fourier transform. By using following formulas for the Laplace transform:

## Eq. 28

$${L}^{-1}\left[\frac{\mathrm{exp}(-k\sqrt{p})}{\sqrt{p}}\right]=\frac{1}{\sqrt{\pi t}}\phantom{\rule{0.2em}{0ex}}\mathrm{exp}(-{k}^{2}\u22154t),$$## Eq. 29

$$L\left[{e}^{-at}G\left(t\right)\right]=G(p+a),\phantom{\rule{1em}{0ex}}G\left(p\right)=L\left\{G\left(t\right)\right\},$$## Eq. 30

$$G({s}_{1},{s}_{2},z,t)={\left(\frac{c}{4D\pi t}\right)}^{1\u22152}\phantom{\rule{0.2em}{0ex}}\mathrm{exp}\{-[c{\mu}_{a}+Dc({s}_{1}^{2}+{s}_{2}^{2})]t\}\mathrm{exp}(-\frac{{z}^{2}}{4Dct}).$$Finally, the 2-D inverse Fourier transform is employed using the formula (and its corresponding equation with respect to $y$ and ${s}_{2}$ ):

## Eq. 31

$$F\left[\mathrm{exp}(-\frac{b{x}^{2}}{2})\right]=\frac{\sqrt{2\pi}}{\sqrt{b}}\phantom{\rule{0.2em}{0ex}}\mathrm{exp}(-\frac{{s}^{2}}{2b}),$$^{13}

#### Appendix B: Two-Layered Turbid Medium

Similar to the last subsection we started using the Laplace and Fourier transform to obtain an ordinary partial differential equation in the variable
$z$
. Then the boundary conditions for the two-layered turbid medium were applied. The resulting equation system was solved, as presented in the accompanying paper.^{1} There, we got for the fluence rate from a two-layered turbid medium having an infinitely thick second layer in the transformed space, assuming that
$D={D}_{1}={D}_{2}$
and
${n}_{1}={n}_{2}$
,

## Eq. 33

$$G(s,z,p)=\frac{\mathrm{sinh}\left[{\alpha}_{1}({z}_{0}+{z}_{b})\right]}{D{\alpha}_{1}}\frac{{\alpha}_{1}\phantom{\rule{0.2em}{0ex}}\mathrm{cosh}\left[{\alpha}_{1}({l}_{1}-z)\right]+{\alpha}_{2}\phantom{\rule{0.2em}{0ex}}\mathrm{sinh}\left[{\alpha}_{1}({l}_{1}-z)\right]}{{\alpha}_{1}\phantom{\rule{0.2em}{0ex}}\mathrm{cosh}\left[{\alpha}_{1}({l}_{1}+{z}_{b})\right]+{\alpha}_{2}\phantom{\rule{0.2em}{0ex}}\mathrm{sinh}\left[{\alpha}_{1}({l}_{1}+{z}_{b})\right]}-\frac{\mathrm{sinh}\left[{\alpha}_{1}({z}_{0}-z)\right]}{D{\alpha}_{1}},$$We make use of the following auxiliary equation

## Eq. 34

$$\frac{{\alpha}_{1}\phantom{\rule{0.2em}{0ex}}\mathrm{cosh}\left[{\alpha}_{1}({l}_{1}-z)\right]+{\alpha}_{2}\phantom{\rule{0.2em}{0ex}}\mathrm{sinh}\left[{\alpha}_{1}({l}_{1}-z)\right]}{{\alpha}_{1}\phantom{\rule{0.2em}{0ex}}\mathrm{cosh}\left[{\alpha}_{1}({l}_{1}+{z}_{b})\right]+{\alpha}_{2}\phantom{\rule{0.2em}{0ex}}\mathrm{sinh}\left[{\alpha}_{1}({l}_{1}+{z}_{b})\right]}=\frac{({\alpha}_{1}+{\alpha}_{2})\mathrm{exp}\left[{\alpha}_{1}({l}_{1}-z)\right]+({\alpha}_{1}-{\alpha}_{2})\mathrm{exp}-{\alpha}_{1}({l}_{1}-z)}{({\alpha}_{1}+{\alpha}_{2})\mathrm{exp}\left[{\alpha}_{1}({l}_{1}+{z}_{b})\right]+({\alpha}_{1}-{\alpha}_{2})\mathrm{exp}[-{\alpha}_{1}({l}_{1}+{z}_{b})]}=\frac{\mathrm{exp}\left[{\alpha}_{1}({l}_{1}-z)\right]+({\alpha}_{1}-{\alpha}_{2}\u2215{\alpha}_{1}+{\alpha}_{2})\mathrm{exp}[-{\alpha}_{1}({l}_{1}-z)]}{\mathrm{exp}\left[{\alpha}_{1}({l}_{1}+{z}_{b})\right]+({\alpha}_{1}-{\alpha}_{2}\u2215{\alpha}_{1}+{\alpha}_{2})\mathrm{exp}[-{\alpha}_{1}({l}_{1}+{z}_{b})]}=\mathrm{exp}[-{\alpha}_{1}(z+{z}_{b})]\frac{1+{D}_{a}\phantom{\rule{0.2em}{0ex}}\mathrm{exp}[-2{\alpha}_{1}({l}_{1}-z)]}{1+{D}_{a}\phantom{\rule{0.2em}{0ex}}\mathrm{exp}[-2{\alpha}_{1}({l}_{1}+{z}_{b})]}=\mathrm{exp}[-{\alpha}_{1}(z+{z}_{b})]\{1+{D}_{a}\phantom{\rule{0.2em}{0ex}}\mathrm{exp}[-2{\alpha}_{1}({l}_{1}-z)]\}\sum _{m=0}^{\infty}{(-{D}_{a})}^{m}\phantom{\rule{0.2em}{0ex}}\mathrm{exp}[-2m{\alpha}_{1}({l}_{1}+{z}_{b})],$$## Eq. 35

$$G(s,z,p)=\frac{\mathrm{exp}\left[{\alpha}_{1}({z}_{0}-z)\right]-\mathrm{exp}[-{\alpha}_{1}({z}_{0}+z+2{z}_{b})]}{2D{\alpha}_{1}}(1+\{1-\mathrm{exp}\left[2{\alpha}_{1}(z+{z}_{b})\right]\}\sum _{m=1}^{\infty}{(-{D}_{a})}^{m}\phantom{\rule{0.2em}{0ex}}\mathrm{exp}[-2m{\alpha}_{1}({l}_{1}+{z}_{b})])-\frac{\mathrm{exp}\left[{\alpha}_{1}({z}_{0}-z)\right]-\mathrm{exp}[-{\alpha}_{1}({z}_{0}-z)]}{2D{\alpha}_{1}}.$$## Eq. 36

$$G(s,z,p)=\frac{\mathrm{exp}\left[{\alpha}_{1}({z}_{0}-z)\right]-\mathrm{exp}[-{\alpha}_{1}({z}_{0}+z+2{z}_{b})]}{2D{\alpha}_{1}}+\frac{1}{2D{\alpha}_{1}}\sum _{m=1}^{\infty}\left[(\mathrm{exp}\{-{\alpha}_{1}[{z}_{0}-z+2m(l+{z}_{b})]\}-\mathrm{exp}\{-{\alpha}_{1}[{z}_{0}+z+2{z}_{b}+2m(l+{z}_{b})]\}+\mathrm{exp}\{-{\alpha}_{1}[-{z}_{0}+z+2m(l+{z}_{b})]\}-\mathrm{exp}\{-{\alpha}_{1}[-{z}_{0}-z-2{z}_{b}+2m(l+{z}_{b})]\}){(-{D}_{a})}^{m}\right].$$Subsequently, we perform the inverse transforms. Employing the Laplace transform formulas^{14}

## Eq. 37

$$\frac{\mathrm{exp}[-a{(p+b)}^{1\u22152}]}{{(p+b)}^{1\u22152}},\phantom{\rule{1em}{0ex}}a>0\iff \frac{1}{\sqrt{\pi t}}{e}^{-bt}\phantom{\rule{0.2em}{0ex}}\mathrm{exp}(-\frac{{a}^{2}}{4t})\epsilon \left(t\right),$$## Eq. 38

$${\left[\frac{{(p+a)}^{1\u22152}-{(p+b)}^{1\u22152}}{{(p+a)}^{1\u22152}+{(p+b)}^{1\u22152}}\right]}^{m},\phantom{\rule{1em}{0ex}}m>0\iff \frac{m}{t}\phantom{\rule{0.2em}{0ex}}\mathrm{exp}[-\frac{1}{2}(a+b)t]{I}_{m}\left(\frac{a-b}{2}t\right)\epsilon \left(t\right),$$## Eq. 39

$$X\left(p\right)\cdot Y\left(p\right)\iff X\left(t\right)\ast Y\left(t\right)={\int}_{0}^{t}X\left(\tau \right)Y(t-\tau )\mathrm{d}\tau ,$$## Eq. 40

$$\mathrm{exp}(-Dct\cdot {s}_{1}^{2})\iff \frac{1}{{\left(4\pi Dct\right)}^{1\u22152}}\phantom{\rule{0.2em}{0ex}}\mathrm{exp}(-\frac{{x}^{2}}{4Dct}),$$## Eq. 41

$$G(\mathbf{r},t)=\frac{c}{{\left(4\pi Dct\right)}^{3\u22152}}\phantom{\rule{0.2em}{0ex}}\mathrm{exp}(-{\mu}_{a1}ct)\mathrm{exp}(-\frac{{\rho}^{2}}{4Dct})\{\mathrm{exp}[-\frac{{(z-{z}_{0})}^{2}}{4Dct}]-\mathrm{exp}[-\frac{{(z+{z}_{0}+2{z}_{b})}^{2}}{4Dct}]\}+\frac{1}{4\pi Dct}\phantom{\rule{0.2em}{0ex}}\mathrm{exp}(-\frac{{\rho}^{2}}{4Dct}){\int}_{0}^{t}\frac{1}{2}{\left[\frac{c}{D\pi (t-\tau )}\right]}^{1\u22152}\phantom{\rule{0.2em}{0ex}}\mathrm{exp}[-{\mu}_{a1}c(t-\tau )]\mathrm{exp}[-\frac{1}{2}({\mu}_{a1}+{\mu}_{a2})c\tau ]\times \sum _{m=-\infty}^{\infty}(\mathrm{exp}\{-\frac{{[z+{z}_{0}+2{z}_{b}+2m({l}_{1}+zb)]}^{2}}{4Dc(t-\tau )}\}-\mathrm{exp}\{-\frac{{[{z}_{0}-z+2m({l}_{1}+zb)]}^{2}}{4Dc(t-\tau )}\})\frac{\left|m\right|}{\tau}{I}_{\left|m\right|}\left(\frac{{\mu}_{a2}-{\mu}_{a1}}{2}c\tau \right)\mathrm{d}\tau .$$Equation 41 was evaluated by Gauss integration (see Fig. 3). The first term of Eq. 41 is identical to the solution of the diffusion equation for a semi-infinite homogeneous geometry.

To accelerate the calculations we derive an alternative formula for $G(\mathbf{r},t)$ . We start with Eq. 36, the first part of which is inverse Laplace and Fourier transformed as in the preceding, resulting in the first term of $G(\mathbf{r},t)$ given in Eq. 41. The second part of Eq. 36 is transformed as follows:

## Eq. 42

$$\frac{1}{2D{\alpha}_{1}}\{\mathrm{exp}[-{\alpha}_{1}({z}_{0}-z)]-\mathrm{exp}[-{\alpha}_{1}({z}_{0}+z+2{z}_{b})]+\mathrm{exp}\left[{\alpha}_{1}({z}_{0}-z)\right]-\mathrm{exp}\left[{\alpha}_{1}({z}_{0}+z+2{z}_{b})\right]\}\sum _{m=1}^{\infty}{\{-{D}_{a}\phantom{\rule{0.2em}{0ex}}\mathrm{exp}[-2{\alpha}_{1}(l+{z}_{b})]\}}^{m}=\frac{1}{2D{\alpha}_{1}}\{\mathrm{cosh}\left[{\alpha}_{1}({z}_{0}-z)\right]-\mathrm{cosh}\left[{\alpha}_{1}({z}_{0}+z+2{z}_{b})\right]\}\{\frac{1}{1+{D}_{a}\phantom{\rule{0.2em}{0ex}}\mathrm{exp}[-2{\alpha}_{1}(l+{z}_{b})]}-1\}=\frac{1}{2D{\alpha}_{1}}\{\mathrm{cosh}\left[{\alpha}_{1}({z}_{0}+z+2{z}_{b})\right]-\mathrm{cosh}\left[{\alpha}_{1}({z}_{0}-z)\right]\}\frac{{D}_{a}}{{D}_{a}+\mathrm{exp}\left[2{\alpha}_{1}(l+{z}_{b})\right]}.$$Then, Eq. 42 is inverse Fourier transformed with respect to the $x$ and $y$ coordinates, whereas the inverse transform with respect to the time variable is not analytically evaluated. For the calculation of the inverse Fourier transform, we apply the corresponding shifting rule [see Eq. 29] with $a=Dc{s}^{2}$ . In the inverse Laplace transform with respect to the time variable we used $p=i\omega $ , because the imaginary axis is within the region of convergence. Thus, we obtain