1.

Introduction

In the last decade, considerable attention has been given to the investigation of biological tissue by using novel photonic-based imaging and diagnostic modalities.¹ Optical techniques offer a strong advantage given that diagnostic procedures can be performed noninvasively on biotissues in situ. Optical diagnostic systems illuminate biological tissues with UV, visible, or infrared light and attempt to analyze the resulting reflectance, fluorescence emission, or Raman spectra.² Recently, for the same purpose, the polarization of light scattered in biological tissues has gained much attention.³ Various aspects of polarized light scattered in tissues as well as in other tissue-like scattering media have been extensively studied in the past.^3–6 Nevertheless, the development and practical implementation of polarization-based optical systems for day-to-day diagnostic routines, as well as interpretation and quantitative analysis of the experimental results, require an accurate description of polarized light propagation in biological tissues. Due to the complex structure of biological tissues and high scattering of light therein, the analytical approaches describing radiative transfer are impractical and cannot be applied. Stochastic alternatives, such as Monte Carlo (MC) methods,⁷ are ideally suited for modeling of optical radiation propagation in complex scattering media.

The conventional MC approach is comprehensively described elsewhere (see, e.g., Refs. 8 and 9), and is based on the modeling of energy or intensity transfer through the medium. It does not take into account the wave nature of light, i.e., polarization, coherence, phase retardation, interference, and associated wave phenomena. Recently, a new generation of stochastic MC techniques has been developed and widely employed in the studies of coherent effects of multiple scattering, such as enhancement of coherent backscattering (CBS) and decay of the temporal correlation function of scattered light.^10,11 MC has been used to study the depolarization of the backscattered electromagnetic field¹² and for calculation of low-coherence enhanced backscattering.¹³ The amplitude scattering matrix for CBS has been computed with the electric field MC,¹⁴ and by using the Jones $N$-matrix formalism, the implementation of birefringent properties of the medium has been performed.¹⁵ MC approach has been extended and extensively used for imitation of two-dimensional images of human skin obtained by optical coherence tomography.^16–19

In the current report, we present further development of the unified MC approach²⁰ for modeling of coherent polarized light propagation in a complex multiple scattering medium. The developed model utilizes Jones formalism specially adapted for the online object oriented MC (O3MC) approach developed earlier.^20,21 The model takes into account the wave properties of light, including temporal coherence, polarization, interference, and reflection/refraction at the medium boundary, as well as the helicity flip of the backscattered circularly and/or elliptically polarized light.

2.

GPU Accelerated MC Modeling of Photon Migration in Scattering Media

The conventional MC model of photon migration in a random scattering medium^8,9 is based on simulation of a large number of photon packet trajectories ($\sim {10}^6$ to ${10}^9$) consisting of $N_{\mathrm{s}}$ scattering events $(N_{\mathrm{s}}\sim {10}^2–{10}^3)$. Typically, MC is used to simulate diffuse reflectance, fluence rate distribution, sampling volume, reflectance spectra, etc.⁷ Depending on the computational performance of available facilities, this procedure can be extremely time-consuming (in a range from hours to several days).

Recently, it has been demonstrated that MC approach can be effectively sped up, allowing the results of simulations to be achieved in near real-time.^20,21 The developed O3MC code utilizes the compute unified device architecture (CUDA) parallel computing platform (recently introduced by NVIDIA Corporation, Santa Clara, California) that allows simultaneous execution of thousands of lightweight parallel threads.²² Each thread is responsible for simulating the propagation of one photon packet in the scattering medium. The block diagram showing the memory allocation scheme for storing photon packet trajectories is presented in Fig. 1.

Fig. 1

Schematic presentation of the online object oriented Monte Carlo (O3MC) graphics processing unit (GPU) memory allocation for storing the photon packet trajectories.

