## 2 Materials and Methods
### Diffusion Theory
#### 2.1.2 Diffuse Light in Two Dimensions  
For the later purpose of a numerical experiment, we here consider light propagation in two dimensions. Let us suppose that the two-dimensional half space of positive \(y\) is occupied with a random medium in the \(x\)-\(y\) plane. The energy density \(u\) of light at position \(\boldsymbol{\rho}=(x,y)\) in the medium due to the incident beam on the boundary (the \(x\)-axis) at \(\boldsymbol{\rho}_{\mathrm{s}}^{i}\) (\(i=1,\ldots,M_{\mathrm{s}}\)) obeys the following diffusion equation. \[\left(\frac{1}{c}\frac{\partial}{\partial t}-D\Delta+\mu_{a}(\boldsymbol{\rho} )\right)u(\boldsymbol{\rho},t;\boldsymbol{\rho}_{\mathrm{s}}^{i})=0, \tag{7}\]  
with the boundary condition  
\[-\ell\frac{\partial}{\partial y}u(\boldsymbol{\rho},t;\boldsymbol{\rho}_{ \mathrm{s}}^{i})+u(\boldsymbol{\rho},t;\boldsymbol{\rho}_{\mathrm{s}}^{i})= \delta(x-x_{\mathrm{s}}^{i})\delta(t), \tag{8}\]and the initial condition \(u(\mathbf{\rho},0;\mathbf{\rho}_{s}^{i})=0\). The incident laser beam on the right-hand side of the boundary condition was assumed to be a pulse at the position \(\mathbf{\rho}_{s}^{i}=(x_{s}^{i},0)\). We suppose that the diffusion coefficient \(D\) is a positive constant but \(\mu_{a}\) varies in space. Below we develop the Rytov approximation. Let us write \(\mu_{a}(\mathbf{\rho})\geq 0\) as  
\[\mu_{a}(\mathbf{\rho})=\mu_{a0}+\delta\mu_{a}(\mathbf{\rho}), \tag{9}\]  
where \(\mu_{a0}\) is a constant and the perturbation \(\delta\mu_{a}(\mathbf{\rho})\) spatially varies. We note the relation,  
\[u(\mathbf{\rho},t;\mathbf{\rho}_{s}^{i})=u_{0}(\mathbf{\rho},t;\mathbf{\rho}_{s}^{i})-\int_{0} ^{t}\int_{0}^{\infty}\int_{-\infty}^{\infty}G(\mathbf{\rho},t;\mathbf{\rho}^{\prime},s )\delta\mu_{a}(\mathbf{\rho}^{\prime})u(\mathbf{\rho}^{\prime},s;\mathbf{\rho}_{s}^{i})\, dx^{\prime}dy^{\prime}ds, \tag{10}\]  
where \(u_{0}(\mathbf{\rho},t;\mathbf{\rho}_{s}^{i})\) is the solution to the diffusion Equation (7) with the zeroth-order coefficient replaced by \(b_{0}\). We introduce the Green's function \(G(\mathbf{\rho},t;\mathbf{\rho}^{\prime},s)\), which satisfies  
\[\left(\frac{1}{c}\frac{\partial}{\partial t}-D\Delta+\mu_{a0}\right)G=\delta( \mathbf{\rho}-\mathbf{\rho}^{\prime})\delta(t-s), \tag{11}\]  
with the boundary condition \(-\ell\frac{\partial G}{\partial y}+G=0\) at \(y=0\), and the initial condition \(G=0\) at \(t=0\). The above relation (10) can be directly verified. By recursive substitution, we obtain \(u\) as  
\[u(\mathbf{\rho},t;\mathbf{\rho}_{s}^{i})=u_{0}(\mathbf{\rho},t;\mathbf{\rho}_{s}^{i})-\int_{0} ^{t}\int_{0}^{\infty}\int_{-\infty}^{\infty}G(\mathbf{\rho},t;\mathbf{\rho}^{\prime},s )\delta\mu_{a}(\mathbf{\rho}^{\prime})u_{0}(\mathbf{\rho}^{\prime},s;\mathbf{\rho}_{s}^{i })\,dx^{\prime}dy^{\prime}ds+O((\delta\mu_{a})^{2}). \tag{12}\]  
By neglecting higher-order terms assuming that \(\delta\mu_{a}\) is small, we arrive at the (first) Born approximation [31], in which \(u\) is given by  
\[u(\mathbf{\rho},t;\mathbf{\rho}_{s}^{i})=u_{0}(\mathbf{\rho},t;\mathbf{\rho}_{s}^{i})-c\int_{ 0}^{t}\int_{0}^{\infty}\int_{-\infty}^{\infty}G(\mathbf{\rho},t;\mathbf{\rho}^{\prime },s)\delta\mu_{a}(\mathbf{\rho}^{\prime})u_{0}(\mathbf{\rho}^{\prime},s;\mathbf{\rho}_{s}^ {i})\,dx^{\prime}dy^{\prime}ds. \tag{13}\]  
We note that the Green's function is obtained as [27; 28; 29; 30]  
\[ \begin{split} G(\mathbf{\rho},t;\mathbf{\rho}^{\prime},s)&= \frac{e^{-\mu_{a0}c(t-s)}}{4\pi D(t-s)}e^{-\frac{(x-y^{\prime})^{2}}{4 Dc(t-s)} }\Bigg{[}e^{-\frac{(y-y^{\prime})^{2}}{4 Dc(t-s)}}+e^{-\frac{(y+y^{\prime})^{2 }}{4 Dc(t-s)}}\\ &-\frac{\sqrt{4\pi Dc(t-s)}}{\ell}e^{-\frac{(y+y^{\prime})^{2}}{ 4 Dc(t-s)}}e^{\frac{(y+y^{\prime}+2 Dc(t-s)/\ell}{\sqrt{4 Dc(t-s)}}}\Bigg{)}^{2} \text{erfc}\left(\frac{y+y^{\prime}+2 Dc(t-s)/\ell}{\sqrt{4 Dc(t-s)}}\right) \Bigg{]},\end{split}  \tag{14}\]  
for \(t>s\), and otherwise \(G(\mathbf{\rho},t;\mathbf{\rho}^{\prime},s)=0\).