1. Introduction

Surface tissue layers are the regions where pre-cancerous alterations occur and thus provide the most important diagnostic information for early cancer detection. Optical techniques specifically probing the surface layers have been under intensive development for over a decade. Among these, polarization gating [2,3] and low-coherence enhanced backscattering (LEBS) [4,5] are two powerful means to selectively detect the signals from the superficial tissue and reject those from the underlying bulks. A combination of polarization and EBS would reveal even more properties of the non-diffusive layer just under the surface. The incident polarization predefines a specific direction of space, so certain spatial symmetry, such as the cylindrical symmetry, is broken. The exposed spatial/angular variation of EBS may carry extra information about the subject under investigation. Experimental data employing linear polarization demonstrate that the co-polarized EBS cone for Rayleigh scatterers possess different widths at two orthogonal detection planes, while such difference is absent for Mie scatterers [6,7].

EBS is a manifestation of time-reversal symmetry in light scattering in turbid medium, where the photons traveling along exact opposite paths interfere constructively in the backward direction. Typical theoretical modeling of EBS utilizes a Fourier Transform on incoherent backscattering signal calculated either using diffusion approximation [8] or Monte Carlo [9,10]. Such approach is valid for scalar waves, but cannot handle the depolarization process for polarized light. Recent development of electric field Monte Carlo method tracks the electric fields of photons along time-reversed paths and explicitly calculates their interference [11]. The simulation conducted on a turbid medium containing 100 nm diameter polystyrene spheres under linearly polarized illumination at 514.5 nm wavelength predicts different patterns for the co-polarized and cross-polarized EBS cones. Moreover, numerical analysis using superposition T-matrix method also reproduces the azimuthal asymmetry in co-polarized EBS cone for Rayleigh scatterer clusters even at very high volume concentrations [12].

However, the azimuthal isotropy observed in turbid medium containing Mie scatterers still needs to be explained [6,7]. Finite difference time domain (FDTD) [13] and its variation pseudo-spectral time domain (PSTD) [14,15] can provide exact numerical solutions to Maxwell’s equation of arbitrary material configuration. Further, recent development of PSTD on staggered grid and local Fourier basis (SLPSTD) can achieve full parallelization of PSTD in a distributed supercomputing environment and enables simulation of light scattering at tissue level [1]. Previous analysis of polarization-independent EBS disclosed that the ultra-wide LEBS cone originates from the low-order scattered photons from the non-diffusive layer [1,5,16]. In this paper, a further study is conducted on the EBS cone of Mie scatterer clusters under linearly polarized incidence, in which the backscattered signal is collected in either co-polarized or cross-polarized configuration. We discover both EBS cones from the superficial layer (corresponding to low-order backscattering) exhibit azimuthal anisotropy patterns which are extremely sensitive to the photons’ penetration depth. These patterns quickly blur out as the order of scatterings increase.

Yet another prominent polarization effect associated with EBS is the polarization memory effect in Mie scattering clusters when circular polarization incidence is employed [17]. Unlike the linear polarization which totally depolarizes after several scattering events, circularly polarization can survive much longer pathlengths. In this paper, we present an exact numerical investigation of this effect using SLPSTD. The curve of the difference signal between the co-polarized and cross-polarized backscatterings versus the penetration depth is obtained. This curve has a dip at small depth, indicating a helicity-flipping effect associated with top layer backscattering; and then rises linearly as the photons propagate deeper, indicating a lack of depolarization as the difference signal does not attenuate. The signal from the deep layers consists in photons undergone multiple small-angle scatterings which largely preserves their helicity states. From the curve the thickness of the helicity-flipping layer can be identified, and the exact boundary between the helicity flipping and helicity preserving regions is located. Our simulation method can be applied to evaluate the performance of subsurface imaging techniques that explore the circular polarization memory effect [18–20].

2. Brief description of the SLPSTD method

SLPSTD is a renovation over traditional PSTD. Details of SLPSTD can be found elsewhere [1]. Major new features of SLPSTD include Fourier transforms on local basis, allocation of field positions on Yee lattice and overlapping domain decomposition of the model space. The outcome of these modifications is high-efficiency parallelization on memory-distributed computer clusters. In this algorithm the spatial derivative of field ψ in Maxwell equations is reformatted as ${\psi }^{\prime }\equiv F^{-1}\left[jk{e}^{jk\mathrm{\Delta }/2}F\left(\psi \right)\right]$, where k is the discrete wave vector, Δ the grid spacing, and F and $F^{-1}$ the FFT and inverse FFT, respectively. This revision of the derivative operator allows the locations of electric fields and magnetic fields to be shifted by half a grid, so a Yee lattice can be adopted. The staggered-gird configuration can naturally handle spatial discontinuities, such as point sources, in the computation space without resorting to point-spreading tricks [21,22].

