Exploring the evolution of circular polarized light backscattered from turbid tissue-like disperse medium utilizing generalized Monte Carlo modeling approach with a combined use of Jones and Stokes–Mueller formalisms

Abstract. Significance Phase retardation of circularly polarized light (CPL), backscattered by biological tissue, is used extensively for quantitative evaluation of cervical intraepithelial neoplasia, presence of senile Alzheimer’s plaques, and characterization of biotissues with optical anisotropy. The Stokes polarimetry and Mueller matrix approaches demonstrate high potential in definitive non-invasive cancer diagnosis and tissue characterization. The ultimate understanding of CPL interaction with tissues is essential for advancing medical diagnostics, optical imaging, therapeutic applications, and the development of optical instruments and devices. Aim We investigate propagation of CPL within turbid tissue-like scattering medium utilizing a combination of Jones and Stokes–Mueller formalisms in a Monte Carlo (MC) modeling approach. We explore the fundamentals of CPL memory effect and depolarization formation. Approach The generalized MC computational approach developed for polarization tracking within turbid tissue-like scattering medium is based on the iterative solution of the Bethe–Salpeter equation. The approach handles helicity response of CPL scattered in turbid medium and provides explicit expressions for assessment of its polarization state. Results Evolution of CPL backscattered by tissue-like medium at different conditions of observation in terms of source–detector configuration is assessed quantitatively. The depolarization of light is presented in terms of the coherence matrix and Stokes–Mueller formalism. The obtained results reveal the origins of the helicity flip of CPL depending on the source–detector configuration and the properties of the medium and are in a good agreement with the experiment. Conclusions By integrating Jones and Stokes–Mueller formalisms, the combined MC approach allows for a more complete representation of polarization effects in complex optical systems. The developed model is suitable to imitate propagation of the light beams of different shape and profile, including Gaussian, Bessel, Hermite–Gaussian, and Laguerre–Gaussian beams, within tissue-like medium. Diverse configuration of the experimental conditions, coherent properties of light, and peculiarities of polarization can be also taken into account.


Introduction
Recent 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][2][3] including detection of the senile Alzheimer's plaques 4,5 and quantitative evaluation of the cervical intraepithelial neoplasia. 6,7Proper exploration of the CPL-tissue interaction requires accurate self-consistent descriptive simulation tools. 1,8,99][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. 286][37][38][39] Recently, it has been shown on the fundamental level that VRTE-and BS-based approaches are equivalent under certain conditions. 40Advantages 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 MC 14,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. 42Most 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 helicity 35,39,43 [see Fig. 1(a)].5][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.

Theory
2.1 Relation Between Stokes and Jones Formalism Stokes vector is traditionally defined for the fully polarized light in the following form: 43 E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 1 ; 1 1 7 ; 6 1 6 Here, j denotes the imaginary unit, asterisk corresponds to complex conjugation, E x ¼ Ẽ0x e jδ x e jωt , E y ¼ Ẽ0y e jδ y e jωt is a complex electric field of the plane wave propagating along z axis (wave vector k ↑↑ e z ), with Ẽ0x ; Ẽ0y being wave real amplitudes multiplied by complex e −jkr factor with position r, δ x ; δ y corresponding to phases, and E 0x ¼ Ẽ0x e jδ x , E 0y ¼ Ẽ0y e jδ y being wave complex amplitudes.Both complex fields E x , E y can be decomposed into real (R) and imaginary (I) parts: In terms of Jones formalism, it can be written as E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 3 ; 1 1 7 ; 4 1 1 Here, jJi 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 RðjJiÞ and IðjJiÞ containing information on the phase difference δ ¼ δ y − δ x between them.Jones vectors for all of the pure polarization states 42,43 can be represented in this manner.In particular, for linear polarized light along x axis jHi and along y axis jVi, we have It is possible to write down both linear polarization vectors with account for non-zero phase shifts.For example, in case δ x ¼ δ y ¼ π∕4: Similarly, linearly polarized light components along diagonal directions can be expressed as ; t e m p : i n t r a l i n k -; s e c 2 . 1 ; 1 1 7 ; 1 8 6 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 E components are in fact parts of the real-valued linearly polarized e/m waves RðjJiÞ, IðjJiÞ.
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 MCphoton trajectories. 13,19,47Stokes 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.

