# Time-Domain Near-Infrared Spectroscopy and Imaging: A Review
## 2 Theoretical Background of TD-NIRS
### Monte Carlo Simulations  
Instead of directly solving the RTE or the DE with analytical and numerical methods, stochastic methods have been used to obtain statistically equivalent solutions of the RTE or the DE. One of them is random walk theory, which describes light propagation as a sequence of steps on a regular cubic lattice [21]. But this method oversimplifies the phenomenon, and the most common stochastic methods are Monte Carlo (MC) simulations. [PERSON] and [PERSON] [22] first introduced an MC simulation into this field, and [PERSON] et al. [23] provided free software for an MC simulation to calculate light propagation in multi-layered tissue with their published paper explaining the details. Now free software codes for MC simulations are available from the website of the Oregon Laser Center in Portland [24], and these codes are used by many researchers in this field.  
The concept of MC simulations is to track the trajectories of energy packets, which are often called simply photons for brevity, while the trajectories of the photons are determined to statistically satisfy the optical properties of the medium, as shown in the following equations for Figure 4. The injected photon with an energy of \(E_{i}\) is scattered with a scattering path length of \(L\) to the direction of the polar and azimuthal angles of \(\theta\) and \(\varphi\), respectively, and the energy is attenuated to \(E_{i+1}=E_{i}[\mu_{s}/(\mu_{s}+\mu_{a})]\). \(L\), \(\theta\), and \(\varphi\) are determined using uniformly distributed random numbers in the range of \((0,1)\), \(R_{1}\), \(R_{2}\), and \(R_{3}\) by the following equations,  
\[L=-\frac{\ln\left(R_{1}\right)}{\mu_{t}}=-\frac{\ln\left(R_{1}\right)}{\mu_{s }+\mu_{a}},\ \theta=f_{2}^{-1}(R_{2}),\ f_{2}(\theta)=\frac{1}{2}\int_{0}^{\theta}p( \theta t)\sin\theta td\theta,\ \varphi=2\pi R_{3}, \tag{16}\]  
where the phase function \(p(\theta)\) is assumed to be independent of \(\varphi,f_{2}(\theta)\) is the accumulated phase function of \(p(\theta),f_{2}{}^{-1}\) means the inverse of \(f_{2}(\theta)\). When \(p(\theta)\) is constant, meaning that \(g=0\) and the diffusion approximation is applicable, \(f_{2}(\theta)\) becomes as simple as \(f_{2}(\theta)=(1-\cos\theta)/2\), and resultantly \(\theta\) is simply determined by \(\theta=\cos^{-1}(1-R_{2})\).  
Monte Carlo simulations can equivalently provide solutions to the RTE even under the conditions for which the DE breaks down, and due to its flexibility to geometry, [PERSON] et al. [25] succeeded in developing an MC code for predicting 3D and TD light propagation in the adult human head, which is a complex heterogeneous medium including a non-scattering and non-absorbing cerebro-spinal fluid (CSF) layer. However, for obtaining accurate solutions, long computation times are needed because the statistical noises are inversely proportional to the square root of the number of injected photons. Some  
Figure 4: Schematic of a Monte Carlo simulation.  
techniques have been studied to overcome this drawback of MC simulations and will be reviewed in a later section.