Monte Carlo simulations for solving the radiative transfer equation

Light propagation is governed by Maxwell’s equations in terms of electric and magnetic fields. In media with a complex spatial distribution of the refractive index, often classified as random media, Maxwell’s equations can be reduced to a radiative transport equation. In random media, interference effects are neglected and the superposition of intensities rather than fields is considered. Most biological tissues act as random media, and the radiative transfer equation (RTE) gives the most complete description of light propagation in tissue:

\(\hat{s}\bullet\nabla L\left(\mathbf{r},\hat{s}\right)+\mu_t\left(\mathbf{r}\right)L\left(\mathbf{r},\hat{s}\right)=\mu_s\left(\mathbf{r}\right)\int_{4\pi} p\left({\hat{s}}^\prime\rightarrow\hat{s}\right)L\left(\mathbf{r},\hat{s}\prime\right)d{\hat{s}}^\prime+q\left(\mathbf{r},\hat{s}\right)\) 

The radiance \(L\left(\mathbf{r},\hat{s}\right)\) represents the power per unit solid angle and per unit area propagating along direction \(\hat{s}\) at a point \(\mathbf{r}\) inside the medium. The spatial distribution of the radiance is governed by: 1) collision events, described by the term \(\mu_t\left(\mathbf{r}\right)I\left(\mathbf{r},\hat{s}\right)\) (loss term), which accounts for absorption events and scattering events that deflect the light away from direction \(\hat{s}\) (the extinction coefficient, \(\mu_t\left(\mathbf{r}\right)=\mu_a\left(\mathbf{r}\right)+\mu_s\left(\mathbf{r}\right)\), is the sum of the absorption and scattering coefficients); 2) scattering events that redirect light from an arbitrary direction \(\hat{s}\prime\) into the direction \(\hat{s}\),as described by the term \(\mu_s\left(\mathbf{r}\right)\int_{4\pi} p\left({\hat{s}}^\prime\rightarrow\hat{s}\right)L\left(\mathbf{r},\hat{s}\prime\right)d{\hat{s}}^\prime\) (gain term); 3) light sources, described by \(q\left(\mathbf{r},\hat{s}\right)\) (gain term). How can one solve such a complex integro-differential equation? Monte Carlo and random walk methods are stochastic approaches that have been used since the early 1950s for solving partial differential equations (PDE) by using probability models of the governing equations. For the case of the RTE, the probability model uses random numbers to select the pathlength between collision events and the direction after each scattering event. By running many case histories (typical numbers are in the order of 10⁸-10⁹) one can calculate physical parameters of interests such as the fluence rate (obtained by integrating the radiance in all directions) at all points in the medium, or the light collected by a detector at the boundary of the medium. 

Figure 1 shows a Monte Carlo simulation of the spatial distribution of the fluence rate (represented in false colors) generated in a scattering medium by a unit-power laser beam that is incident on its surface. A representative photon path is also illustrated by the white line, showing the broken path resulting from scattering, and the end of the path resulting from absorption.
 

Fig. 1. Fluence rate per unit source power (color map) generated by a laser beam in a scattering medium, and representative photon path (white line) obtained with Monte Carlo simulations.