Degrees of Polarization
In 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 J: 48 E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 6 ; 1 1 4 ; 5 3 6 Here, J xx J yy − J xy J yx ≥ 0. With Eq. ( 6) we have also provided a connection between coherence matrix and Stokes parameters ðI; Q; U; VÞ of the partially polarized light.Brackets hi correspond to the field averaging procedure.Traditionally, time-averaging hFðtÞi ¼ lim T→∞ 1 2T ∫ ∞ −∞ FðtÞdt with respect to the detector finite integration time T is performed, along with spectral and spatial averaging defined by the resolution of the detector. 42,48In this work, brackets hi correspond to the averaging of MC photon intensities.This approach will be covered in Sec.3.3.For partially polarized light, following definitions 43,48 for the degrees of polarization based on the coherence matrix and Stokes approaches are used: E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 7 ; 1 1 4 ; 3 9 1

DoP
E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 8 ; 1 1 4 ; 3 3 3 E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 9 ; 1 1 4 ; 2 8 6 Here, DoP is the total degree of polarization, DoLP is the degree of linear polarization, and DoCP is the degree of circular polarization, DoP 2 ¼ DoLP 2 þDoCP 2 .Partially polarized light can be decomposed into fully polarized and non-polarized parts: 43 E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 1 0 ; 1 1 4 ; 2 1 7 Or, alternatively, partially polarized light can be treated as a superposition of two oppositely polarized waves: 43 E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 1 1 ; 1 1 4 ; 1 2 6 These expressions can be rewritten in more compact form by using Stokes parameters normalized to the intensity of the fully polarized component: E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 1 2 ; 1 1 7 ; 7 1 2 This definition allows to compute the Stokes vector values that are typically provided, e.g., by ThorLabs polarimeters. 49In addition, we can assume that Q n ¼ U n ¼ V n ¼ 0 when DoP ¼ 0 (all Stokes components of the fully depolarized part are equal to zero).Then Eq. ( 10) takes the form E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 1 3 ; 1 1 7 ; 6 3 0 and Eq. ( 11) is written as ; t e m p : i n t r a l i n k -; e 0 1 4 ; 1 1 7 ; 5 5 2 Now in both equations 0 ≤ DoP ≤ 1.
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 E Q -T A R G E T ; t e m p : i n t r a l i n k -; s e c 2 . 2 ; 1 1 7 ; 4 3 6 after terms regroup arriving at E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 1 5 ; 1 1 7 ; 3 5 5 This alternative form of the expression Eq. ( 14) allows to write down expressions for the co-and cross-polarized light components via DoCP: E Q -T A R G E T ; t e m p : i n t r a l i n k -; s e c 2 . 2 ; 1 1 7 ; 2 6 1 Here, I R corresponds to the RCP light and I L corresponds to the LCP light.DoCP value can then be estimated as E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 1 6 ; 1 1 7 ; 2 0 2 We note that this expression has to be treated with care: when I L > I R , 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 DoCP ¼ ðI fully complying with Eq. ( 9).Similar decomposition can be written for the second case when light is linearly polarized: which in turn reduces to E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 1 8 ; 1 1 4 ; 6 4 3 when U ¼ 0. Here, DR ¼ jQj∕I and all polarization degrees are within ½0;1 limits.Intensities of horizontally I k and vertically I ⊥ polarized light can be obtained from Eq. ( 18) to express This expression for DR has been used throughout most of the previous works. 13DoLP also involves intensities of light linearly polarized along þ45 deg and −45 deg axes: 43 Here, . 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 Equation 3.1 Tracking of the Jones Polarization Vector Within the BS-based MC model, 13,19,33 a large amount (N inc > 10 9 ) of MC-photons with predefined statistical weight W j , j ¼ ½1:::N inc is launched from the source oriented under θ i angle to the surface normal, propagates through the turbid medium and statistics is collected from those N ph < N inc arrived on the detector oriented under −θ d 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 μ s , absorption coefficient μ a , anisotropy parameter g, and refractive index n. 18n 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,51With 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 W 0 j , Cartesian coordinates ðx 0 ; y 0 ; 0Þ, propagation direction s 0 (defined both by beam structure and angle θ i between source and surface normal, see Fig. 2) and, most importantly, by the initial polarization state.34]39 By assigning a pair of these vectors P x ¼ ðP xx ; P xy ; P xz Þ, P y ¼ ðP yx ; P yy ; P yz Þ 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 ðx; y; zÞ and that they are orthogonal to the MC-photon unit propagation direction.If photon direction coincides with the z axis, then sum of P x ∼ RðjJiÞ and P y ∼ IðjJiÞ can be interpreted as Jones vector: jJi ¼ P x þ jP y .We emphasize that from P x and P y 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 P x and P y values.
After launch, all MC-photons undergo surface (z ¼ 0) 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 (z > 0) each MC-photon trajectory is modeled as a sequence of the elementary simulations containing limited amount of scattering events N scatt .This procedure has been thoroughly covered in the previous works. 13,19,47At each i'th scattering event, i ¼ ½1:::N scatt , the following computational steps are performed: random path length l i ¼ − ln ξ∕μ s is computed (in this paper, we assume that μ a ≪ μ s and ξ ∈ ð0;1 is a uniformly distributed random number), MC-photon is moved to the next position r i ¼ r i−1 þ s i l i with weight attenuated according to the Beer-Lambert law (W i ¼ W i−1 e −μ a l i ), and the next propagation direction s iþ1 is evaluated via inversion of the Henyey-Greenstein (HG) phase function 52 E Q -T A R G E T ; t e m p : i n t r a l i n k -; s e c 3 .1 ; 1 1 7 ; 2 0 9 where θ 0 is the polar scattering angle in the MC-photon reference plane.Here, we have used the position vector r i ¼ ðx i ; y i ; z i Þ and the unit direction for each scattering event s i ¼ ½s X ; s Y ; s Z i ¼ ½sin θ cos φ; sin θ sin φ; cos θ i , 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,54It should be noted that, basically, any phase function p can be used. 55,56f analytical inversion of p is not possible, e.g., for the case of Mie scattering, then table lookup method is involved to ensure fast computational speed.