Each photon trajectory (Fig. 2) is represented by a one-dimensional (1-D) array of structures and is stored in the GPU’s global memory (see Fig. 1). Each structure consists of three floating point numbers $\mathbf{r}(i)=X_i,Y_i,Z_i$ corresponding to the $i$’th scattering event. The length of the trajectory is limited by the maximum number of scattering events $N_{\mathrm{s}}$, which is in our case ${10}^3$, and depending on the particular task can be extended to ${10}^4$ to ${10}^5$. Tesla M2090 graphics processing unit (GPU) supports simultaneous execution of a maximum 24,576 resident threads per GPU cycle. The memory required for storing these generated trajectories is $\sim 2.74\text{GB}$ in the case of single floating point performance (FPP), and $\sim 5.49\text{GB}$ for double FPP. Thus, the storage of ${10}^9$ photon trajectories to obtain good statistics requires an allocation of $\sim 109\text{PB}$ of global memory, that is, obviously beyond current hardware capabilities.

Fig. 2

Stochastic trajectory of a photon packet traveling from the point of incidence (source) to the point of detection (detector); ${\mathbf{r}}_{i-1}$ and ${\mathbf{r}}_i$ are the $i-1$ and the $i$ scattering events, and $l$ is the photon packet’s pathlength at the $i-1$ scattering event.

To achieve optimal memory allocation, we use a data structure that contains only the parameters of each photon packet that are absolutely necessary for computation of particular quantities of the detected light intensity, e.g., the number of scattering events $N_{\mathrm{s}}$, total path length $L$, state of polarization, number of helicity flips, etc. Note the state of polarization and total pathlength are floating point numbers, whereas the scattering orders and the number of helicity flips are integers. The minimal amount of memory required to store the state of polarization, total pathlength, and number of helicity flips requires 24 bytes for single or 44 bytes for double floating point precision. In order to calculate interference of a large number of photon packets (in our case ${10}^9$ and larger) the allocated memory is approximately 20 GB (${10}^9\times 24$) and 40 GB (${10}^9\times 44$) for single and double FPP, respectively. It is obvious that with the current developments of CUDA technology, it is impossible to allocate such an amount of GPU memory.

The most recent generation of professional GPUs supporting CUDA, which are based on the NVIDIA Kepler architecture (e.g., Tesla $K20X$, $K20$, $K10$), provides a maximum of 8 GB of fast onboard GDDR5 memory. This is still less than the above requirements. A number of approaches, such as tiling and cutoff,²² have been employed to effectively deal with the large data sets and avoid memory bandwidth limitations. For the same reason, in the framework of the current model, we apply another type of memory called page-locked or pinned memory offered by CUDA.²² THE page-locked memory is a specifically allocated RAM that is mapped to the GPU. The main advantage of this special type of memory is the concurrent access between RAM and the GPU’s global memory during a kernel execution. This allows the organization of a so-called memory conveyor.²³ Thus, the GPU thread performs computations over a chunk of data, while the rest is synchronized.²² In addition, the data structure alignment strategy has been utilized to eliminate the negative impact of random memory access across multiple threads, and to ensure that the GPU coalescence of multiple memory loads and stores is effectively working.²²

In older GPUs (i.e., $<2.0$ compute capability), the memory transactions are coalesced within half warps of 16 threads. Since then GPUs have a large $L1$ cache in each multiprocessor with a 128-byte line size, and memory transactions for whole warps of 32 threads are coalesced (e.g., Tesla M2090).²⁴ Thus, access to the global memory has been significantly improved, as well as the memory bandwidth. For instance, Tesla M2090s global memory is GDDR5 and its bandwidth is up to $177\mathrm{GB}/\mathrm{s}$ compared to GeForce 6800 with $30.2\mathrm{GB}/\mathrm{s}$, i.e., nearly six times faster.