In SLPSTD, the whole computation domain is decomposed into subdomains in 3D space, which are one-to-one mapped to different nodes in a parallel computing environment. Inter-node data exchange is performed on the overlapping region between adjacent subdomains. The exchanged data are then processed by an artificial taper function before adding to the local data. FFT is conducted on the local data within each node. This procedure ensures the continuity of the derivatives across subdomain border.

In our previous work, monochromatic continuous wave was used in each run, so multi-wavelength calculations would require multi-runs [1]. In this article, the code has been revised to enable pulsed wave simulations that can finish multi-wavelength calculations in a single run. Accordingly the “on the fly” procedure is adopted to recode the near-to-far field transformation [13].

3. Simulations and discussions

Polarization dependent nondiffusive layer backscattering is simulated using the SLPSTD method. The random medium in the modeling is a water suspension of 2μm diameter polystyrene spheres at 2.9% volume concentration. The refractive indices are 1.33 and 1.59 for water and polystyrene, respectively. The shape of the medium is rectangular with the horizontal dimensions fixed at $100\mu m\times 100\mu m$, while its thickness is subject to variation to control the order of scatterings. Please note the medium is embedded in water, not air. A quasi-monochromatic laser pulse, either linearly or circularly polarized, is incident upon the medium’s top surface. The central frequency $f_0$ of the pulse is $3.82\times {10}^{14}$Hz, corresponding to 785nm wavelength in vacuum. Mie theory predicts that at $f_0$ the scattering mean free path ${\ell }_s=12.5\mu m$ and the anisotropy factor $g=0.93$. The simulation is repeated on a series of medium thickness d ranging from 12.5μm to 75μm, corresponding to optical thickness τ ($\equiv d/{\ell }_s$) from 1 to 6. The numerical realization of the medium is achieved by placing the microspheres randomly yet uniformly in the rectangle. To suppress speckles, the calculated backscattering signal is averaged at two stages. The first one involves an average over 60 frequencies evenly distributed in $\left[0.95f_0,1.05f_0\right]$ for each realization of the medium [23]. Here the multi-frequency data are resolved from the pulsed wave simulation in a single run. The first stage is repeated on 5 different medium realizations and a further mean of the 5 results produces the final value of the backscattering signal.

Unlike the backscattering from deeper diffusion layers where the information about the initial incidence polarization is totally lost, non-diffusive layer backscattering exhibits polarization dependent patterns. Azimuth anisotropy under linearly polarized incidence and polarization memory effect under circularly polarized incidence are investigated in the following. In this study, the order of scatterings experienced by the detected photons is tuned by varying the optical thickness $\tau ={\mu }_{s}d$, i.e., the dimensionless product of the scattering coefficient and the layer thickness. Therefore, an injected photon would in average undergo τ number of scatterings before propagating vertically into depth d. In our simulations the thickness of the superficial layers satisfies$\tau \le 6$, while ${(1-g)}^{-1}\approx 14$ scatterings are required for the photon to become diffusive. So the scatterings in the lateral directions are insignificant. In fact, τ is a commonly used parameter to characterize scattering orders in low-order backscattering studies [1–3].

3.1 Azimuth anisotropy in nondiffusive EBS under linearly polarized incidence on Mie scatterer clusters

The linearly polarized plane wave is incident on the medium’s top surface at a small angle of $5^{°}$ tilted away from the normal to avoid the specular reflection (Fig. 1 ). The backscattering angle θ and the azimuthal angle φ are defined in the observer’s coordinate system $x^{\prime }y{z}^{\prime }$ which is a $5^{°}$ rotation of the medium coordinate system $xyz$ along the y axis. The incidence polarization and the backscattering components (co-polarized $I_{\parallel }$ and cross-polarized $I_{\perp }$) are also shown in Fig. 1. The $z^{\prime }$ axis would be the exact backward direction and the center of the EBS cone. Figures 2a and 2b present the co-polarized and cross-polarized EBS, respectively. Please note Figs. 2a and 2b have different scale of false color. Within each set of plots, the upper row corresponds to $\tau =1,2,3$ and the lower row $\tau =4,5,6$ from left to right. Each color map represents a view of the EBS intensity around the exact backward direction. The radial coordinate is the backscattering angle θ, ranging from ${\text{0}}^{°}$at the center to ${\text{8}}^{°}$at the border. The polar angle is the azimuth angle φ of the backscattering direction, with $\phi =0^{°}$defined in the direction of the incident polarization. Azimuthal anisotropy is evident in both co-polarized and cross-polarized EBS for $\tau \le 3$, but progressively blurs out as τ increases. The co-polarized EBS cones extend in direction parallel to the incidence polarization; while the cross-polarized ones present a fourfold “X” pattern.

 