RCP LCP
5][6] Sample with known optical properties is illuminated with RCP light.Possible MCphoton trajectories with zero, one, and two backscattering events and with photon-surface interactions are presented.Each backscattering event causes a helicity flip represented by the color of the direction arrow.The experimental configuration involves supercontinuum fiber laser source filtered by the acousto-optic tunable filter.The resulting RCP is produced with the half-wave and quarter-wave plates and is focused on the medium surface under θ i angle.The detector is oriented under −θ d angle to the surface normal, collects backscattered light with 20× objective lens and measures Stokes parameters of the registered light with a polarimeter. 49The inset shows simulated sampling volumes for RCP and LCP light components at the relatively large source-detector separation distance ρ (see Sec. 4.3 for more details).
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 P x ; P y can be traced along the MC-photon trajectory r i ; i ¼ ½1:::N scatt via the following procedure, which is obtained from the iterative solution to BS equation: 13,14,19 E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 2 0 ; 1 1 4 ; 6 7 6 Here, Î is the third-rank unit tensor and ⊗ indicates the direct product.Tensor ½ Î − s i ⊗ s i can be explicitly rewritten in the form of 3 × 3 real operator Ûi : 32 Most importantly, operator Ûi guarantees that the electromagnetic field remains transversal experiencing the i'th scattering event.It can be applied to both linear polarization vectors P x ; P y 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 ÛN ÛN−1 ÛN−2 : : : Û2 Û1 of projection operators transforms the initial polarization P x 0 upon a sequence of N scattering events to the final polarization P x N : 19 E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 2 1 ; 1 1 4 ; 4 5 5 The same expression can be used to relate P y N and P y 0 as follows from Eq. ( 2).It is important to note that this procedure always ensures P x i and P y i to be orthogonal to the MC-photon direction s i at each scattering event.It means that if we rewrite P x i and P y i 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 Ûi 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 E values replaced by the corresponding P x i ; P y i components.
We should explicitly note that the approach based on the Bethe-Salpeter equation was rigorously introduced for the case of pure Rayleigh scattering. 32In 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 n r between the scatterer (e.g., cell component such as nucleus, n s ) and the surrounding medium (e.g., cytoplasm, n m ) is typically small (jn r − 1j < 0.1, n r ¼ n s ∕n m ), 18 we conclude that we actually deal with the so-called soft scattering particles. 57,58In this case, particle size d should obey the relation kdjn r − 1j ≪ 1, where k ¼ 2π∕λ.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,58On 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,32We 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. 59At 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 W i < 10 −4 (statistical weight becomes negligible as follows from the Beer-Lambert law) or the amount of scattering events N scatt becomes larger than 10 3 .These limitations ensure proper trajectory tracing cut-off. 19We continue launching MC-photons until the certain amount (no less than N ph ¼ 10 5 ) arrives on the detector.Detection procedure consists of the two checks: MC-photon coordinates should lie within the detector area ) and refracted direction s N should meet the detector numerical aperture (NA) requirements.We would limit those directions by using acosðs N • s d Þ < NA, where s d ¼ ½sinð−θ d Þ; 0; cosð−θ d Þ is the unit vector collinear to the detector axis.Both here and in the subsequent sections, N is considered to be an index of the detection event.