Thus, these new advances in hardware development have been employed to implement the two-step MC simulation of photon packet propagation in scattering medium. In the first step, the MC code generates 24,576 trajectories satisfying the detection conditions. These trajectories are stored in the GPU’s global/page-locked memory and used in the second step to calculate the resulting intensity of scattered light, as presented in Fig. 3. The generation of photon trajectories repeats until the desired number of photons has reached each pixel of the detector (typically, $\sim 5000$ or more). These pixels are arranged a three-dimensional (3-D) dynamic array $[x\times y\times N_{\mathrm{ph}}]$ consisting of the data structures to mimic a real detector (e.g., CCD camera). The $x$ and $y$ are the indices of the array corresponding to an individual pixel (output coordinates of the photon packet trajectory) at the detecting area, $N_{\mathrm{ph}}$ is the total number of detected photon packets that defines the number of the data structures/trajectories preserved at the pixel for further calculation. For the detector $100\times 100\text{pixels}$ pixels with the pixel size 10 μm, the maximum memory required to store the array of ${10}^9$ structures (or ${10}^6$ to ${10}^7$ trajectories) is $\sim 50\text{GB}$ for the single FPP and $\sim 100\text{GB}$ for the double FPP. The current version of CUDA $(5·x)$ provides support only for 1-D arrays, so the standard multidimensional array linearization procedure has been employed to represent the 3-D array in the 1-D memory space.²²

Fig. 3

The block-diagram presents the principles of O3MC GPU-based computations of the resulting intensity of scattered light and its spatial distribution at the detector area. Each pixel at the detection area is associated with $(x,y)$ coordinates and containes $N_{\mathrm{ph}}$ photon trajectories.

Thus, calculation of the resulting intensity at each pixel is governed by a single GPU thread according to a particular problem and associated mathematical expression discussed below. This two-step MC approach permits the reuse of the same trajectories for calculating various quantities of light scattered in the medium (i.e., intensities of linearly and elliptically polarizations, intensity fluctuations, correlation of scattered light, etc.). This avoids the need to regenerate photon packets’ trajectories which is the most time-consuming part of MC modeling.

3.

Model of Coherent Polarized Light Propagation in Scattering Media

Polarization and coherence are the fundamental properties of light with respect to its wave nature.²⁵ Polarization is defined by the electric vector lying in the plane perpendicular to the direction of the light wave propagation. Coherence is a quantitative measure of the degree of phase correlation of the light wave. Various MC models have been developed in the past to simulate the propagation of polarized light in scattering media.^{12,14,15,26–29} To describe and track the polarization changes, recently developed MC approaches typically utilize Jones^30,31 and/or Stokes-Mueller^32,33 formalisms, but struggle with the limitations when modeling of coherent polarized light propagation in highly scattering media is required. In this paper, in the framework of further development of the unified MC model for the needs of biomedical optics,^20,21 we adapt the Jones-based formalism to handle linear or circular polarization of coherent light traveling through a random turbid medium.

3.1.

Linearly Polarized Light

In the Jones approach, linearly polarized light is presented as follows:³⁴

Eq. (1)

$$\overrightarrow{\mathbf{E}}=\left[\begin{array}{c}{E}_x\\ E_y\end{array}\right]=\left[\begin{array}{c}{E}_{0_x}{e}^{i{\varphi }_x}\\ E_{0_y}{e}^{i{\varphi }_y}\end{array}\right]=\left[\begin{array}{c}{E}_{\mathrm{p}}\\ E_{\mathrm{s}}\end{array}\right],$$

where $E_{0_x}$ and $E_{0_y}$ are the instantaneous $x$ and $y$ amplitudes of the wave propagating in the $z$-direction, and coincident with $\mathrm{p}$ and $\mathrm{s}$ polarizations ($E_{\mathrm{p}}$ and $E_{\mathrm{s}}$, respectively), $i$ is the imaginary unit, and $\varphi$ is the phase. The propagation of a wave through the medium is then described by a $2\times 2$ matrix:³³

Eq. (2)

$${\overrightarrow{\mathbf{E}}}^{\mathrm{out}}=\mathbf{J}·{\overrightarrow{\mathbf{E}}}^{\mathrm{in}},$$