Fig. 1 Schematic diagram of the system setup. To avoid the specular reflection, the incident light is slightly tilted from the normal of the medium surface by a $5°$ angle, so the coordinate system $x^{\prime }y{z}^{\prime }$ is a rotation of the coordinate system $xyz$ along the y axis by $5°$. The backscattering angle θ and the azimuthal angle φ are defined in the $x^{\prime }y{z}^{\prime }$ coordinate system, where the $z^{\prime }$ axis is the exact backward direction and the center of the EBS cone. The incidence polarization lies in the $x^{\prime }{z}^{\prime }$ plane, whereas the polarization of the backscattered light is decomposed into a component parallel to $x^{\prime }{z}^{\prime }$ plane ($I_{\parallel }$) and its orthogonal component ($I_{\perp }$).

Download Full Size | PDF

 

Fig. 2 (a) Co-polarized and (b) cross-polarized EBS from random medium consisting of Mie scatterers under linearly polarized incidence. Each color map represents a view of the EBS intensity along the exact backward direction. In both (a) and (b), from left to right, the upper row corresponds to medium of optical thickness 1, 2 and 3, and the lower row 4, 5 and 6, respectively. The backscattering angle θ ranges from $0^{°}$at the center to $8^{°}$at the border. The azimuth angleφis defined from the direction of the incident polarization ($\phi =0°$). Azimuthal anisotropy appears as intensity fluctuation within each color map. The patterns are the most prominent for $\tau =1$ in both (a) and (b), but progressively become blurred as τ increases to 6.

Download Full Size | PDF

To characterize the degree of anisotropy, the backscattering intensity is integrated over θ at each azimuth angle φ as follows

 

Fig. 3 Polar diagrams of the integrated intensities (Eqs. (1)-(2), blue line) and their empirical fitting (Eqs. (3)-(4), red line): (a) $\tau =1$, co-polarized; (b) $\tau =6$, co-polarized; (c) $\tau =1$, cross-polarized; and (d) $\tau =6$, cross-polarized. Because only the curve’s shape is relevant to the degree of anisotropy, each individual curve is normalized to its own intensity at $\phi =0°$.

Download Full Size | PDF

 

Fig. 4 Parameters characterizing the degree of anisotropy: (a) ε for the co-polarized EBS, and (b) γ for the cross-polarized EBS.

Download Full Size | PDF

Actual experiments report an anisotropy existing in the bulk EBS from Rayleigh scattering medium, but non-observable in the case of Mie scattering medium [6]. Such discrepancy can now be explained as follows based on the results of Fig. 3. The low-order EBS from the non-diffusive superficial layer is indeed anisotropic for both Rayleigh and Mie scatterings, but its percentages in the total EBS are different. In the case of small Rayleigh scatterers, the g factor is small and large-angle scattering is significant. Low-order scatterings are more capable to turn the incident photons’ direction backward and make a large contribution to the total signal. On the other hand, each Mie scattering is mainly forward and multiple scatterings are required to turn the incident photon’s direction backward, and inevitably the polarization state is lost along the long path. Consequently, the bulk EBS is dominated by the isotropic signal from deep layers, in which the anisotropic low-order EBS component is buried.

3.2 Polarization memory effect in nondiffusive layer backscattering under circularly polarized incidence

In this section, nondiffusive layer backscattering under circularly polarized incidence is simulated. Illumination of right circular polarization is incident perpendicularly upon the top surface of the turbid medium. Because we intent to compare the absolute intensities of backscatterings from media of various thicknesses (Fig. 5 ), a tilted incidence is inconvenient because the varying sidewall would change the illuminated area. At normal incidence, the specular reflection is determined to be within a tiny angle ${\theta }_0=0.25°$ which is then excluded in all the following analysis.

 

Fig. 5 Difference signal between co-polarized and cross-polarized backscatterings for various media thickness. The region $\theta <{0.25}^{°}$is excluded to remove the specular reflection. The curves are found to be equidistant when$\tau >2$, indicating the two helicity contents cannot be equalized as the photons propagate deeper and the incident circular polarization state is remembered. On the other hand, the signal is negative when $\tau <3$, which means the overall backscattering from this region is helicity-flipped.

Download Full Size | PDF

The most prominent feature involving circular polarization is the polarization memory effect [17,24]. In previous SLPSTD simulation of linear-polarization gating [1], the difference signal between the co-polarized and cross-polarized backscatterings was found to rise to a saturation plateau at $\tau =4$ and consequently the gated region is $\tau \le 4$. Here the same procedure can be applied to investigate the circular polarization scattering effect. The difference signal between co-polarized (helicity preserving) and cross-polarized (helicity flipping) backscatterings is calculated for random media of various optical thicknesses as