Interface Influence
Operator Ûi 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, T P ; T S correspond to the transmission coefficients for P-and S-polarized (or jHi and jVi) waves, and R P ; R S correspond to the reflection coefficients.We have also introduced angle of the incident light θ c , angle of the transmitted light θ t , and medium refractive indices n 1;2 .Fresnel coefficients can be complex-valued, for example, in case of total internal reflection due to Snell law n 1 sin θ c ¼ n 2 sin θ t .As a consequence, these coefficients cannot be separately applied to each linear polarization vector P x;y ; 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. 48It means that at the event of the MC-photon interaction with the surface we have to rewrite both P vectors in the corresponding reference frame ðx 0 ; y 0 ; z 0 Þ, defined by the MC-photon direction and its projection to the interface of the surface, via applying proper transformation matrix.
If i − 1 is the index of the event of the MC-photon interaction with the surface, and i is the index of the next scattering event, account for the Fresnel coefficients can be mathematically expressed in the following form: Here, P 0 are polarization vectors transformed to the reference frame associated with the MC-photon's plane of incidence.In terms of polarization vector components: E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 2 2 ; 1 1 7 ; 2 8 9 For the transmission, it is enough to replace R P ; R S with their counterparts T P ; T S .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 P x , and it is possible to account for polarization changes at the interface via absolute values jT P j 2 ; jT S j 2 ; jR P j 2 ; jR S j 2 of Fresnel coefficients as outlined in the previous works. 13This procedure influences the absolute value of polarization vectors, and, correspondingly, the weight of each MC-photon.After account for the interface influence, both P 0 vectors are transformed back to the global ðx; y; zÞ reference frame.We would further use the notations ðx 0 ; y 0 ; z 0 Þ and P 0 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 Wang 15 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.

Detected Light Intensity Components, Stokes Vector, and Polarization
Degrees Each 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 W N j ¼ W 0 j expð−μ a P N j i¼1 l i Þ, where 0 < N j < N scatt is index of the detection event for j'th MC-photon, and l i 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 E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 2 3 ; 1 1 4 ; 6 2 4 E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 2 4 ; 1 1 4 ; 5 6 7 For completeness, we also provide expressions for all intensities of the linearly polarized light: E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 2 5 ; 1 1 4 ; 5 1 4 E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 2 6 ; 1 1 4 ; 4 5 7 E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 2 7 ; 1 1 4 ; 4 1 8 E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 2 8 ; 1 1 4 ; 3 7 8 Here, Γ R ¼ 2 is the Rayleigh factor derived from the optical theorem in Born approximation and cos 2 θ is the square cosine of the scattering angle weighted by the single scattering cross-section. 13,19,32,33For an arbitrary orientation of the detector (see Fig. 2), both P x and P y are supposed to be rewritten in the new Cartesian basis with z 0 axis being collinear to the detector axis.Stokes parameters are related to the light intensity components as E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 2 9 ; 1 1 4 ; 2 6 9 Final expressions for the Stokes parameters withing the BS-based MC are E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 3 0 a ; 1 1 4 ; 2 3 2 E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 3 0 b ; 1 1 4 ; 1 7 5 E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 3 0 c ; 1 1 4 ; 1 3 6 E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 3 0 d ; 1 1 4 ; 9 6 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 P 0 in the selected reference frame ðx 0 ; y 0 ; z 0 Þ can be used.The obtained P 0 values can be directly substituted into Eqs.( 23)-( 30) providing appropriate intensity, Stokes, or degree of polarization values.