where $\mathbf{J}$ is the Jones matrix, and ${\overrightarrow{\mathbf{E}}}^{\mathrm{in}}$ and ${\overrightarrow{\mathbf{E}}}^{\mathrm{out}}$ define, respectively, the polarization of incident beam and the beam emerging from a polarizing element that changes the state of polarization. If the polarized light propagates through multiple polarizing elements, the resulting Jones matrix is given by the product of each element (i.e., by multiple $2\times 2$ matrix multiplications³⁵). An application of the alternative Stokes-Mueller formalism for the same problem requires multiplication of $4\times 4$ matrices. The implementation of these approaches into MC modeling requires considerable allocation of memory and time. For instance, to store the 16-element Mueller matrices for 24,576 photon packets sampled in the scattering medium (considering the limited number of scattering orders, $N_{\mathrm{s}}={10}^3$) requires an allocation of $\sim 14\text{GB}$ of the GPU’s global memory. The modeling of ${10}^9$ photon packets would require an allocation of $\sim 600\text{PB}$ of global memory. To the authors’ knowledge, no currently available GPUs provide that much memory.

In the MC code, presented here, the propagation of linear polarized light in a scattering medium is modeled an analogy to the iterative procedure of the solution of Bethe-Salpeter equation.¹⁰ The incident light polarized along the $x$-axis is defined by the unit vector ${\overrightarrow{\mathbf{P}}}^{\mathrm{in}}=[P_x,P_y,P_z]=[1,0,0]$. The output vector ${\overrightarrow{\mathbf{P}}}^{\mathrm{out}}$ is determined at the end of each photon trajectory by the sequential transformations of the initial polarization:^10–11

$${\overrightarrow{\mathbf{P}}}^{\mathrm{out}}=\prod _{j=1}^{N_{\mathrm{s}}}[\widehat{I}-\frac{(r_{j+1}-r_j)\otimes (r_{j+1}-r_j)}{{|r_{j+1}-r_j|}^{-2}}]{\overrightarrow{\mathbf{P}}}^{(\mathrm{in})},$$

where $\widehat{I}$ is the unit forth-rank tensor and $\otimes$ indicates a direct product of two vectors (see Fig. 2). Finally, the corresponding Jones vector components $E_x$ and $E_y$ are defined as

Eq. (3)

$$\begin{gathered} E_x^2=W_j{P}_x^2{\mathrm{\Gamma }}_{\mathrm{R}}^{N_{\mathrm{s}}}, \\ {E}_y^2=W_j{P}_y^2{\mathrm{\Gamma }}_{\mathrm{R}}^{N_{\mathrm{s}}}, \end{gathered}$$

where $W_j$ is the statistical weight of $j$’th photon packet (the initial weight of the photon packet is $W_0=1$), and ${\mathrm{\Gamma }}_{\mathrm{R}}=2{(1+\overline{{\cos }^2\theta })}^{-1}$ is the Rayleigh factor. This relationship is based on the optical theorem,^36,37 which relates the scattering cross section and the scattering pathlength $l_{\mathrm{s}}$.

Total internal reflection and/or refraction at the medium boundary are taken into account by splitting the photon packet into the transmitted and reflected parts.³⁸ Thus, the statistical weight of the detected $j$’th photon packet that experienced $M$ number of reflections/refractions at the medium boundary is defined as

Eq. (4)

$$W_j=W_0[1-R^{\mathrm{in}}(\alpha )][\prod _{f=1}^{M_j}{R}_f(\alpha )][1-R^{\mathrm{out}}(\alpha )],$$

where $R^{\mathrm{in}}$ and $R^{\mathrm{out}}$ are the Fresnel reflection coefficients at the medium/air boundary for the photon packets entering and leaving the medium, respectively; $R(\alpha )$ is the Fresnel reflection coefficient that accounts for all other reflections at the surface of the medium. The Fresnel reflection coefficients (Fig. 4) are defined for $\mathrm{s}$ and $\mathrm{p}$ polarizations as³⁹