Helicity flipping shows up in Fig. 5 as the negative difference signal when $\tau <3$. In this region, which is also associated with low-order backscattering, the cross-polarized signal exceeds the co-polarized one. To determine the exact depth of helicity-flipping, we conducted extra simulations on refined optical thicknesses $\tau =$0.6, 0.8, 1.2, 1.4 and 1.6. The difference signal (Eq. (7)) is further integrated, i.e.

The integrated intensity is plotted in Fig. 6 . The data points within $0.6<\tau <2$ are easily fitted to a third polynomial and the minimum of the fitting is located at optical depth 0.95. Therefore, we can conclude that the helicity-flipping is mainly due to the single large-angle scattering within $\tau \le 0.95$. In region $2<\tau <6$ the data points fall onto a straight line which crosses the 0-intensity line at $\tau =2.4$, indicating the overall signal from regions above this plane still favors helicity flipping.

 

Fig. 6 Integrated intensity vs. optical thicknesses, (a) the $0.6<\tau <2$ portion is a third order polynomial fit to emphasize the helicity flipping; (b) the $2<\tau <6$ portion is a linear fit to reflect the polarization memory effect. In (a) the minimum is determined at $\tau =0.95$. The descending part ($\tau <0.95$) indicates the first single large-angle scattering is responsible for helicity flipping; the curve is still negative for $0.95<\tau <2.\text{4}$, indicating the overall effect still favors helicity flipping.

Download Full Size | PDF

Existence of the circular polarization memory effect and a $\tau =0.95$ helicity-flipping top layer can be explained as follows: in Mie scattering, the dominant forward scattering preserves helicity while the large angle scattering flips it; the signal component from deep layers consists in photons undergone multiple small angle scatterings and the helicity information is remembered along the long path. The helicity-flipping component originates mostly from the single large-angle backscattering within the layer ~1 optical thickness beneath the surface, but the rising curve within $1<\tau <2$ indicates the competition from multiple small angle scatterings starts to build up.

4. Improvement on code execution efficiency

All the computations were performed on the Dawning 5000A supercomputer at the Shanghai Supercomputer Center (SSC). The supercomputing environment remains the same as Ref [1], including 8 hardware nodes connected via Infiniband. Each hardware node contains four 2.0GHz 4-core AMD Barcelona CPUs, resulting in 32 CPUs with totally 128 cores. However, in this work pulsed excitation is incorporated into the code and simulations at multi-frequencies can be conducted in a single run. This greatly facilitates the suppression of speckles. For a direct comparison, we consider a medium of size $100\mu m\times 100\mu m\times 50\mu m$. The monochromatic wave approach in the previous work takes 5.9 hours per frequency [1], so an average over 60 frequencies would require 60 runs, or 354 hours. On the other hand, the 60-frequency simulations can be finished in a single run within 9.7 hours using the pulsed excitation implementation in this paper.

5. Summary

In this paper, polarization dependent characteristics of nondiffusive layer backscattering at Mie particle scale is investigated using the SLPSTD method. The progression of depolarization in the scattering process is simulated by varying the optical thickness of the medium. For linear polarization incidence, the EBS from the nondiffusive layer shows azimuthal anisotropy both for co-polarized and cross-polarized detection. The degree of anisotropy can be characterized by two model parameters and is found to decrease as the optical thickness increases. For Mie scatterers, the low-order component of EBS from the surface layer is in fact azimuthally anisotropic, but is buried in the isotropic component from the deep layers. The competition between the two components results in different EBS cone shapes observed in experiments for Rayleigh and Mie scatterer clusters. An EBS technique capable of resolving penetration depth, such as LEBS [4,5], should be able to unveil the azimuthal anisotropy in Mie scattering medium. For circular polarization studies, details about the polarization memory effect and the helicity-flipping phenomenon are for the first time disclosed by exact numerical solutions. A thin layer ~1 optical thickness under the surface of the Mie scattering medium is responsible for the helicity-flipping. Deeper than $\tau =2$, the difference signal between the co-polarized and cross-polarized backscatterings shows a linear dependence on the penetration depth. In other words, beyond the helicity-flipping region, the increased amount of signal is proportional to the increased penetration depth and thus the increased number of scatterers. Therefore, the helicity state is preserved as the photons propagate into successively deeper layers. Our numerical technique can be directly applied to study the behavior of circular polarization in turbid medium of complex geometry. The simulations can serve as guidelines for subsurface imaging that explores the circular polarization memory effect.