Computation of Mueller Matrix Components
We 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); O is the non-polarized light, H corresponds to I k , V to I ⊥ , P to I þ45 deg , M to I −45 deg , R to I R , and L to I L .In terms of our model, Mueller matrix M of the single detected photon can be expressed as E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 3 2 ; 1 1 7 ; 3  2 ; 0Þ.Circular polarization states jRCPi and jLCPi are accounted for as superpositions of jHi and jVi according to Eq. ( 4).M 31 ¼ M rot 12 means that this element can be obtained via rotation of M 12 by −π∕4. 60Matrix 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): E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 3 3 ; 1 1 7 ; 1 4 3 Here, M N j corresponds to the Mueller matrix of the j'th photon, which was detected at the N j scattering event, and all Mueller matrix elements are independently multiplied by the scalar term W N j Γ N j R 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.

Setup Configuration
Our 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. 615][6] In this setup, we employ multiwavelength 450 to 650 nm light source with 15 μm diameter incident under θ i on the tissue-like surface characterized by μ s ; μ a ; g, and n.In the following, these values are selected to closely match the properties of real tissues or tissue phantoms. 62Incident light is right circularly polarized.We collect the scattered depolarized signal in the detector with 50 μm diameter oriented under θ d 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 jRCPi beam is simulated as a plane wave (uniform distribution of MC-photons, direction defined solely by θ i ) with λ ¼ 640 nm and polarization vectors defined as P 0 x 0 ¼ ð1;0; 0Þ, P 0 y 0 ¼ ð0;1; 0Þ 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: P x 0 ¼ ðcos θ i ; 0; sin θ i Þ, P y 0 ¼ ð0;1; 0Þ.In the model, we consider two source-detector configurations: with the angular incidence and collection of light (θ i ¼ 55 deg, θ d ¼ 30 deg), and with the vertically positioned source and detector (θ i ¼ θ d ¼ 0).The ρ value is scaled to the transport mean free path l Ã ¼ μ −1 s ð1 − gÞ −1 representing the average distance that light propagates before its direction of propagation is randomized. 58,63,64We 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.