Eq. (5)

$$\begin{gathered} R_{\mathrm{s}}={\left|\frac{n_1\cos {\alpha }_i-n_2\cos {\alpha }_t}{n_1\cos {\alpha }_i+n_2\cos {\alpha }_t}\right|}^2, \\ {R}_{\mathrm{p}}={\left|\frac{n_1\cos {\alpha }_t-n_2\cos {\alpha }_i}{n_1\cos {\alpha }_t+n_2\cos {\alpha }_i}\right|}^2. \end{gathered}$$

Here, ${\alpha }_i$ and ${\alpha }_t$ are the angles of incidence and transmittance of the photon packet at the surface of the medium, respectively, $n_1$ and $n_2$ are the refractive indices.

Fig. 4

Fresnel reflection coefficients $R_{\mathrm{s}}$ (circles) and $R_{\mathrm{p}}$ (triangles) for light traveling: (a) from air ($n_1=1$) to the medium ($n_2=1.5$), and (b) from the medium ($n_1=1.5$) to air ($n_2=1$).

Absorption is taken into account according to the microscopic Beer-Lambert law³⁸ by attenuating the statistical weight $W_j$ of the photon packet by the factor $\exp (-{\mu }_{\mathrm{a}}L)$, where ${\mu }_{\mathrm{a}}$ is the absorption coefficient and $L$ is the photon packet pathlength.

Finally, for the detected photon packets (typically, $N_{\mathrm{ph}}={10}^7-{10}^9$), taking into account the interference and the coherent properties of light, the resulting intensities of polarized $I_{\parallel }$ and depolarized $I_{\perp }$ light are defined as

Eq. (6)

$$\begin{gathered} I_{\parallel }=\frac{1}{N_{\mathrm{ph}}}\sum _{i=1}^{N_{\mathrm{ph}}}{E}_{x_i}^2+\frac{2}{N_{\mathrm{ph}}}\sum _{i=1}^{N_{\mathrm{ph}}}\sum _{j=i+1}^{N_{\mathrm{ph}}}{E}_{x_i}{E}_{x_j}\times \cos [\frac{2\pi }{\lambda }(L_i-L_j)+\pi (m_i-m_j)]\exp [-{\left(\frac{L_i-L_j}{l_{\mathrm{c}}}\right)}^2], \\ {I}_{\perp }=\frac{1}{N_{\mathrm{ph}}}\sum _{i=1}^{N_{\mathrm{ph}}}{E}_{y_i}^2+\frac{2}{N_{\mathrm{ph}}}\sum _{i=1}^{N_{\mathrm{ph}}}\sum _{j=i+1}^{N_{\mathrm{ph}}}{E}_{y_i}{E}_{y_j}\times \cos [\frac{2\pi }{\lambda }(L_i-L_j)]\exp [-{\left(\frac{L_i-L_j}{l_{\mathrm{c}}}\right)}^2], \end{gathered}$$

where $L_i$ and $L_j$ are the individual total pathlengths of $i$’th and $j$’th photon packets, respectively; $m$ is the number of displacements of $E_x$ giving a correct representation of the nodality⁴⁰ for the $i$’th and $j$’th photon packets, $l_{\mathrm{c}}$ is the temporal coherence length of the incident radiation, and $\lambda$ is the wavelength of the incident light.

3.2.

Circularly Polarized Light

A similar approach can be applied for modeling of circularly and/or elliptically polarized light propagation in turbid scattering media. Circularly polarized light is characterized by two orthogonal vectors, i.e., for the incident light with right-handed circular polarization (RCP) each photon packet injected into the medium is assigned with two orthogonal vectors ${\overrightarrow{\mathbf{P}}}_x^{\mathrm{in}}=[1,0,0]$ and ${\overrightarrow{\mathbf{P}}}_y^{\mathrm{in}}=[0,1,0]$ with a phase shift $({\varphi }_x-{\varphi }_y=\pi /2)$. Thus, the resulting components of RCP and left-handed circular polarization (LCP) at the end of each photon trajectory are, respectively, defined as follows:

Eq. (7)

$$\begin{gathered} \left[\begin{array}{c}{E}_{xx}^2\\ E_{yy}^2\end{array}\right]=\left[\begin{array}{c}{W}_j{P}_{x_x}^2{\mathrm{\Gamma }}_{\mathrm{R}}^{N_{\mathrm{s}}}\\ W_j{P}_{y_y}^2{\mathrm{\Gamma }}_{\mathrm{R}}^{N_{\mathrm{s}}}\end{array}\right], \\ \left[\begin{array}{c}{E}_{xy}^2\\ E_{yx}^2\end{array}\right]=\left[\begin{array}{c}{W}_j{P}_{x_y}^2{\mathrm{\Gamma }}_{\mathrm{R}}^{N_{\mathrm{s}}}\\ W_j{P}_{y_x}^2{\mathrm{\Gamma }}_{\mathrm{R}}^{N_{\mathrm{s}}}\end{array}\right]. \end{gathered}$$

It should also be pointed out here that when circularly or elliptically polarized light is reflected at the medium surface, the state of polarization undergoes a helicity flip.⁴¹ The same goes for backscattered light (Fig. 5); i.e., if the photon packet is scattered backwards [where $\theta >\pi /2$, $\theta$ is the angle between the incident $k$ and scattered $k^{\prime }$ wave vectors ($\theta ={\widehat{kk}}^{\prime }$), see Fig. 5], this results in an RCP photon packet changing to LCP; if the photon packet scatters forward ($\theta <\pi /2$), the helicity remains the same. This helicity flip phenomenon is of fundamental importance.⁴² Linear polarization possesses no such sense of the direction in which it travels.

Fig. 5

Schematic presentation of forward and backscattering for circularly polarized light. When right-handed circular polarization (RCP) light is scattered in forward direction (${\widehat{kk}}^{\prime }<\pi /2$) its helicity is preserved, whereas for backscattered light (${\widehat{kk}}^{\prime }>\pi /2$), the phase shift $\delta ={\varphi }_x-{\varphi }_y$ between $E_x$ and $E_y$ is reversed, causing the light to change from RCP to left-handed circular polarization (LCP).

Therefore, considering the propagation of circularly polarized light through a turbid medium where multiple scattering events occur, the photon packets that have backscattered an odd number of times will correspond to a reversal in helicity, and thus, contribute the cross-polarized portion of the detected signal ($I_{\mathrm{cross}}$). Correspondingly, the photon packets that have experienced an even number of backscattering events contribute to the copolarized signal ($I_{\mathrm{co}}$) and the handedness of light is unchanged with respect to the incident polarization. Thus, for $m$ total number of helicity flips, either due to reflections on the medium boundary or backscattering events ($\theta >\pi /2$), and taking into account interference, the copolarized light (that has preserved its helicity) is given by

Eq. (8)

