Propagation of Cylindrical Vector Laser Beams in Turbid Tissue-Like Scattering Media

: We explore the propagation of the cylindrical vector beams (CVB) in turbid tissue-like scattering medium in comparison with the conventional Gaussian laser beam. The study of propagation of CVB and Gaussian laser beams in the medium is performed utilizing the uniﬁed electric ﬁeld Monte Carlo model. The implemented Monte Carlo model is a part of a generalized on-line computational tool and utilizes parallel computing, executed on the NVIDIA Graphics Processing Units (GPUs) supporting Compute Uniﬁed Device Architecture (CUDA). Using extensive computational studies, we demonstrate that after propagation through the turbid tissue-like scattering medium, the degree of fringe contrast for CVB becomes at least twice higher in comparison to the conventional linearly polarized Gaussian beam. The results of simulations agree with the results of experimental studies. Both experimental and theoretical results suggest that there is a high potential of the application of CVB in the diagnosis of biological tissues.


Introduction
For a number of years, considerable attention has been given to the probing of biological tissues with photonic-based imaging and optical diagnostic techniques [1]. The interdisciplinary field of biomedical optical diagnosis aims to use light-based tools for non-invasive assessment of structural alterations in biological tissues and changes in their physiological and subsurface scattering properties (e.g., due to cancers, wounds, specific conditions, etc.) [1]. Photonic-based modalities provide strong advantages such as the ability to perform diagnostic procedures non-invasively on bio-tissues in vivo. Optical diagnostic systems utilize different types of illumination of biological tissues with visible, UV or infrared light, and attempt to analyze detected reflected, scattered, emitted or otherwise modulated optical radiation [2]. A number of prominent methods based on optical spectroscopy to measure tissue properties such as oxygenation and blood content have been developed since the 1930s [1]. In the seminal works, measurement of the fluorescence spectra has been utilized to differentiate between cancerous and normal tissues [3,4]. Subsequently, Raman spectroscopy was introduced as a robust detection tool for cancers in breast tissues [5] forming the basis for the so-called Raman 'optical biopsy' [6] that later underwent several considerable practical improvements [7][8][9].
On the other hand, polarization or spin angular momentum (SAM) are the light's quite salient features along with the spectral or coherence properties. Multiple aspects of polarized light scattering and interaction with biological tissues, as well as other turbid media have been extensively studied in the past [10][11][12][13]. A large proportion of these studies associated with polarized light are focused on non-invasive characterization of biological tissues with the special attention to earlier detection of malignancies [14]. Recently, it has been demonstrated that polarization and SAM are able to distinguish different grades of cancer [15].
In fact, the structure of light can be more complex, i.e., light beams can be radially or azimuthally polarized, related to their spatial structure (shown in Figure 1) [16]. It has been demonstrated, that arbitrary vector vortex beams can be generated and detected at longer distances and show substantially lower scintillation comparing to traditional scalar beams while propagating in scattering environments [17,18]. The non-separability of Vector beams were suggested to be used to encode information in optical communication [19]. In recent studies, the propagation of complex vector laser beams (known as the Laguerre-Gaussian (LG) beams) in mouse brain tissue was investigated demonstrating varying transmission in the different functional areas of the brain [20]. A higher-order Poincare sphere (HOPS) concept was specifically introduced for quantitative assessment of the Stokes vectors distribution [21,22]. Owing to the fact that the transformation of polarization of probing light during its transfer through the turbid tissue-like scattering medium is independent of the state of polarization of incident light [23], further improvements of the SAM-based techniques are required, especially for screening of superficial biological tissues. The development and practical implementation of complex vector laser beams, such as CVB and LG beams, in diagnostic practice, as well as interpretation and quantitative analysis of the experimental results, require extensive modelling in order to provide an ultimate understanding of their propagation in a random tissue-like scattering medium. Due to the complex and highly inhomogeneous multi-layered structure of biological tissues, no practical analytical solutions describing radiative transfer therein exists. Fortunately, a stochastic alternative, such as the Monte Carlo (MC) method, is well suited for simulation of the effects of multiple scattering and polarized light propagation in turbid tissue-like scattering media [24]. Numerous MC models have been developed in the past to simulate the propagation of polarized light [25][26][27][28][29][30][31][32].
In the current paper, within the framework of a further development of a unified MC modeling platform [33] we present a computational model for imitation of CVB's propagation through a turbid tissue-like scattering media with a primary aim to prove the concept of applicability of CVB in biomedical optical diagnosis. The MC model utilizes an analogue to Jones formalism, specially adapted for the online object oriented MC approach developed earlier [33,34]. The model takes into account the wave properties of light, including temporal coherence, polarization memory [35], interference and reflection/refraction at the medium boundary. To access quantitatively the degree of randomization of the probing light along its scattering trajectory within the medium the conventional metrics, known as depolarization ratio (DR) is utilized.