Depolarization of the CPL Backscattered by Turbid Tissue-Like Medium
We 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,17ain 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 ρ∕l Ã .As one can see, for the short separation distances (ρ∕l Ã < 1 for the vertical setup and ρ∕l Ã < 0.8 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. 44The 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 (ρ > l Ã ).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 (I R ¼ I L ) 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 (ρ > l Ã ), 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 Q n ; U n ; V n 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 I L and preserved I R components of the backscattered CPL are equalized [see Figs.2(b) and 2(c)], the DoCP 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 l Ã .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 (NA > 70 deg).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 NA ¼ 30 deg ensuring that only light photons meeting the condition acosðs N • s d Þ < NA (see Sec. 3.1) are collected from the sample surface.5][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 (μ s ¼ 4 mm −1 , μ a ¼ 0.05 mm −1 , g ¼ 0.8, n ¼ 1.46 at λ ¼ 640 nm) 62 .We observe that limitation of the NA in the model led to the shift of the helicity flip location toward the source [ρ∕l Ã ∼ 0.6 for NA ¼ 30 deg in Fig. 4(a) as opposed to ρ∕l Ã ∼ 0.8 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 (ρ∕l Ã ∼ 1), whereas angular source-detector placement causes helicity flip position to shift toward the source (ρ∕l Ã ∼ 0.8).In other words, the larger θ i and θ d are, the closer helicity flip is to the source.Alternatively, this can be interpreted in terms of the medium refractive index n, which modifies the effective incident and detection angles θ i ; θ d 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 density 1,8,25,64,65 .

In-Depth Spatial Distribution of the CPL Components and Polarization
Memory Besides 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 sourcedetector separation distances ρ∕l Ã .With these maps, we demonstrate that I R and I L light portions statistically propagate at different depths within the sample, as suggested in previous works of Sridhar and da Silva. 66This 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 sourcedetector configuration.It can be clearly seen that prior to the helicity flip point I L > I R [Fig.3(a) for ρ∕l Ã ¼ 0.4, Fig. 3(f) for ρ∕l Ã ¼ 0.4; 0.8], and after the flip I L < I R [Figs.3(a) and 3(f) for ρ∕l Ã ¼ 1.2] in agreement with the results discussed in previous section.This proves the selfconsistency 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 with 16,17 E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 3 4 ; 1 1 4 ; 1 0 1 ][6] Here, I N j corresponds to the detected intensity of the j'th-th MC-photon defined by the expression under the sum sign, i.e., in Eqs. ( 23) and ( 24), N ph is the amount of detected photons, L j ðr 0 Þ is a path length of the j'th-th MC-photon within a voxel centered at r 0 , and L 0 is linear size of the voxel.Evaluation of Eq. ( 34) provides us with a 3D array SVðx; y; zÞ depicting detector depth sensitivity within each voxel.2D maps shown in Figs.2(a) and 2(f) are computed as SV R ðx; 0; zÞ − SV L ðx; 0; zÞ with SV R ; SV L defined via corresponding I R ; I L intensities.To the best of our knowledge, this is the first time when the discussed phenomena of right-and leftpolarized 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 P ¼ I L ∕I R .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 (ρ < l Ã ) 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 (NA > 70 deg) and anisotropy of scattering g result in a broad range of scattering angles of photons and their path length distribution, leading to an asymmetry of the indepth spatial distribution, which is strengthened when both source and detector are not oriented along the normal to the surface of the turbid medium.

Mueller Matrix Evaluation
Finally, 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 1 × 1 cm area, ρ ¼ 0. 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,69As 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., where l Ã is the transport length.Scale step on the colorbar for regions with preserved helicity (blue) is chosen differently from the scale for regions with flipped helicity (red) in order to make the distribution details visible.

Conclusion
We 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 (ρ < l Ã ).A transition from LCP to RCP is revealed for longer distances (ρ > l Ã ) 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-Fig.6 Mueller matrix elements obtained by MC modeling for turbid scattering medium with the following optical properties: μ s ¼ 1 mm −1 , μ a ¼ 0.01 mm −1 , g ¼ 0.74, n 1.33.Here, the detector registers the signal transmitted through medium with 4 mm thickness.The dimension of each image is 1 × 1 cm, which is equal to the detector size.The individual images are represented by a two-letter combination that denotes the input polarization and the output analyzer orientation as defined in Eq. (31).
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.

Disclosures
No conflicts of interest, financial or otherwise, are declared by the authors.

Fig. 1
Fig. 1 (a) Physics of the helicity flip: when RCP light is scattered in forward direction its helicity is preserved, whereas for backscattered light its polarization state is changed to LCP.(b) Degenerate polarization states jHi; jV i; jL þ45 deg i; jL −45 deg i; jRCPi; jLCPi (defined in Sec.2.1) and helicity flip (polarization state crossing the equator) depicted on the Poincaré sphere.
Lopushenko et al.: Exploring the evolution of circular polarized light backscattered. . .

m 11 m 12 m 13 m 14 m 21 m 22 m 23 m 24 m 31 m 32 m 33 m 34 m 41 m 42 m 43 m 44
e m p : i n t r a l i n k -; s e c 3 .4

Fig. 3
Fig. 3 (a) Difference between sampling volumes for the intensity of cross-polarized I L (red) and co-polarized I R (blue) light arriving on the detector for the θ i ¼ 55 deg, θ d ¼ 30 deg setup configuration with the variable source-detector separation distance ρ expressed in terms of transport length l Ã , (b) I L ; I R as functions of the source-detector separation for the θ i ¼ 55 deg, θ d ¼ 30 deg setup, (c) the same for the θ i ¼ θ d ¼ 0 deg setup, (d) degrees of polarization DoP (red), DoCP (blue), DoLP (green), and corresponding normalized Stokes vector components Q n ; U n ; V n on the Poincaré sphere for the θ i ¼ 55 deg, θ d ¼ 30 deg setup, (e) the same for the θ i ¼ θ d ¼ 0 deg setup, (f) difference between I L ; I R sampling volumes for the θ i ¼ θ d ¼ 0 deg setup and the same source-detector separation distances ρ∕l Ã as on (a).In these simulations detector with open NA has been considered.Points on the Poincaré spheres are colored gradually from red to yellow, which corresponds to the increase of ρ∕l Ã distance.

Fig. 5
Fig. 5 Polarization memory P ¼ I L ∕I R for (a) the angular setup with θ i ¼ 55 deg, θ d ¼ 30 deg and for (b) the vertical setup with θ i ¼ θ d ¼ 0 deg as a function of the dimensionless penetration depth z∕l Ã and source-detector separation ρ∕l Ã, where l Ã is the transport length.Scale step on the colorbar for regions with preserved helicity (blue) is chosen differently from the scale for regions with flipped helicity (red) in order to make the distribution details visible.