$$\begin{gathered} I_{xx}=\frac{1}{N_{\mathrm{ph}}}\sum _{i=1}^{N_{\mathrm{ph}}}(E_{x{x}_i}^2\left|\cos \frac{\pi m_i}{2}\right|+E_{y{x}_i}^2\left|\sin \frac{\pi m_i}{2}\right|) \\ +\frac{2}{N_{\mathrm{ph}}}\sum _{i=1}^{N_{\mathrm{ph}}}\sum _{j=i+1}^{N_{\mathrm{ph}}}(E_{x{x}_i}\left|\cos \frac{\pi m_i}{2}\right|+E_{y{x}_i}\left|\sin \frac{\pi m_i}{2}\right|)(E_{x{x}_j}\left|\cos \frac{\pi m_j}{2}\right|+E_{y{x}_j}\left|\sin \frac{\pi m_j}{2}\right|) \\ \times \cos [\frac{2\pi }{\lambda }(L_i-L_j)+\pi (m_i-m_j)]\exp [-{\left(\frac{L_i-L_j}{l_{\mathrm{c}}}\right)}^2], \\ {I}_{yy}=\frac{1}{N_{\mathrm{ph}}}\sum _{i=1}^{N_{\mathrm{ph}}}(E_{y{y}_i}^2\left|\cos \frac{\pi m_i}{2}\right|+E_{x{y}_i}^2\left|\sin \frac{\pi m_i}{2}\right|) \\ +\frac{2}{N_{\mathrm{ph}}}\sum _{i=1}^{N_{\mathrm{ph}}}\sum _{j=i+1}^{N_{\mathrm{ph}}}(E_{y{y}_i}\left|\cos \frac{\pi m_i}{2}\right|+E_{x{y}_i}\left|\sin \frac{\pi m_i}{2}\right|)(E_{y{y}_j}\left|\cos \frac{\pi m_j}{2}\right|+E_{x{y}_j}\left|\sin \frac{\pi m_j}{2}\right|) \\ \times \cos [\frac{2\pi }{\lambda }(L_i-L_j)+\pi (m_i-m_j)]\exp [-{\left(\frac{L_i-L_j}{l_{\mathrm{c}}}\right)}^2]. \end{gathered}$$

Similar equations can be written for depolarized light, defined by $I_{yx}$ and $I_{xy}$, respectively.

4.

Results and Discussion

Figure 5 presents the results of simulation of spatial distribution of $I_{\parallel }$, $I_{\perp }$, and the depolarization ratio $\mathrm{DR}=(I_{\parallel }-I_{\perp })/(I_{\parallel }+I_{\perp })$ at the surface of the medium with various values of scattering anisotropy. These results both quantitatively and qualitatively are well agreed with the results of semianalytic MC model,¹⁵ a comparison of two MC models presented in Ref. 43.

It can be clearly seen from Fig. 6 that $I_{\parallel }$ significantly prevails over the intensity of the depolarized light $I_{\perp }$. These results are in good agreement with the results of an analytical solution and other models (Table 1). The exact Milne solution for the Rayleigh $(g=0)$ noncoherent $(l_{\mathrm{c}}=0)$ scattering gives $I_{\parallel }/I_{\perp }\approx 1.92$^44–47 and $\mathrm{DR}=0.31$,^44–47 while the results of MC modeling suggests $I_{\parallel }/I_{\perp }=1.93$ and $\mathrm{DR}=0.3172$; a similar value of $\mathrm{DR}=0.33$ has been obtained elsewhere.¹²

Fig. 6

Surface distributions of $I_{\parallel }$, $I_{\perp }$, and depolarization ratios of backscattered linearly polarized light for various values of scattering anisotropy: from the top $g=0$, 0.5, 0.9 and 0.98. MC simulations were performed with noncoherent light ($l_{\mathrm{c}}=0\mathrm{mm}$) normally incident on the semi-infinite scattering medium with ${\mu }_{\mathrm{s}}=30{\mathrm{mm}}^{-1}$, ${\mu }_{\mathrm{a}}=0.001{\mathrm{mm}}^{-1}$, and detected over entire medium surface.

Table 1

Comparison of the results of MC simulation with the exact solution by Milne and the results of alternative modeling (Refs. 44–47) Here, I∥ and I⊥ are counted as total intensity of I∥(x,y) and I⊥(x,y) profiles presented in Fig. 6.

Correspondingly, Fig. 7 presents the effective probing depth distribution (sampling volume) for $I_{\parallel }$, $I_{\perp }$, and their difference, counted by analogy to the depolarization ratio. As one can see for the medium with high anisotropy of scattering ($g=0.98$) the effective probing depth of the linearly polarized component $I_{\parallel }$ exceeds the propagation depth for the depolarized one $I_{\perp }$. Therefore, with a significant increase of the photon packets’ pathlengths within the medium, the residual polarization in the backscattered light becomes preserved due to a large relative contribution of the photon packets’ trajectories with low orders of scattering (see Fig. 6).