Online GPU-Accelerated Monte Carlo Modeling of Photon Migration in Scattering Media
MC modeling of energy transfer through the medium has been widely studied and its principles are described in detail elsewhere [36]. Several practical implementations of the MC approach for simulation of light propagation in turbid tissue-like scattering medium have been developed [28,[37][38][39][40] and used extensively in biophotonics and biomedical optics [41]. These models are based on the so-called 'scalar' approach, i.e., when the incident light is assumed to be incoherent and non-polarized. In analogy to these MC techniques mentioned above in our MC model a photon packet is assigned with an initial statistical weight and injected into a semi-infinite medium. Subsequently, the photon packet undergoes a sequence of events representing light-tissue interaction such as scattering, absorption, reflection/refraction at the medium boundary, until it is either fully absorbed or leaves the medium. At each event, the position of photon packets are updated and its statistical weight is reduced by a factor e −(µ a +µ s )L . A new scattering direction of the photon packet is chosen using Henyey-Greenstein phase function [42]: where θ is the polar scattering angle (θ ∈ [0, π]) and g is the anisotropy factor, defined as the average cos θ = 1 −1 cos θ F HG (cos θ) d cos θ for any given phase function (g ∈ [−1, 1]). The MC method of the photon migration through turbid tissue-like scattering medium is based on a simulation of a large number of photon packet trajectories (∼10 6 -10 9 ) consisting of N scattering events (N ∼ 10 2 -10 3 ). Depending on the computational performance of the available facilities this procedure could be extremely time-consuming and might take hours or even several days to complete. We demonstrated that the MC approach can be effectively sped up allowing obtaining the results nearly in real-time [33,34] using General-purpose computing on graphics processing units (GPGPU). This was achieved by utilizing the Compute Unified Device Architecture (CUDA) parallel computing platform created by NVIDIA Corporation for its hardware, allowing simultaneous execution of thousands of lightweight parallel threads [43].
The block diagram showing the memory allocation scheme for storing photon packet trajectories and the calculation of final intensities is shown in Figure 2. Each photon trajectory is represented by a linearized 1D-array of structures and is stored in the GPU's global memory. A lightweight parallel thread is responsible for simulating the propagation of one photon trajectory in the scattering medium. A photon trajectory consists of a set of points in 3D Cartesian coordinate space defined by three floating point numbers r(i) = X i , Y i , Z i corresponding to coordinates at the i-th scattering event. The length of the trajectory is limited by the maximum number of scattering events N, which is in our case 10 3 , and depending on the particular task can be extended to 10 4 -10 5 .
Certain sets of parameters of the photon packets, such as the number of scattering events N i , path length L i , polarization vector − → P i , scattering angle θ, angle of detection α i , etc. are saved for the subsequent processing of the detected spatial distribution of the intensity of simulated light. Therefore, a two-pass MC simulation of photon packet propagation in scattering medium is used. In the first pass, the MC code generates photon packet trajectories satisfying the detection conditions, e.g., size and position of the detector, numerical aperture, etc. The trajectories are stored in the GPU's global memory and used in the second pass to calculate the resulting intensity of scattered light, as presented in Figure 2, left. The generation of photon trajectories repeats until the maximum number of photon packets have been reached in a pixel of the detector (typically, 10 5 or more). The pixels are arranged in a 3D dynamic array [x × y × N ph ] to mimic a real detector (e.g., CCD camera). The x and y are the indices of the array corresponding to an individual pixel bins (calculated from the r(N s ) of the photon packet trajectory) at the detecting area, N ph is the total number of detected photon packets that defines the unique index corresponding to the trajectories preserved at the pixel for further calculation. In the current scenario, the amount of GPU memory required for storing these generated trajectories is ∼3 Gb in case of single floating point performance (FPP), and ∼6 Gb for double FPP. This two-step MC approach permits the 'recycling' of the same trajectories for calculating various quantities of light scattered in the medium, including polarization, phase shift and others. It also avoids a need to re-generate photon packet trajectories which is the most time consuming procedure of MC modeling. By integrating fast parallel computations on graphics cards with modern web technologies, such as HTML5 and JavaScript, a computational platform for simulations of light transport in turbid media has been developed [33,34]. The key idea behind this computational platform is the creation of a universal tool to imitate the results of typical experiments utilized in major applications in biomedical optics. A recently updated version of the user interface of the tool where an application that can be selected by clicking on a particular icon is shown in Figure 2, right. The platform is available to the worldwide community through the cloud-based interface at http://www.lighttransport.net/.
To describe and track polarization changes, Jones [44,45] and/or Stokes-Mueller [46,47] formalism are required. In the current study, in frame of further development of a unified MC model, the Jones-analogues-based formalism is adopted to handle propagation of each partial component of CVB traveling through a turbid tissue-like scattering medium. The concept of iterative procedure of the solution of Bethe-Salpeter equation [48] is used to describes the transfer of a pair of two complex-conjugated fields, incident into the point of source R S (x, y) with the wave vector k i and outgoing in the detecting point R D (x, y) with the wave vector k D . Thus, the propagation of polarized light through the turbid medium is presented in terms of scattering orders, as a a sequence of geometric transformations after each scattering event. Here, the trajectories of the photon packet are weighted (W) in accordance with the polarization state, and the polarization vector of the scattered wave − → P i is transformed upon the i-th scattering event, defined as [24]: where − → e i is the unit vector aligned along the trajectory element of a photon packet after the ith scattering event. Thus, tensorÛ i = [Î − − → e i ⊗ − → e i ] is presented as: iX −e iX e iY −e iX e iZ −e iX e iY 1 − e 2 iY −e iY e iZ −e iX e iZ e iX e iZ and, thus, the chain of projection operatorsÛ i transforms the initial polarization − → P 0 upon a sequence of n scattering events to the final polarization − → P n , i.e.,

− →
Consequently, propagation of co-and cross-polarized components of the electromagnetic field in the medium is described along the same trajectories obtained for the scalar field. To link the scalar and vector nature of electromagnetic fields according to the optical theorem [49,50] the Rayleigh factor is taken into account at every scattering event for electromagnetic field [24,51]: The reflection and/or refraction at the medium boundary is taken into account by splitting the photon packet into the transmitted and reflected parts [52]. Thus, the statistical weight of the detected photon packet for the Transverse Magnetic (TM) and Transverse Electric (TE) components is defined as: where R TM (α i ) and R TE (α i ) are the Fresnel coefficients, M is the number of reflections on the medium boundary.
Finally, the co-I and cross-I ⊥ polarized intensities of scattered light at the detection point (2D array indexes) are counted as a superposition of all polarization states detected at the pixel (x, y): where N is the number of scattering events experienced by i-th photon packet along its trajectory L within the medium. In our study the Stokes formalism is applied for definition of the degree of polarization DR [53]: where S 0 , S 1 , S 2 and S 3 are the Stokes vector components. As the co-I and cross-I ⊥ polarized components are, respectively, oriented along X and Y axes of the Cartesian system that associated with the scattering medium, for the horizontal/vertical linearly polarized light both S 2 and S 3 are 0. Thus, the spatial distribution of DR(x, y) at the surface of the medium can be presented as: where I tot (x, y) = I (x, y) + I ⊥ (x, y). Thus, the DR can be considered as an analogy of fringe contrast, which is typically used as a measure of the visibility quality in the interference measurements [54], known also as interference visibility and/or fringe visibility. Figure 3 shows the results of simulation of the spatial distribution of co-I , cross-I ⊥ polarized components, and DR for CVB propagated through the scattering medium of various optical density. As one can see for the highly optically dense medium, the CVBs is almost totally depolarized. The developed MC model has been verified by a direct comparison of the results of modeling of CVB and Gaussian beam propagation through the turbid scattering medium. In terms of accuracy of the modelling, we emphasize that early, it has been shown that MC approach is able to imitate the results of actual experiment [55]. The obtained results are also in a good agreement with the results of alternative studies of deep transmission of LG vortex beams through the turbid media in comparison with the Gaussian beam [56].

Results and Discussion
As one can see, the range of the degree of fringe contrast (DR) for CVB becomes at least twice higher compared to the conventional linearly polarized Gaussian beam (see Figure 4). In addition, the modeling of polarization transfer through the scattering medium was verified against the results of exact analytical solution and alternative MC models, presented in Table 1.  The exact Milne solution for the Rayleigh (g = 0), non-coherent (temporal coherence length of laser radiation l c = 0) scattering gives I /I ⊥ ≈ 1.92 [57][58][59][60], and DR = 0.31 [57][58][59][60], while the results of MC modelling suggests I /I ⊥ = 1.93 and DR = 0.3172; a similar value of DR = 0.33 has been obtained elsewhere [25].
To confirm the results of the theoretical predictions, we built an experimental setup for the generation of CVB that is based on the design of previous developments [61] and depicted in Figure 5. A collimated 690 nm laser beam, produced by Crystalaser (DL690-030 with nominal coherence length less than 10 m) propagates through a half wave plate (HW) to obtain horizontal polarization, a pinhole for mode cleaning and then passes through a 50:50 beam splitter (BS) reaching the Pluto NIR II reflecting spacial light modulator (SLM) from Holoeye (see Figure 5). The SLM imprints the phase (e i(lφ) + e −i(lφ+t) ) in the wave front, where l is the topological charge associated with the beam and t is a prism phase that separates spatially both contrary topological charges. The beam resultant from this prism phase is reflected in M3 and its polarization is rotated 90 degrees with respect to the beam that is reflected from the BS. Both beams are recombined in the polarizing beams splitter (PBS), and then go through a quarter wave plate resulting in the superposition of two contrary circular and topological polarizations, which results in the desired CVB. A telescopic system (TS) is used to expand the beam at the sample, which consist of a square cross-section quartz cuvette (4 × 1 × 1 cm 3 ) containing water with polystyrene microspheres (by Kisker Biotech) of 3 µm diameter. The beam then goes to a lens for a 2X magnification and then to a CMOS camera for visualization purposes. A combination of a quarter waveplate (QW), a HW and a PBS are used for polarization analysis. Figure 6 shows the measurements results for a CVBs transmitting through turbid tissue-like phantom scattering media. For the medium with isotropic scattering, the CVBs is almost immediately depolarized. On the contrary strong forward scattering preserves polarization within the tissue. When scattering is highly anisotropic, which is common for many biological tissues [1], CVBs can be detected and analyzed at larger distances from the light source.  Figure 7 shows the results of experimental measurements of CVB propagation through the turbid medium of various OD. Both the experimental and modelling results are presented in terms of fringe contrast for the media with different OD (OD = µ s (1 − g)d), where d is the thickness of the medium/cuvette. To variate scattering properties/OD of the medium the solutions of micro-spheres in water of different concentration were used. As one can see the modeling results are well agreed with the results of experimental measurements. The MC model has been compiled, tested and proven to work on Windows 10/Ubuntu GNU Linux 18.04 LTS platforms. The time required to simulate 10 9 photon packet trajectories is dependent on the parameters of the scattering medium, which are in a range from µ s = 30 mm −1 , g = 0.9, n = 1.4 to µ a = 0.001 mm −1 , and takes approximately 20 min for the optical density of 1.08.

Summary and Conclusions
In the current study, we have introduced the electric field MC approach specifically developed for the fast modeling of propagation of the complex structured light through a turbid tissue-like scattering medium. The model takes into account the coherent properties of light, the influence of total reflection and refraction at the medium boundary, and mutual interference of polarization components of the wave front of CVB due to the scattering in the medium. We show that after propagation through the turbid tissue-like scattering medium, the degree of fringe contrast for CVB becomes at least double in comparison to the conventional linearly polarized Gaussian beam. The results of simulations agree with the results of experimental studies. Both experimental and theoretical results suggest that there is a high potential in application of CVB in the diagnosis of biological tissues. The developed MC approach is extremely flexible to operate with different optical properties of the medium, including scattering, absorption, refractive index, and anisotropy of scattering, as well as with different laser beam structures and complexity of the medium composition. The presented MC approach is part of the online model, and can be easily extended to model propagation of complex structured light with angular momentum in random, optically active and/or birefringent scattering media, and potentially quantum entanglements therein.