1.IntroductionRecent advances of the biomedical polarimetry have clearly demonstrated that circularly polarized light (CPL) can be effectively used for overall characterization of biological tissues with optical anisotropy,1–3 including detection of the senile Alzheimer’s plaques4,5 and quantitative evaluation of the cervical intraepithelial neoplasia.6,7 Proper exploration of the CPL-tissue interaction requires accurate self-consistent descriptive simulation tools.1,8,9 Monte Carlo (MC) based approaches are widely recognized as efficient tools for analyzing light scattering by biological tissues and turbid medium.10–14 In biophotonics, MC methods, such as MCML,15 created by L. Wang and S. Jacques, were originally designed to simulate scalar light transport within turbid scattering medium16,17 and were fundamentally relying on the radiative transfer equation (RTE).18–20 As significant role of polarized light in extending diagnostic capabilities of biomedical tools became apparent,21,22 MC methods evolved accordingly resulting in many practical and popular tools particularly developed by Ramella-Roman, Prahl, and Jacques,23,24 Hielscher,25,26 Wang,27 and Xu.28 Fundamental ground for these polarized MC approaches was established by the vector radiative transfer equation (VRTE), which represents a system of equations for each Stokes parameter and can be rigorously derived from the Maxwell electromagnetic theory.29–31 At the same time, an approach based on the iterative solution to Bethe–Salpeter (BS) equation 19,32–34 utilizing Jones vector formalism has been demonstrated to be effective for polarization tracking of MC-photons within turbid tissue-like medium and simulation of coherent backscattering.13,14,35–39 Recently, it has been shown on the fundamental level that VRTE- and BS-based approaches are equivalent under certain conditions.40 Advantages of the BS-based approach involve a direct relation to the analytic Milne solution and intuitive physical interpretation of the multiple scattering process via ladder diagrams. Modern implementations of the polarization-resolved MC14,41 aim to provide a comprehensive description of polarized light scattering with either Jones or Mueller formalism, depending on the representation of the polarization state.42 Most interest is shown in CPL which, unlike linearly polarized light, possesses a unique sense of directional awareness allowing to determine if light was forward or backscattered due to its intrinsic angular momentum associated with helicity35,39,43 [see Fig. 1(a)]. This peculiar property of CPL is a manifestation of anisotropy of scattering8 and is also known as polarization memory.44–46 Stokes vector polarimetry approach with the Poincaré sphere as a graphical tool is viewed as one of the most fitting instruments for light characterization with account for helicity [see Fig. 1(b)]. In this work, we address the conservation of the polarization memory and penetration depth of the CPL scattered in turbid tissue-like medium. We introduce an MC modeling approach specially developed to unify and generalize BS-based simulation of linearly, circularly, and/or elliptically polarized light propagation. For the first time we express the BS-based MC model in terms of the Stokes–Mueller formalism and show that our approach efficiently allows to compute Jones and Stokes vectors, Mueller matrix components, and all degrees of polarization. We explore the evolution of the CPL depolarization while propagating within turbid tissue-like scattering medium and consider the dynamic binding of circular polarization memory with the helicity flips occuring along the light path length within the medium. 2.Theory2.1.Relation Between Stokes and Jones FormalismStokes vector is traditionally defined for the fully polarized light in the following form:43 Here, denotes the imaginary unit, asterisk corresponds to complex conjugation, , is a complex electric field of the plane wave propagating along axis (wave vector ), with being wave real amplitudes multiplied by complex factor with position , corresponding to phases, and , being wave complex amplitudes. Both complex fields , can be decomposed into real () and imaginary () parts: In terms of Jones formalism, it can be written as Here, is a polarization state described by the non-normalized Jones vector. We emphasize that expression Eq. (3) implies that an arbitrarily polarized electromagnetic field can be considered as a superposition of two linearly polarized fields and containing information on the phase difference between them. Jones vectors for all of the pure polarization states42,43 can be represented in this manner. In particular, for linear polarized light along axis and along axis , we have Here, . It is possible to write down both linear polarization vectors with account for non-zero phase shifts. For example, in case : Similarly, linearly polarized light components along diagonal directions can be expressed as In the following, we will mostly consider Jones vectors for the right circular polarization (RCP) and left circular polarization (LCP): By substituting field components Eqs. (2) and (3) into the definition Eq. (1) and performing some straightforward algebra, we arrive at the following expressions for the Stokes vector: It is important to note that here all variables are real-valued and that components are in fact parts of the real-valued linearly polarized e/m waves , . Established relation Eq. (5) is the fundamental one to relate Stokes formalism with the existing BS technique developed to trace evolution of Jones polarization vector along MC-photon trajectories.13,19,47 Stokes formalism enables to immediately recognize the CPL helicity flips appearing as the Stokes vector locus crossing equator on the Poincaré sphere [see Fig. 1(b)]. We note that Eqs. (1)–(5) are written in the local reference frame of the wave. 2.2.Degrees of PolarizationIn order to consider partially polarized light field averaging procedures are commonly used. This can clearly be seen on the example of the Wolf’s coherence matrix :48 Here, . With Eq. (6) we have also provided a connection between coherence matrix and Stokes parameters of the partially polarized light. Brackets correspond to the field averaging procedure. Traditionally, time-averaging with respect to the detector finite integration time is performed, along with spectral and spatial averaging defined by the resolution of the detector.42,48 In this work, brackets correspond to the averaging of MC photon intensities. This approach will be covered in Sec. 3.3. For partially polarized light, following definitions43,48 for the degrees of polarization based on the coherence matrix and Stokes approaches are used: Here, DoP is the total degree of polarization, DoLP is the degree of linear polarization, and DoCP is the degree of circular polarization, . Partially polarized light can be decomposed into fully polarized and non-polarized parts:43 Or, alternatively, partially polarized light can be treated as a superposition of two oppositely polarized waves:43 These expressions can be rewritten in more compact form by using Stokes parameters normalized to the intensity of the fully polarized component: This definition allows to compute the Stokes vector values that are typically provided, e.g., by ThorLabs polarimeters.49 In addition, we can assume that when (all Stokes components of the fully depolarized part are equal to zero). Then Eq. (10) takes the form and Eq. (11) is written asNow in both equations . Important specific cases of the expressions Eqs. (13) and (14) include decomposition of the CPL into the fully polarized right- and left-handed parts and decomposition of the linearly polarized light into orthogonal components. For the first case, we rewrite Eq. (13) as after terms regroup arriving atThis alternative form of the expression Eq. (14) allows to write down expressions for the co- and cross-polarized light components via DoCP: Here, corresponds to the RCP light and corresponds to the LCP light. DoCP value can then be estimated as We note that this expression has to be treated with care: when , we supposedly arrive at negative DoCP values. However, this does not actually contradict the definition Eq. (9), because expression Eq. (16) is derived under the assumption that RCP intensity is always larger than LCP one, as follows from Eq. (15). Otherwise, we should appropriately rewrite these equations, arriving at , which generally results in fully complying with Eq. (9). Similar decomposition can be written for the second case when light is linearly polarized: which in turn reduces to when . Here, and all polarization degrees are within limits. Intensities of horizontally and vertically polarized light can be obtained from Eq. (18) to express . This expression for has been used throughout most of the previous works.13 DoLP also involves intensities of light linearly polarized along and axes:43Here, . Now, we have established theoretical background and can proceed with the description of the developed MC approach. 3.Monte Carlo Based on the Bethe–Salpeter Equation3.1.Tracking of the Jones Polarization VectorWithin the BS-based MC model,13,19,33 a large amount () of MC-photons with pre-defined statistical weight , is launched from the source oriented under angle to the surface normal, propagates through the turbid medium and statistics is collected from those arrived on the detector oriented under angle to the surface normal (see Fig. 2). Here, the minus sign corresponds to the opposite direction of the detector to the surface normal as compared with the direction of the source. Turbid medium is defined by scattering coefficient , absorption coefficient , anisotropy parameter , and refractive index .18 In addition, tissue-like medium implies low contrast between refractive indices of the host medium and scatterers (e.g., cellular components, organelles, extracellular matrices, and other microstructures). In this work, we consider a uniform distribution of MC-photons, noting that in general our approach allows to simulate spatial and phase distributions for a wide variety of light beams, including Gaussian, Bessel, Hermite–Gaussian, and Laguerre–Gaussian beams with complex shape carrying orbital angular momentum (OAM). To account for these beam types, it is necessary both to ensure the appropriate initial distribution of the MC-photons relevant to the beam intensity and phase profiles and to set the correct initial directions of the MC-photons according to the Poynting vector trajectories that render energy transfer within the beam.50,51 With the next development, we plan to implement this technique in our model to investigate the conservation of OAM in tissue-like medium. Each MC-photon at the source is characterized by the initial statistical weight , Cartesian coordinates , propagation direction (defined both by beam structure and angle between source and surface normal, see Fig. 2) and, most importantly, by the initial polarization state. We introduce a real-valued vector that corresponds to the direction of the linearly polarized field.13,19,32–34,39 By assigning a pair of these vectors , to each MC-photon, we are able to define two separate independent linear polarization states similarly to Eq. (3). It is important to note that here both polarization vectors are written in the global Cartesian coordinate system and that they are orthogonal to the MC-photon unit propagation direction. If photon direction coincides with the axis, then sum of and can be interpreted as Jones vector: . We emphasize that from and we can always compute the Jones vector associated with the MC-photon and vice versa; by knowing the polarization state (Jones vector) of the MC-photon we can always reconstruct and values. After launch, all MC-photons undergo surface () interaction and are transmitted to the turbid medium layer with account for the Snell’s law and the appropriate Fresnel coefficients influencing MC-photon weights, directions, and polarization (see Sec. 3.2). In turbid medium () each MC-photon trajectory is modeled as a sequence of the elementary simulations containing limited amount of scattering events . This procedure has been thoroughly covered in the previous works.13,19,47 At each ’th scattering event, , the following computational steps are performed: random path length is computed (in this paper, we assume that and is a uniformly distributed random number), MC-photon is moved to the next position with weight attenuated according to the Beer–Lambert law (), and the next propagation direction is evaluated via inversion of the Henyey–Greenstein (HG) phase function52 where is the polar scattering angle in the MC-photon reference plane. Here, we have used the position vector and the unit direction for each scattering event , with as azimuthal and polar angles that correspond to the global Cartesian coordinates. HG function has been traditionally employed in the MC simulations as a substitute to the rigorous Mie phase function due to its high performance and the ability to provide realistic results complying with the experimental tissue measurements.15,53,54 It should be noted that, basically, any phase function can be used.55,56 If analytical inversion of is not possible, e.g., for the case of Mie scattering, then table lookup method is involved to ensure fast computational speed.At each step, we check if MC-photon path crosses the medium boundary and invoke surface refraction–transmission and detection procedures if this is the case (see Sec. 3.2). Evolution of each linearly polarized state can be traced along the MC-photon trajectory via the following procedure, which is obtained from the iterative solution to BS equation:13,14,19 Here, is the third-rank unit tensor and indicates the direct product. Tensor can be explicitly rewritten in the form of real operator :32 Most importantly, operator guarantees that the electromagnetic field remains transversal experiencing the ’th scattering event. It can be applied to both linear polarization vectors simultaneously as follows from Eq. (2), and it accounts for the helicity flips when considering pair of the polarization vectors that correspond to the circularly or elliptically polarized MC-photon [see Fig. 1(a)]. Eventually, the chain of projection operators transforms the initial polarization upon a sequence of scattering events to the final polarization :19 The same expression can be used to relate and as follows from Eq. (2). It is important to note that this procedure always ensures and to be orthogonal to the MC-photon direction at each scattering event. It means that if we rewrite and in terms of the MC-photon local reference frame using the appropriate transformation matrix, we will obtain Jones vectors with third component equal to zero. This peculiarity can be verified, e.g., numerically, but, most importantly, polarization tracing Eq. (21) does not inherently require reference frame tracking and allows one to avoid computation of the scattering and rotation matrices as proposed by the VRTE-based approaches,23,28 leading to the computational demand of polarization-enabled MC to be comparable to the demand of scalar MC. Tensor ensures that each individual MC-photon always remains fully polarized. Then Stokes vector values can be obtained for each MC-photon at any scattering event via Eq. (5) with values replaced by the corresponding components. We should explicitly note that the approach based on the Bethe–Salpeter equation was rigorously introduced for the case of pure Rayleigh scattering.32 In case of biotissues, we deal with scatterers with the size comparable to or a few times higher than the wavelength . Keeping in mind that within biological media fluctuation of the relative refractive index between the scatterer (e.g., cell component such as nucleus, ) and the surrounding medium (e.g., cytoplasm, ) is typically small (, ),18 we conclude that we actually deal with the so-called soft scattering particles.57,58 In this case, particle size should obey the relation , where . Then Rayleigh–Gans–Debye (RGD) approach can be applied to describe scattering by soft particles characterized by the non-isotropic scattering phase function.32,57,58 On these grounds, the proposed BS-based MC polarization tracing can be treated as the first-order approximation to RGD and applied to simulate polarized light scattering in biological media.19,32 We also note that in this paper non-birefringent and non-optically active medium is considered; while birefringence is known to be an important feature of biological tissues, it has been reported that, e.g., for skin it is almost impossible to observe the phase changes occurring due to birefringence at normal conditions.59 At the same time, account for birefringence can be added into the developed model by properly implementing account for the ordinary and extraordinary optical pathlengths of MC-photons influencing the phase shift and polarization state. We repeat the outlined computational steps for each scattering event until one of the following conditions is met: either (statistical weight becomes negligible as follows from the Beer–Lambert law) or the amount of scattering events becomes larger than . These limitations ensure proper trajectory tracing cut-off.19 We continue launching MC-photons until the certain amount (no less than ) arrives on the detector. Detection procedure consists of the two checks: MC-photon coordinates should lie within the detector area (, ) and refracted direction should meet the detector numerical aperture (NA) requirements. We would limit those directions by using , where is the unit vector collinear to the detector axis. Both here and in the subsequent sections, is considered to be an index of the detection event. 3.2.Interface InfluenceOperator allows us to trace the polarization evolution at each scattering event within the turbid medium, as shown by Eq. (21). However, it does not account for the phenomena occurring at the medium boundaries. In this case, the well-known Fresnel coefficients have to be applied to polarized light:48 , , , . Here, correspond to the transmission coefficients for P- and S-polarized (or and ) waves, and correspond to the reflection coefficients. We have also introduced angle of the incident light , angle of the transmitted light , and medium refractive indices . Fresnel coefficients can be complex-valued, for example, in case of total internal reflection due to Snell law . As a consequence, these coefficients cannot be separately applied to each linear polarization vector ; instead, the complex counterpart of Eq. (3) has to be reconstructed from the pair of vectors Eq. (2) prior to applying Fresnel coefficients. After that, the new reflected or transmitted vectors can be decomposed back into two separate linear polarization states, and polarization tracing procedure from Sec. 3.1 can be continued. We also have to keep in mind that Fresnel coefficients are derived in the wave’s plane of incidence.48 It means that at the event of the MC-photon interaction with the surface we have to rewrite both vectors in the corresponding reference frame , defined by the MC-photon direction and its projection to the interface of the surface, via applying proper transformation matrix. If is the index of the event of the MC-photon interaction with the surface, and is the index of the next scattering event, account for the Fresnel coefficients can be mathematically expressed in the following form: , , . Here, are polarization vectors transformed to the reference frame associated with the MC-photon’s plane of incidence. In terms of polarization vector components: For the transmission, it is enough to replace with their counterparts . At the same time, in the specific case of linearly polarized light where phase information is not usually relevant, the field has only one polarization vector , and it is possible to account for polarization changes at the interface via absolute values of Fresnel coefficients as outlined in the previous works.13 This procedure influences the absolute value of polarization vectors, and, correspondingly, the weight of each MC-photon. After account for the interface influence, both vectors are transformed back to the global reference frame. We would further use the notations and in order to emphasize that non-laboratory reference frame is used; in addition to the plane of incidence, this could be either source or detector reference frame, or local reference frame of the MC-photon. We also note that it is necessary to properly select transmitted or reflected MC-photons in multilayered medium. It can be done via implementing selection procedure following Wang15 at each interface between medium layers, adding proper account for the polarization state of the MC-photon. In this work, we consider homogeneous scattering medium with single layer. 3.3.Detected Light Intensity Components, Stokes Vector, and Polarization DegreesEach MC-photon that arrived on the detector is fully polarized and its polarization state is known from Eq. (21) with account for reflections/refractions by Eq. (22). Every detected MC-photon also possesses weight attenuated with respect to the Beer–Lambert law , where is index of the detection event for ’th MC-photon, and is the path length between two neighboring scattering events. If detector plane is parallel to the medium surface, then averaging of the MC-photon ensemble intensity components is performed as follows:34,39 For completeness, we also provide expressions for all intensities of the linearly polarized light: Here, is the Rayleigh factor derived from the optical theorem in Born approximation and is the square cosine of the scattering angle weighted by the single scattering cross-section. 13,19,32,33 For an arbitrary orientation of the detector (see Fig. 2), both and are supposed to be rewritten in the new Cartesian basis with axis being collinear to the detector axis. Stokes parameters are related to the light intensity components as Final expressions for the Stokes parameters withing the BS-based MC are Degrees of polarization can then be computed either via definitions Eqs. (7)–(9) or, equivalently, via expressions for intensity components Eqs. (16) and (19). Depending on the detection conditions, it might be necessary to compute any of the given parameters in the reference frame other than the global one, e.g., in the detector reference frame or in the local reference frame of each MC-photon. For this purpose, transformation matrix providing in the selected reference frame can be used. The obtained values can be directly substituted into Eqs. (23)–(30) providing appropriate intensity, Stokes, or degree of polarization values. 3.4.Computation of Mueller Matrix ComponentsWe have demonstrated that within the proposed MC approach, such parameters as Jones vector Eq. (21), Stokes vector for partially polarized light Eq. (30), Wolf coherence matrix Eq. (6), and degrees of polarization Eqs. (7)–(9), can be evaluated. We also stress that it is possible to compute Mueller matrix elements. We consider Mueller matrix in its general form: Mueller matrix elements are usually measured with the following setup configurations:60 Here, the first letter corresponds to the source polarization, and the second letter corresponds to the measured intensity (with analyzer); is the non-polarized light, corresponds to , to , to , to , to , and to . In terms of our model, Mueller matrix of the single detected photon can be expressed as Here, and are the real-valued vectors introduced in Sec. 3.1 and computed via Eq. (21) for incident linear polarizations , . Similarly, and are vectors computed for incident diagonal linear polarizations and . Circular polarization states and are accounted for as superpositions of and according to Eq. (4). means that this element can be obtained via rotation of by .60 Matrix Eq. (32) is valid when the detector plane coincides with the medium surface, as outlined in Sec. 3.3. Mueller matrix of the detected signal can then be obtained via the ensemble averaging procedure following the Eqs. (23)–(28): Here, corresponds to the Mueller matrix of the ’th photon, which was detected at the scattering event, and all Mueller matrix elements are independently multiplied by the scalar term for each photon. Now our formulation of the generalized BS-based polarization MC is complete. We emphasize that with Eqs. (32)–(33) we can compute Mueller matrix within one simulation, so it is not required to launch separate MC-photons with different polarization states. This factor, along with the remarks made in Sec. 3.1 [see Eq. (21)], contributes to the high computational performance of our approach. 4.Results and Discussion4.1.Setup ConfigurationOur theoretical model is oriented toward the most common experimental setups employed to study both forward (transmission) scattering and backscattering by biotissues with non-invasive diagnostic purposes.61 In particular, we verify the obtained simulation results against measurements performed with the backscattering setup, which has been thoroughly described in our previous works.4–6 In this setup, we employ multiwavelength 450 to 650 nm light source with diameter incident under on the tissue-like surface characterized by , and . In the following, these values are selected to closely match the properties of real tissues or tissue phantoms.62 Incident light is right circularly polarized. We collect the scattered depolarized signal in the detector with diameter oriented under with respect to surface normal and separated from the source by distance (see Fig. 2). In order to properly study the evolution of CPL, we use an infinity-corrected objective in the detection arm ensuring that polarimeter registers Stokes parameters that correspond to the MC-photon local reference frames. In this paper, incident beam is simulated as a plane wave (uniform distribution of MC-photons, direction defined solely by ) with and polarization vectors defined as , in the reference frame of the source. In the global reference frame, which is further employed in the scattering simulation, these vectors take the following form: , . In the model, we consider two source–detector configurations: with the angular incidence and collection of light (, ), and with the vertically positioned source and detector (). The value is scaled to the transport mean free path representing the average distance that light propagates before its direction of propagation is randomized.58,63,64 We collect detector statistics Eqs. (23)–(30) via evaluating polarization vectors in the local reference frame for each MC-photon, which corresponds to the experimental detection conditions. 4.2.Depolarization of the CPL Backscattered by Turbid Tissue-Like MediumWe investigate the process of CPL depolarization in terms of the Stokes vector and light intensity components both via processing surface signal registered by the detector (see Sec. 3.1) and via analyzing in-depth distribution of the detected response represented by sampling volume.16,17 Main results are summarized in Fig. 3. We begin the analysis by studying the intensity components of the scattered light. Figures 3(b) and 3(c) show an interplay of the oppositely polarized RCP (blue) and LCP (red) intensities upon increase of the source–detector separation . As one can see, for the short separation distances ( for the vertical setup and for the angular setup), the helicity of incident RCP light is flipped due to backscattering, and the flipped LCP light is inversely related to the emerging RCP component. The LCP light is formed due to odd number of the helicity flips occurred along the consecutive scattering events within the medium between points of incidence and detection, whereas the appearance of RCP is based on the even number of flips.44 The decrease of LCP with the increase of source–detector separation is compensated with the proportional increase of RCP light, clearly illustrating predictions Eq. (15). The RCP stream becomes dominating over the LCP at larger source–detector separation (). This allows us to conclude that the angular momentum of light is preserved and that multiple scattering maintains the helicity of incident CPL (RCP). At the flip point (demarcated by red and blue background colors), the intensities of two streams of light with opposite helicities are equalized () and their superposition originates linear polarization. The polarization memory is revealed as a flip of the backscattered CPL at the source–detector separation over the transport length (), tailing the helicity of incident RCP light. The resulting superposition of the scattered RCP and LCP light is registered by the detector as elliptically polarized light. It should be noted that elliptical polarization can be observed with any non-zero phase of the incident CPL if the plane of observation is not parallel or perpendicular to the original vibration direction of the field, which is accounted for in the developed model. We proceed with the analysis of light depolarization by comparing DoP, DoLP, and DoCP versus source–detector separation. Corresponding plots are presented in Figs. 3(d) and 3(e) along with the normalized Stokes vector components are depicted on the Poincaré sphere. DoCP represents the fraction of the CPL that is preserved or retained after the multiple scattering. With the increase of source–detector separation the DoCP is decreased due to reduction of low scattering orders contribution to the backscattered light. At a particular source–detector separation where flipped and preserved components of the backscattered CPL are equalized [see Figs. 2(b) and 2(c)], the reaches a minimum value. The depolarization minimum represents the point at which the components of scattered circularly light with opposite helicity, LCP and RCP, are superimposed. The depolarization minimum is coincided with the demarcation line between non-diffusive and diffusing path lengths of scattering photons characterized by . This phenomenon is well pronounced when utilizing the angular source–detector configuration [see Figs. 2(b) and 2(d)]. These results significantly contribute to our understanding of the depolarization processes within tissues and prove to be useful, e.g., for the advanced alignment of the experimental setup with a conventional polarimeter employed to measure Stokes parameters and degrees of polarization of the backscattered elliptically polarized light. All data present in Fig. 3 have been computed with open NA of the detector (). In order to both explore the aperture influence and validate the results toward experimental data, another set of simulations was performed with aperture limited to ensuring that only light photons meeting the condition (see Sec. 3.1) are collected from the sample surface. From Fig. 4, we find good agreement of the MC simulations with experimental measurements performed with the setup described in previous works.4–6 Our simulation parameters provided in the beginning of the results section are already adjusted to approximately match the experimental setup configuration. In the experiment, we have carried out polarization measurements of RCP light scattered by thick phantom with known optical properties (, , , at )62. We observe that limitation of the NA in the model led to the shift of the helicity flip location toward the source [ for in Fig. 4(a) as opposed to for open NA in Fig. 3(b)]. We also notice that, as shown in Figs. 3(b) and 3(c), vertical source–detector setup leads to the helicity flip position being shifted away from the source (), whereas angular source–detector placement causes helicity flip position to shift toward the source (). In other words, the larger and are, the closer helicity flip is to the source. Alternatively, this can be interpreted in terms of the medium refractive index , which modifies the effective incident and detection angles according to Snell refraction law. It should be also pointed out that depolarization composition of the backscattered CPL varies depending on the properties of turbid tissue-like disperse medium, such as its scattering characteristics, the size and composition of scattering particles implying different scattering phase functions, and the overall optical density1,8,25,64,65. 4.3.In-Depth Spatial Distribution of the CPL Components and Polarization MemoryBesides analysis of the surface response presented in the previous section, computer simulation provides an important insight on the in-depth light–tissue interaction. Sampling volume is a traditional parameter characterizing the detector depth sensitivity. Figures 3(a) and 3(f) show two-dimensional (2D) maps computed as difference between sampling volumes (SV) of the oppositely polarized RCP (blue) and LCP (red) light for several selected dimensionless source–detector separation distances . With these maps, we demonstrate that and light portions statistically propagate at different depths within the sample, as suggested in previous works of Sridhar and da Silva.66 This result is well pronounced in the angular source–detector configuration [see Fig. 3(a)]. An important outcome is the possibility to tune the penetration depth of both left- and right-polarized components of light via adjusting angle and position of the source–detector configuration. It can be clearly seen that prior to the helicity flip point [Fig. 3(a) for , Fig. 3(f) for ], and after the flip [Figs. 3(a) and 3(f) for ] in agreement with the results discussed in previous section. This proves the self-consistency of the proposed MC model and supports the capability of the model to operate with depolarized light through considering fully polarized orthogonal states. In this work, sampling volumes have been computed with16,17 Here, corresponds to the detected intensity of the ’th-th MC-photon defined by the expression under the sum sign, i.e., in Eqs. (23) and (24), is the amount of detected photons, is a path length of the ’th-th MC-photon within a voxel centered at , and is linear size of the voxel. Evaluation of Eq. (34) provides us with a 3D array depicting detector depth sensitivity within each voxel. 2D maps shown in Figs. 2(a) and 2(f) are computed as with defined via corresponding intensities. To the best of our knowledge, this is the first time when the discussed phenomena of right- and left-polarized light components possessing different sampling volumes is both quantitatively and qualitatively described with the MC simulations. To conclude this section, we point out that within our model it is possible to extensively study the distribution of polarized light within tissue in terms of polarization extinction ratio (PER):67 . PER characterizes the extent of polarization cross talk between flipped and preserved components of the backscattered CPL. Figure 5 shows the in-depth spatial distribution of the polarization memory, presented by analogy to the photon-measurement density function,68 in terms of gradient of PER computed similarly to the sampling volume in Eq. (34). PER refers to the relative intensities of LCP and RCP components and describes the mixing of flipped polarization with the orthogonal one as a result of multiple scattering interactions. Figure 5 shows a strong localization of LCP component in relation to the incident polarization state at the short () source–detector distances for both setup configurations. The linear polarization, emerged as a superposition of LCP and RCP components, demarcates areas of their localization. The wide aperture of the detector () and anisotropy of scattering result in a broad range of scattering angles of photons and their path length distribution, leading to an asymmetry of the in-depth spatial distribution, which is strengthened when both source and detector are not oriented along the normal to the surface of the turbid medium. 4.4.Mueller Matrix EvaluationFinally, in Fig. 6, we present an example of Mueller matrix elements computed by Eqs. (32) and (33). This data were obtained for the vertically positioned source and detector. Here, the detector registers the transmitted signal in area, . These results demonstrate that our developed approach is inherently capable of carrying out Mueller matrix computations. The ability to simulate Mueller matrix numerically is especially relevant because most of the experimental research on interaction of the polarized light with tissues employs Stokes–Mueller formalism as a standard.61,69 As outlined in Sec. 3.4, one of the main advantages of our approach is the ability to evaluate Mueller matrix without the need to launch multiple simulations for different incident polarization states. By presenting the established model in this paper, we aim to further develop our Mueller matrix MC with respect to applications in the course of the subsequent research. 5.ConclusionWe introduce an MC modeling approach that provides combined Jones and Stokes–Mueller formalism. Our model utilizes the polarization tracing framework based on the iterative solution to Bethe–Salpeter equation. The reflection and refraction of the linearly, elliptical, and/or CPL at the medium surface are generalized and properly included in the model. Self-consistency of the proposed model is ensured by the developed theoretical framework and confirmed by both numerical experiments and phantom measurements. One of the main advantages of the proposed approach is the ability to evaluate Mueller matrix elements, as well as other characteristics, such as sampling volumes or degrees of polarization, with single simulation. The results of modeling studies reveal the origins of the experimentally observed helicity flip that depends both on the configuration of the source–detector setup and turbid medium properties. First, we have shown that for the CPL backscattered from the turbid medium the flipped helicity survival is prevailed at the short source–detector separation (). A transition from LCP to RCP is revealed for longer distances () resulting in preservation of the helicity of incident light. Second, we have demonstrated that backscattered CPL within MC is appropriately decomposed into two fully polarized orthogonal components with opposite helicities, and their polarization state is fully defined. Third, we have reported on the different penetration depth of RCP and LCP light as demonstrated by the sampling volume simulations. And finally, we have addressed the in-depth binding of circular polarization memory with the helicity flips occurring within the medium. It should be pointed out that developed MC framework is suitable to imitate light beams of different shapes, such as traditional point sources, plane waves, Gaussian and Bessel beams, as well as complex laser beams carrying OAM (e.g., Laguerre–Gaussian) via appropriate definition of the initial MC-photon intensity and direction distributions. In addition, diverse source–detector configurations, coherent properties of incident light and arbitrary polarization states can be taken into account without further modifications of the code core components. In summary, the combined use of Jones and Stokes–Mueller formalisms in MC modeling offers benefits, such as comprehensive polarization modeling, flexibility in simulating different optical elements, accurate representation of complex optical systems, validation against experimental data, and enhanced understanding of polarization phenomena. These advantages make this approach valuable in a wide range of fields, including biomedical optics, remote sensing, atmospheric optics, and more. Code and Data AvailabilityData underlying the results are available from the corresponding author upon reasonable request. AcknowledgmentsThe authors acknowledge the support from ATTRACT II META-HiLight project funded by the European Union’s Horizon 2020 research and innovative programme under Grant Agreement No. 101004462, the Academy of Finland (Grant Project 325097), the Leverhulme Trust and The Royal Society (Ref. No.: APX111232 APEX Awards 2021). ReferencesM. Singh and I. Vitkin,
“Spatial helicity response metric to quantify particle size and turbidity of heterogeneous media through circular polarization imaging,”
Sci. Rep., 13 2231 https://doi.org/10.1038/s41598-023-29444-9 SRCEC3 2045-2322
(2023).
Google Scholar
M. Peyvasteh et al.,
“Two-point Stokes vector diagnostic approach for characterization of optically anisotropic biological tissues,”
J. Phys. D: Appl. Phys., 53 395401 https://doi.org/10.1088/1361-6463/ab9571 JPAPBE 0022-3727
(2020).
Google Scholar
A. Ushenko et al.,
“Stokes-correlometry analysis of biological tissues with polycrystalline structure,”
IEEE J. Sel. Top. Quantum Electron., 25 7101612 https://doi.org/10.1109/JSTQE.2018.2865443 IJSQEN 1077-260X
(2019).
Google Scholar
M. Borovkova et al.,
“Screening of Alzheimer’s disease with multiwavelength Stokes polarimetry in a mouse model,”
IEEE Trans. Med. Imaging, 41 977
–982 https://doi.org/10.1109/TMI.2021.3129700 ITMID4 0278-0062
(2022).
Google Scholar
M. Borovkova et al.,
“Evaluating β-amyloidosis progression in Alzheimer’s disease with Mueller polarimetry,”
Biomed. Opt. Express, 11 4509
–4519 https://doi.org/10.1364/BOE.396294 BOEICL 2156-7085
(2020).
Google Scholar
D. Ivanov et al.,
“Colon cancer detection via Poincaré sphere representation and 2D polarimetric mapping of ex vivo tissue samples,”
J. Biophotonics, 13 e202000082 https://doi.org/10.1002/jbio.202000082
(2020).
Google Scholar
B. Kunnen et al.,
“Application of circularly polarized light for non-invasive diagnosis of cancerous tissues and turbid tissue-like scattering media,”
J. Biophotonics, 8 317
–323 https://doi.org/10.1002/jbio.201400104
(2015).
Google Scholar
C. M. Macdonald, S. L. Jacques and I. V. Meglinski,
“Circular polarization memory in polydisperse scattering media,”
Phys. Rev. E, 91 033204 https://doi.org/10.1103/PhysRevE.91.033204
(2015).
Google Scholar
J. C. Ramella-Roman, I. Saytashev and M. Piccini,
“A review of polarization-based imaging technologies for clinical and preclinical applications,”
J. Opt., 22 123001 https://doi.org/10.1088/2040-8986/abbf8a
(2020).
Google Scholar
Q. Fang, F. Martelli and L. Lilge,
“Special section guest editorial: Introduction to the special section celebrating 30 years of open source Monte Carlo codes in biomedical optics,”
J. Biomed. Opt., 27 083001 https://doi.org/10.1117/1.JBO.27.8.083001 JBOPFO 1083-3668
(2022).
Google Scholar
N. Nishizawa and T. Kuchimaru,
“Depth estimation of tumor invasion in early gastric cancer using scattering of circularly polarized light: Monte Carlo simulation study,”
J. Biophotonics, 15 e202200062 https://doi.org/10.1002/jbio.202200062
(2022).
Google Scholar
V. Periyasamy and M. Pramanik,
“Advances in Monte Carlo simulation for light propagation in tissue,”
IEEE Rev. Biomed. Eng., 10 122
–135 https://doi.org/10.1109/RBME.2017.2739801
(2017).
Google Scholar
A. Doronin et al.,
“Backscattering of linearly polarized light from turbid tissue-like scattering medium with rough surface,”
J. Biomed. Opt., 21 071117 https://doi.org/10.1117/1.JBO.21.7.071117 JBOPFO 1083-3668
(2016).
Google Scholar
A. Doronin, C. Macdonald and I. Meglinski,
“Propagation of coherent polarized light in highly scattering turbid media,”
J. Biomed. Opt., 19 025005 https://doi.org/10.1117/1.JBO.19.2.025005 JBOPFO 1083-3668
(2014).
Google Scholar
L. Wang, S. L. Jacques and L. Zheng,
“MCML-Monte Carlo modeling of light transport in multi-layered tissues,”
Comput. Methods Prog. Biomed., 47 131
–146 https://doi.org/10.1016/0169-2607(95)01640-F CMPBEK 0169-2607
(1995).
Google Scholar
I. V. Meglinski and S. J. Matcher,
“Modelling the sampling volume for skin blood oxygenation measurements,”
Med. Biol. Eng. Comput., 39 44
–50 https://doi.org/10.1007/BF02345265 MBECDY 0140-0118
(2001).
Google Scholar
I. Meglinski et al.,
“Study of the possibility of increasing the probing depth by the method of reflection confocal microscopy upon immersion clearing of near-surface human skin layers,”
Quantum Electron., 32 875
–882 https://doi.org/10.1070/QE2002v032n10ABEH002309 QUELEZ 1063-7818
(2002).
Google Scholar
V. V. Tuchin, Tissue Optics: Light Scattering Methods and Instruments for Medical Diagnostics, 3rd ed.SPIE Press, Bellingham, Washington
(2015). Google Scholar
I. Meglinski, A. Doronin, Advanced Biophotonics: Tissue Optical Sectioning, 1
–72 CRC Press, Boca Raton
(2013). Google Scholar
A. Sassaroli and F. Martelli,
“Equivalence of four Monte Carlo methods for photon migration in turbid media,”
J. Opt. Soc. Am. A, 29 2110
–2117 https://doi.org/10.1364/JOSAA.29.002110 JOAOD6 0740-3232
(2012).
Google Scholar
T. Novikova et al.,
“Polarized light for biomedical applications,”
J. Biomed. Opt., 21 071001 https://doi.org/10.1117/1.JBO.21.7.071001 JBOPFO 1083-3668
(2016).
Google Scholar
I. Meglinski, T. Novikova and K. Dholakia,
“Polarization and orbital angular momentum of light in biomedical applications: feature issue introduction,”
Biomed. Opt. Express, 12 6255
–6258 https://doi.org/10.1364/BOE.442828 BOEICL 2156-7085
(2021).
Google Scholar
J. C. Ramella-Roman, S. A. Prahl and S. L. Jacques,
“Three Monte Carlo programs of polarized light transport into scattering media: part I,”
Opt. Express, 13 4420
–4438 https://doi.org/10.1364/OPEX.13.004420 OPEXFF 1094-4087
(2005).
Google Scholar
J. C. Ramella-Roman, S. A. Prahl and S. L. Jacques,
“Three Monte Carlo programs of polarized light transport into scattering media: part II,”
Opt. Express, 13 10392
–10405 https://doi.org/10.1364/OPEX.13.010392 OPEXFF 1094-4087
(2005).
Google Scholar
A. H. Hielscher, J. R. Mourant and I. J. Bigio,
“Influence of particle size and concentration on the diffuse backscattering of polarized light from tissue phantoms and biological cell suspensions,”
Appl. Opt., 36 125
–135 https://doi.org/10.1364/AO.36.000125 APOPAI 0003-6935
(1997).
Google Scholar
S. Bartel and A. H. Hielscher,
“Monte Carlo simulations of the diffuse backscattering Mueller matrix for highly scattering media,”
Appl. Opt., 39 1580
–1588 https://doi.org/10.1364/AO.39.001580 APOPAI 0003-6935
(2000).
Google Scholar
X. Wang and L. V. Wang,
“Propagation of polarized light in birefringent turbid media: time-resolved simulations,”
Opt. Express, 9 254
–259 https://doi.org/10.1364/OE.9.000254 OPEXFF 1094-4087
(2001).
Google Scholar
M. Xu,
“Electric field Monte Carlo simulation of polarized light propagation in turbid media,”
Opt. Express, 12 6530
–6539 https://doi.org/10.1364/OPEX.12.006530 OPEXFF 1094-4087
(2004).
Google Scholar
M. I. Mishchenko,
“Vector radiative transfer equation for arbitrarily shaped and arbitrarily oriented particles: a microphysical derivation from statistical electromagnetics,”
Appl. Opt., 41 7114
–7134 https://doi.org/10.1364/AO.41.007114 APOPAI 0003-6935
(2002).
Google Scholar
M. J. Raković et al.,
“Light backscattering polarization patterns from turbid media: theory and experiment,”
Appl. Opt., 38 3399
–3408 https://doi.org/10.1364/AO.38.003399 APOPAI 0003-6935
(1999).
Google Scholar
H. H. Tynes et al.,
“Monte Carlo and multicomponent approximation methods for vector radiative transfer by use of effective Mueller matrix calculations,”
Appl. Opt., 40 400
–412 https://doi.org/10.1364/AO.40.000400 APOPAI 0003-6935
(2001).
Google Scholar
E. Akkermans et al.,
“Theoretical study of the coherent backscattering of light by disordered media,”
J. Phys. France, 49 77
–98 https://doi.org/10.1051/jphys:0198800490107700
(1988).
Google Scholar
I. Meglinski et al.,
“Monte Carlo simulation of coherent effects in multiple scattering,”
Proc. R. Soc. A, 461 43
–53 https://doi.org/10.1098/rspa.2004.1369 PRLAAZ 1364-5021
(2005).
Google Scholar
V. L. Kuz’min and I. V. Meglinski,
“Backscattering of linearly and circularly polarized light in randomly inhomogeneous media,”
Opt. Spectrosc., 106 257
–267 https://doi.org/10.1134/S0030400X09020180 OPSUA3 0030-400X
(2009).
Google Scholar
I. Meglinski and V. L. Kuz’min,
“Coherent backscattering of circularly polarized optical radiation from a disperse random medium,”
Prog. Electromagn. Res. M, 21 1972
–1977 https://doi.org/10.2528/PIERM10102106
(2011).
Google Scholar
A. Doronin et al.,
“Two electric field Monte Carlo models of coherent backscattering of polarized light,”
J. Opt. Soc. Am. A, 31 2394
–2400 https://doi.org/10.1364/JOSAA.31.002394 JOAOD6 0740-3232
(2014).
Google Scholar
S. V. Gangnus, S. J. Matcher and I. Meglinski,
“Monte Carlo modeling of polarized light propagation in biological tissues,”
Laser Phys., 14
(6), 886
–891
(2004).
Google Scholar
V. L. Kuzmin and I. Meglinski,
“Anomalous polarization phenomena during light scattered in random media,”
J. Exp. Theor. Phys., 137 742
–753 https://doi.org/10.1134/S1063776110050031 JTPHES 1063-7761
(2010).
Google Scholar
V. L. Kuzmin and I. V. Meglinski,
“Dependence of the circular polarization of backscattered light in random media on anisotropy of scatterers,”
Opt. Spectrosc., 108 99
–106 https://doi.org/10.1134/S0030400X10010157 OPSUA3 0030-400X
(2010).
Google Scholar
A. Doicu and M. I. Mishchenko,
“An overview of methods for deriving the radiative transfer theory from the Maxwell equations. II: Approach based on the Dyson and Bethe–Salpeter equations,”
J. Quantum Spectrosc. Radiat. Transfer, 224 25
–36 https://doi.org/10.1016/j.jqsrt.2018.10.032
(2019).
Google Scholar
S. Yan et al.,
“Graphics-processing-unit-accelerated Monte Carlo simulation of polarized light in complex three-dimensional media,”
J. Biomed. Opt., 27 083015 https://doi.org/10.1117/1.JBO.27.8.083015 JBOPFO 1083-3668
(2022).
Google Scholar
H. G. Akarçay et al.,
“Monte Carlo modeling of polarized light propagation: Stokes vs Jones - Part I,”
Appl. Opt., 53 7576
–7585 https://doi.org/10.1364/AO.53.007576 APOPAI 0003-6935
(2014).
Google Scholar
D. Goldstein, Polarized Light, Marcel Dekker, New York
(2003). Google Scholar
F. C. MacKintosh et al.,
“Polarization memory of multiply scattered light,”
Phys. Rev. B, 40 9342
–9345 https://doi.org/10.1103/PhysRevB.40.9342
(1989).
Google Scholar
M. Xu and R. R. Alfano,
“Circular polarization memory of light,”
Phys. Rev. E, 72 065601 https://doi.org/10.1103/PhysRevE.72.065601
(2005).
Google Scholar
Y. L. Kim et al.,
“Circular polarization memory effect in low-coherence enhanced backscattering of light,”
Opt. Lett., 31 2744
–2746 https://doi.org/10.1364/OL.31.002744 OPLEDP 0146-9592
(2006).
Google Scholar
A. Bykov, A. Doronin, I. Meglinski, Deep Imaging in Tissue and Biomedical Materials: Using Linear and Nonlinear Optical Methods, 295
–322 Pan Stanford Publishing, Singapore
(2017). Google Scholar
M. Born and E. Wolf, Principles of Optics, Cambridge University Press(
(2019). Google Scholar
, “Compact polarimeter PAX1000 operation manual,”
https://www.thorlabs.com/_sd.cfm?fileName=MTN007790-D02.pdf&partNumber=PAX1000VIS
(2022).
Google Scholar
M. V. Berry and K. T. McDonald,
“Exact and geometrical optics energy trajectories in twisted beams,”
J. Opt. A: Pure Appl. Opt., 10
(3), 035005 https://doi.org/10.1088/1464-4258/10/3/035005
(2008).
Google Scholar
I. Lopushenko et al.,
“Propagation of shaped light carrying orbital angular momentum through turbid tissue-like scattering medium,”
in CLEO/Eur.-EQEC Conf. Abstr.,
(2023). https://doi.org/10.1109/CLEO/Europe-EQEC57999.2023.10231893 Google Scholar
L. G. Henyey and J. L. Greenstein,
“Diffuse radiation in the Galaxy,”
Astrophys. J., 93 70
–83 https://doi.org/10.1086/144246 ASJOAB 0004-637X
(1941).
Google Scholar
D. Toublanc,
“Henyey–Greenstein and Mie phase functions in Monte Carlo radiative transfer computations,”
Appl. Opt., 35 3270
–3274 https://doi.org/10.1364/AO.35.003270 APOPAI 0003-6935
(1996).
Google Scholar
Y. Fu and S. L. Jacques,
“Monte Carlo simulation study on phase function,”
Proc. SPIE, 6084 608408 https://doi.org/10.1117/12.657490 PSISDG 0277-786X
(2006).
Google Scholar
Handbook of Optical Biomedical Diagnostics, Second Edition, Volume 1: Light-Tissue Interaction, SPIE Press, Bellingham, Washington
(2016). Google Scholar
V. Kuz’min, A. Val’kov and L. Zubkov,
“Photon diffusion in random media and anisotropy of scattering in the Henyey–Greenstein and Rayleigh–Gans models,”
J. Exp. Theor. Phys., 128 396
–406 https://doi.org/10.1134/S1063776119020109 JTPHES 1063-7761
(2019).
Google Scholar
C. F. Bohren and D. R. Huffman, Absorption and Scattering of Light by Small Particles, Wiley(
(1983). Google Scholar
A. Ishumaru, Wave Propagation and Scattering in Random Media, Academic Press, New York
(1978). Google Scholar
M. A. Borovkova et al.,
“Role of scattering and birefringence in phase retardation revealed by locus of Stokes vector on Poincaré sphere,”
J. Biomed. Opt., 25 057001 https://doi.org/10.1117/1.JBO.25.5.057001 JBOPFO 1083-3668
(2020).
Google Scholar
A. H. Hielscher et al.,
“Diffuse backscattering Mueller matrices of highly scattering media,”
Opt. Express, 1 441
–453 https://doi.org/10.1364/OE.1.000441 OPEXFF 1094-4087
(1997).
Google Scholar
Polarized Light in Biomedical Imaging and Sensing: Clinical and Preclinical Applications, Springer(
(2022). Google Scholar
O. Sieryi et al.,
“Tissue-mimicking phantoms for biomedical applications,”
Proc. SPIE, 11363 1136312 https://doi.org/10.1117/12.2560174 PSISDG 0277-786X
(2020).
Google Scholar
A. D. Kim and M. Moscoso,
“Influence of the relative refractive index on the depolarization of multiply scattered waves,”
Phys. Rev. E, 64 026612 https://doi.org/10.1103/PhysRevE.64.026612
(2001).
Google Scholar
L. F. Rojas-Ochoa et al.,
“Depolarization of backscattered linearly polarized light,”
J. Opt. Soc. Am. A, 21 1799
–1804 https://doi.org/10.1364/JOSAA.21.001799 JOAOD6 0740-3232
(2004).
Google Scholar
D. Bicout et al.,
“Depolarization of multiply scattered waves by spherical diffusers: influence of the size parameter,”
Phys. Rev. E, 49 1767
–1770 https://doi.org/10.1103/PhysRevE.49.1767
(1994).
Google Scholar
S. Sridhar and A. da Silva,
“Enhanced contrast and depth resolution in polarization imaging using elliptically polarized light,”
J. Biomed. Opt., 21 071107 https://doi.org/10.1117/1.JBO.21.7.071107 JBOPFO 1083-3668
(2016).
Google Scholar
R. M. Azzam and N. M. Bashara, Ellipsometry and Polarized Light, Elsevier, New York
(1989). Google Scholar
S. R. Arridge,
“Photon-measurement density functions. Part I: Analytical forms,”
Appl. Opt., 34 7395
–7409 https://doi.org/10.1364/AO.34.007395 APOPAI 0003-6935
(1995).
Google Scholar
T. Novikova and J. C. Ramella-Roman,
“Is a complete Mueller matrix necessary in biomedical imaging?,”
Opt. Lett., 47
(21), 5549
–5552 https://doi.org/10.1364/OL.471239 OPLEDP 0146-9592
(2022).
Google Scholar
|
Polarization
Polarized light
Sensors
Turbidity
Mueller matrices
Light scattering
Scattering