Fig. 7

Corresponding sampling volumes for $I_{\parallel }$ and $I_{\perp }$, and the spatial distributions of the depolarization ratio DR simulated for the same parameters as in Fig. 5.

Figure 8 presents the surface distribution of $I_{\parallel }$, $I_{\perp }$, and $\mathrm{DR}$ for backscattered linearly polarized light with various coherence lengths. The mutual interference between individual photon packets results in the formation of speckle patterns.

Fig. 8

Intensity distributions of $I_{\parallel ,\perp }$ and the depolarization ratios DR at the surface of the medium for backscattered linearly polarized light of various coherence lengths $l_{\mathrm{c}}$, from top: $l_{\mathrm{c}}=0.01\mu \text{m}$, 0.1 mm, 1 mm, 1 m, counted for a highly anisotropic scattering medium with $g=0.98$, ${\mu }_{\mathrm{s}}=30{\mathrm{mm}}^{-1}$, ${\mu }_{\mathrm{a}}=0.001{\mathrm{mm}}^{-1}$.

Figure 9 shows the surface distribution of the portions of copolarized $I_{\mathrm{co}}$ (RCP) and cross-polarized $I_{\mathrm{cross}}$ (LCP) intensities of backscattered light, and the helicity flip is clearly observed.

Fig. 9

Surface distributions of $I_{\mathrm{co}}$ (top) and $I_{\mathrm{cross}}$ (bottom) polarized components of backscattered circularly polarized light, $d$ is the distance at the surface, $l^{*}=1/{\mu }_{\mathrm{s}}^{\prime }$; counted for a semi-infinite scattering medium (${\mu }_{\mathrm{s}}=30{\mathrm{mm}}^{-1}$, ${\mu }_{\mathrm{a}}=0.001{\mathrm{mm}}^{-1}$, $g=0.98$, and $n=1.5$).

The results of simulation are in a good agreement with the analytical results and the results of alternative modeling; see Table 2 for details.

Table 2

Comparison of the results of MC simulation with the results of Milne’s solution and the results of alternative modeling (Refs. 44–47).

The model has been compiled, tested, and proven to work on Windows 8.1/Ubuntu GNU Linux 13.10 platforms. The time required to simulate ${10}^{11}$ photon packet trajectories is dependent on the parameters of the scattering medium and source-detector configuration. For instance, computing the surface distribution for the abovementioned detector parameters for a semi-infinite scattering medium (${\mu }_{\mathrm{s}}=30{\mathrm{mm}}^{-1}$, $g=0.98$, $n=1.5$ and ${\mu }_{\mathrm{a}}=0.001{\mathrm{mm}}^{-1}$), illuminated by noncoherent light with $l_{\mathrm{c}}=0$, takes $\approx 2\mathrm{h}$, and $\approx 8\mathrm{h}$ for coherent polarized light ($l_{\mathrm{c}}\ge 0.01\mathrm{mm}$) with the same medium parameters.

5.

Conclusions

In the current report, we have introduced the electric field MC approach specifically developed for fast modeling of linearly and circularly polarized light propagation in highly scattering media. The model takes into account the coherent properties of light, the influence of total reflection and refraction at the medium boundary, helicity flip of circularly polarized light and interference. The results of modeling for both linearly and circularly polarized light are in good agreement with the results of the exact analytical solution by Milne and the results of alternative modeling. The developed MC approach is extremely flexible to operate with different optical properties of the medium, including scattering, absorption, refractive index, anisotropy of scattering, etc. The presented MC approach is a part of the O3MC model²⁰ and can be easily extended to model light propagation in optically active and/or birefringent scattering media, as well as for scattering of complex vector laser beams in turbid tissue-like media.