## 2 Materials and Methods  
### Diffusion Theory  
#### 2.1.1 Diffuse Light in Three Dimensions  
Let us suppose that a random medium occupies the three-dimensional half space. We assume that the random medium is characterized by the diffusion coefficient \(D\) and absorption coefficient \(\mu_{a}\). We have \(D=1/(3\mu_{s}^{\prime})\), where \(\mu_{s}^{\prime}\) is the reduced scattering coefficient. Position in the half space (\(-\infty<x<\infty\), \(-\infty<y<\infty\), \(0<z<\infty\)) is denoted by \(\mathbf{r}=(\boldsymbol{\rho},z)\), \(\boldsymbol{\rho}=(x,y)\). Let \(t\) denote time. Let \(c\) be the speed of light in the medium. We assume an incident beam on the boundary at the origin in the \(x\)-\(y\) plane. The energy density \(u\) of light in the medium obeys the following diffusion equation.  
\[\left(\frac{1}{c}\frac{\partial}{\partial t}-D\Delta+\mu_{a}\right)u(\mathbf{ r},t)=0, \tag{1}\]with the boundary condition  
\[-\ell\frac{\partial}{\partial z}u(\mathbf{r},t)+u(\mathbf{r},t)=\delta(x)\delta(y )q(t), \tag{2}\]  
and the initial condition \(u(\mathbf{r},0)=0\). The right-hand side of the boundary condition is the incident laser beam which illuminates the medium at the origin \((0,0,0)\) with the temporal profile \(q(t)\). In the phantom experiment described below, the incident light enters the phantom in the positive \(z\) direction. The extrapolation distance \(\ell\) is a nonnegative constant. If we consider the diffuse surface reflection, we have [25]  
\[\ell=2D\frac{1+r_{d}}{1-r_{d}}, \tag{3}\]  
where the internal reflection \(r_{d}\) is given by [26]\(r_{d}=-1.4399n^{-2}+0.7099n^{-1}+0.6681+0.0636n\). Let us consider the corresponding surface Green's function \(G_{\mathrm{s}}(\mathbf{r},t;\boldsymbol{\rho}^{\prime},s)\), which satisfies Equation (1) and the boundary condition  
\[\left(-\ell\frac{\partial}{\partial z}+1\right)G_{\mathrm{s}}(\mathbf{r},t; \boldsymbol{\rho}^{\prime},s)=\delta(x-x^{\prime})\delta(y-y^{\prime})\delta (t-s), \tag{4}\]  
together with the initial condition \(G_{\mathrm{s}}(\mathbf{r},0;\boldsymbol{\rho}^{\prime},s)=0\). We obtain [27; 28; 29; 30]  
\[\begin{split} G_{\mathrm{s}}(\mathbf{r},t;\boldsymbol{\rho}^{ \prime},s)&=\frac{cD}{\ell}\Bigg{[}\frac{2e^{-\mu_{\mathrm{d}}c(t -s)}}{\left(4\pi Dc(t-s)\right)^{3/2}}e^{-\frac{(x-y^{\prime})^{2}+(y-y^{ \prime})^{2}}{4 Dc(t-s)}}e^{-\frac{2}{4 Dc(t-s)}}\\ &-\frac{e^{-\mu_{\mathrm{d}}c(t-s)}}{4\pi\ell Dc(t-s)}e^{-\frac{ (x-y^{\prime})^{2}+(y-y^{\prime})^{2}}{4 Dc(t-s)}}e^{\frac{2}{4}}e^{\frac{Dc(t -s)}{\ell^{2}}}\operatorname{erfc}\left(\frac{z+2 Dc(t-s)/\ell}{\sqrt{4 Dc(t-s) }}\right)\Bigg{]},\quad t>s,\end{split} \tag{5}\]  
with \(G_{\mathrm{s}}(\mathbf{r},t;\boldsymbol{\rho}^{\prime},s)=0\) for \(t\leq s\). The complementary error function is defined by \(\operatorname{erfc}(x)=(2/\sqrt{\pi})\int_{x}^{\infty}\exp(-t^{2})\,dt\). Let \(\mathbf{r}_{d}\) be a point in the \(x\)-\(y\) plane with \(|\mathbf{r}_{d}|>0\). We obtain  
\[\begin{split}& u(\mathbf{r}_{d},t)=\int_{0}^{t}G_{\mathrm{s}}( \mathbf{r}_{d},t;\mathbf{0},s)q(s)\,ds\\ &=\int_{0}^{t}e^{-\mu_{\mathrm{d}}c(t-s)}e^{-\frac{|\mathbf{r}_{ d}|^{2}}{4 Dc(t-s)}}\frac{2q(s)}{\left(4\pi Dc(t-s)\right)^{3/2}}\left[1-\frac{ \sqrt{4\pi Dc(t-s)}}{2\ell}e^{\frac{Dc(t-s)}{\ell^{2}}}\operatorname{erfc} \left(\frac{\sqrt{Dc(t-s)}}{\ell}\right)\right]\,ds.\end{split} \tag{6}\]  
This \(u(\mathbf{r}_{d},t)\) in Equation (6) is used for the calculation in Section 3.1.