2.1.

Relation Between Stokes and Jones Formalism

Stokes vector is traditionally defined for the fully polarized light in the following form:⁴³

Here, $j$ denotes the imaginary unit, asterisk corresponds to complex conjugation, $E_x={\tilde{E}}_{0x}{e}^{j{\delta }_x}{e}^{j\omega t}$, $E_y={\tilde{E}}_{0y}{e}^{j{\delta }_y}{e}^{j\omega t}$ is a complex electric field of the plane wave propagating along $z$ axis (wave vector $\mathbf{k}\uparrow \uparrow {\mathbf{e}}_z$), with ${\tilde{E}}_{0x},{\tilde{E}}_{0y}$ being wave real amplitudes multiplied by complex $e^{-j\mathbf{kr}}$ factor with position $\mathbf{r}$, ${\delta }_x,{\delta }_y$ corresponding to phases, and $E_{0x}={\tilde{E}}_{0x}{e}^{j{\delta }_x}$, $E_{0y}={\tilde{E}}_{0y}{e}^{j{\delta }_y}$ being wave complex amplitudes. Both complex fields $E_x$, $E_y$ can be decomposed into real ($\mathfrak{R}$) and imaginary ($\mathfrak{I}$) parts:

In terms of Jones formalism, it can be written as

Here, $|J⟩$ 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 $\mathfrak{R}(|J⟩)$ and $\mathfrak{I}(|J⟩)$ containing information on the phase difference $\delta ={\delta }_y-{\delta }_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 $|H⟩$ and along $y$ axis $|V⟩$, we have

Here, ${\delta }_x={\delta }_y=0$. It is possible to write down both linear polarization vectors with account for non-zero phase shifts. For example, in case ${\delta }_x={\delta }_y=\pi /4$:

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 $\mathbf{E}$ components are in fact parts of the real-valued linearly polarized e/m waves $\mathfrak{R}(|J⟩)$, $\mathfrak{I}(|J⟩)$.

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 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 $\mathbf{J}$:⁴⁸

Here, $J_{xx}{J}_{yy}-J_{xy}{J}_{yx}\ge 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 $⟨⟩$ correspond to the field averaging procedure. Traditionally, time-averaging $⟨F(t)⟩=\underset{T\to \infty }{\lim }\frac{1}{2T}{\int }_{-\infty }^{\infty }F(t)\mathrm{d}t$ 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,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 definitions^43,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, ${\mathrm{DoP}}^2={\mathrm{DoLP}}^2{+\mathrm{DoCP}}^2$. Partially polarized light can be decomposed into fully polarized and non-polarized parts:⁴³

Or, alternatively, partially polarized light can be treated as a superposition of two oppositely polarized waves:⁴³

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.⁴⁹ In addition, we can assume that $Q_n=U_n=V_n=0$ when $\mathrm{DoP}=0$ (all Stokes components of the fully depolarized part are equal to zero). Then Eq. (10) takes the form

Now in both equations $0\le \mathrm{DoP}\le 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

This alternative form of the expression Eq. (14) allows to write down expressions for the co- and cross-polarized light components via DoCP:

Here, $I_R$ corresponds to the RCP light and $I_L$ corresponds to the LCP light. DoCP value can then be estimated as

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 $\mathrm{DoCP}=(I_L-I_R)/(I_L+I_R)$, which generally results in $\mathrm{DoCP}=|I_R-I_L|/(I_R+I_L)$ fully complying with Eq. (9).

Similar decomposition can be written for the second case when light is linearly polarized:

Here, $I=I_{\parallel }+I_{\perp }=I_{+45\mathrm{deg}}+I_{-45\mathrm{deg}}=I_R+I_L$. Now, we have established theoretical background and can proceed with the description of the developed MC approach.

3.1.

Tracking of the Jones Polarization Vector

Within the BS-based MC model,^13,19,33 a large amount ($N_{\mathrm{inc}}>{10}^9$) of MC-photons with pre-defined statistical weight $W_j$, $j=[1...N_{\mathrm{inc}}]$ is launched from the source oriented under ${\theta }_i$ angle to the surface normal, propagates through the turbid medium and statistics is collected from those $N_{\mathrm{ph}}<N_{\mathrm{inc}}$ arrived on the detector oriented under $-{\theta }_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 ${\mu }_s$, absorption coefficient ${\mu }_a$, anisotropy parameter $g$, and refractive index $n$.¹⁸ 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).

Fig. 2

Illustration of the backscattering model with schematically depicted elements of the experimental setup.^4–6 Sample with known optical properties is illuminated with RCP light. Possible MC-photon 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 ${\theta }_i$ angle. The detector is oriented under $-{\theta }_d$ angle to the surface normal, collects backscattered light with $20\times$ objective lens and measures Stokes parameters of the registered light with a polarimeter.⁴⁹ The inset shows simulated sampling volumes for RCP and LCP light components at the relatively large source–detector separation distance $\rho$ (see Sec. 4.3 for more details).

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 $W_{0_j}$, Cartesian coordinates $(x_0,y_0,0)$, propagation direction ${\mathbf{s}}_0$ (defined both by beam structure and angle ${\theta }_i$ between source and surface normal, see Fig. 2) and, most importantly, by the initial polarization state. We introduce a real-valued vector $\mathbf{P}$ that corresponds to the direction of the linearly polarized $\mathbf{E}$ field.^{13,19,32–34,39} By assigning a pair of these vectors ${\mathbf{P}}_x=(P_{xx},P_{xy},P_{xz})$, ${\mathbf{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 ${\mathbf{P}}_x\sim \mathfrak{R}(|J⟩)$ and ${\mathbf{P}}_y\sim \mathfrak{I}(|J⟩)$ can be interpreted as Jones vector: $|J⟩={\mathbf{P}}_x+j{\mathbf{P}}_y$. We emphasize that from ${\mathbf{P}}_x$ and ${\mathbf{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 ${\mathbf{P}}_x$ and ${\mathbf{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_{\text{scatt}}$. This procedure has been thoroughly covered in the previous works.^13,19,47 At each $i$’th scattering event, $i=[1...N_{\text{scatt}}]$, the following computational steps are performed: random path length $l_i=-\ln \xi /{\mu }_s$ is computed (in this paper, we assume that ${\mu }_a\ll {\mu }_s$ and $\xi \in (\mathrm{0,1}]$ is a uniformly distributed random number), MC-photon is moved to the next position ${\mathbf{r}}_i={\mathbf{r}}_{i-1}+{\mathbf{s}}_i{l}_i$ with weight attenuated according to the Beer–Lambert law ($W_i=W_{i-1}{e}^{-{\mu }_a{l}_i}$), and the next propagation direction ${\mathbf{s}}_{i+1}$ is evaluated via inversion of the Henyey–Greenstein (HG) phase function⁵²

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 ${\mathbf{P}}_x,{\mathbf{P}}_y$ can be traced along the MC-photon trajectory ${\mathbf{r}}_i,i=[1...N_{\text{scatt}}]$ via the following procedure, which is obtained from the iterative solution to BS equation:^13,14,19

Here, $\widehat{I}$ is the third-rank unit tensor and $\otimes$ indicates the direct product. Tensor $[\widehat{I}-{\mathbf{s}}_i\otimes {\mathbf{s}}_i]$ can be explicitly rewritten in the form of $3\times 3$ real operator ${\widehat{\mathbf{U}}}_i$:³²

Most importantly, operator ${\widehat{\mathbf{U}}}_i$ guarantees that the electromagnetic field remains transversal experiencing the $i$’th scattering event. It can be applied to both linear polarization vectors ${\mathbf{P}}_x,{\mathbf{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 ${\widehat{\mathbf{U}}}_N{\widehat{\mathbf{U}}}_{N-1}{\widehat{\mathbf{U}}}_{N-2}\dots {\widehat{\mathbf{U}}}_2{\widehat{\mathbf{U}}}_1$ of projection operators transforms the initial polarization ${\mathbf{P}}_{x_0}$ upon a sequence of $N$ scattering events to the final polarization ${\mathbf{P}}_{x_N}$:¹⁹

The same expression can be used to relate ${\mathbf{P}}_{y_N}$ and ${\mathbf{P}}_{y_0}$ as follows from Eq. (2). It is important to note that this procedure always ensures ${\mathbf{P}}_{x_i}$ and ${\mathbf{P}}_{y_i}$ to be orthogonal to the MC-photon direction ${\mathbf{s}}_i$ at each scattering event. It means that if we rewrite ${\mathbf{P}}_{x_i}$ and ${\mathbf{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 ${\widehat{\mathbf{U}}}_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 $\mathbf{E}$ values replaced by the corresponding ${\mathbf{P}}_{x_i},{\mathbf{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.³² In case of biotissues, we deal with scatterers with the size comparable to or a few times higher than the wavelength $\lambda$. 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 ($|n_r-1|<0.1$, $n_r=n_s/n_m$),¹⁸ we conclude that we actually deal with the so-called soft scattering particles.^57,58 In this case, particle size $d$ should obey the relation $kd|n_r-1|\ll 1$, where $k=2\pi /\lambda$. 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.⁵⁹ 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 $W_i<{10}^{-4}$ (statistical weight becomes negligible as follows from the Beer–Lambert law) or the amount of scattering events $N_{\text{scatt}}$ becomes larger than ${10}^3$. These limitations ensure proper trajectory tracing cut-off.¹⁹ We continue launching MC-photons until the certain amount (no less than $N_{\mathrm{ph}}={105}^{}$) arrives on the detector. Detection procedure consists of the two checks: MC-photon coordinates should lie within the detector area ($-r_d+\rho \le x_N\le r_d+\rho$, $-r_d\le y_N\le r_d,z_N=0$) and refracted direction ${\mathbf{s}}_N$ should meet the detector numerical aperture (NA) requirements. We would limit those directions by using $\mathrm{acos}({\mathbf{s}}_N·{\mathbf{s}}_d)<\mathrm{NA}$, where ${\mathbf{s}}_d=[\sin (-{\theta }_d),0,\cos (-{\theta }_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.

3.2.

Interface Influence

Operator ${\widehat{\mathbf{U}}}_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:⁴⁸ $T_P=\frac{2n_1\cos {\theta }_c}{n_2\cos {\theta }_c+n_1\cos {\theta }_t}$, $T_S=\frac{2n_1\cos {\theta }_c}{n_1\cos {\theta }_c+n_2\cos {\theta }_t}$, $R_P=\frac{n_2\cos {\theta }_c-n_1\cos {\theta }_c}{n_2\cos {\theta }_c+n_1\cos {\theta }_t}$, $R_S=\frac{n_1\cos {\theta }_c-n_2\cos {\theta }_t}{n_1\cos {\theta }_c+n_2\cos {\theta }_t}$. Here, $T_P,T_S$ correspond to the transmission coefficients for P- and S-polarized (or $|H⟩$ and $|V⟩$) waves, and $R_P,R_S$ correspond to the reflection coefficients. We have also introduced angle of the incident light ${\theta }_c$, angle of the transmitted light ${\theta }_t$, and medium refractive indices $n_{\mathrm{1,2}}$. Fresnel coefficients can be complex-valued, for example, in case of total internal reflection due to Snell law $n_1\sin {\theta }_c=n_2\sin {\theta }_t$. As a consequence, these coefficients cannot be separately applied to each linear polarization vector ${\mathbf{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.⁴⁸ It means that at the event of the MC-photon interaction with the surface we have to rewrite both $\mathbf{P}$ vectors in the corresponding reference frame $(x^{\prime },y^{\prime },z^{\prime })$, 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: ${(P_x^{\prime })}_i={(P_x^{\prime })}_{i-1}·R_P$, ${(P_y^{\prime })}_i={(P_y^{\prime })}_{i-1}·R_S$, ${(P_z^{\prime })}_i={(P_z^{\prime })}_{i-1}·R_P$. Here, ${\mathbf{P}}^{\prime }$ 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 $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 ${\mathbf{P}}_x$, and it is possible to account for polarization changes at the interface via absolute values $|T_P{|}^2,|T_S{|}^2,|R_P{|}^2,|R_S{|}^2$ of Fresnel coefficients as outlined in the previous works.¹³ This procedure influences the absolute value of polarization vectors, and, correspondingly, the weight of each MC-photon. After account for the interface influence, both ${\mathbf{P}}^{\prime }$ vectors are transformed back to the global $(x,y,z)$ reference frame. We would further use the notations $(x^{\prime },y^{\prime },z^{\prime })$ and ${\mathbf{P}}^{\prime }$ 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¹⁵ 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 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 (-{\mu }_a\sum _{i=1}^{N_j}{l}_i)$, where $0<N_j<N_{\text{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

For completeness, we also provide expressions for all intensities of the linearly polarized light:

Here, ${\mathrm{\Gamma }}_R=\frac{2}{1+\overline{{\cos }^2\theta }}$ is the Rayleigh factor derived from the optical theorem in Born approximation and $\overline{{\cos }^2\theta }$ 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 ${\mathbf{P}}_x$ and ${\mathbf{P}}_y$ are supposed to be rewritten in the new Cartesian basis with $z^{\prime }$ 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 ${\mathbf{P}}^{\prime }$ in the selected reference frame $(x^{\prime },y^{\prime },z^{\prime })$ can be used. The obtained ${\mathbf{P}}^{\prime }$ values can be directly substituted into Eqs. (23)–(30) providing appropriate intensity, Stokes, or degree of polarization values.

3.4.

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:⁶⁰

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_{\parallel }$, $V$ to $I_{\perp }$, $P$ to $I_{+45\hspace{0.17em}\mathrm{deg}}$, $M$ to $I_{-45\hspace{0.17em}\mathrm{deg}}$, $R$ to $I_R$, and $L$ to $I_L$. In terms of our model, Mueller matrix $\mathcal{M}$ of the single detected photon can be expressed as

Here, $P_x=(P_{xx},P_{xy},P_{xz})$ and $P_y=(P_{yx},P_{yy},P_{yz})$ are the real-valued vectors introduced in Sec. 3.1 and computed via Eq. (21) for incident linear polarizations $|H⟩=P_{x_0}=(\mathrm{1,0},0)$, $|V⟩=P_{y_0}=(\mathrm{0,1},0)$. Similarly, ${\mathcal{P}}_x=({\mathcal{P}}_{xx},{\mathcal{P}}_{xy},{\mathcal{P}}_{xz})$ and ${\mathcal{P}}_y=({\mathcal{P}}_{yx},{\mathcal{P}}_{yy},{\mathcal{P}}_{yz})$ are vectors computed for incident diagonal linear polarizations $|L_{+45\mathrm{deg}}⟩={\mathcal{P}}_{x_0}=(\frac{\sqrt{2}}{2},\frac{\sqrt{2}}{2},0)$ and $|L_{-45\mathrm{deg}}⟩={\mathcal{P}}_{y_0}=(\frac{\sqrt{2}}{2},-\frac{\sqrt{2}}{2},0)$. Circular polarization states $|RCP⟩$ and $|LCP⟩$ are accounted for as superpositions of $|H⟩$ and $|V⟩$ according to Eq. (4). ${\mathcal{M}}_{31}={\mathcal{M}}_{12}^{\mathrm{rot}}$ means that this element can be obtained via rotation of ${\mathcal{M}}_{12}$ by $-\pi /4$.⁶⁰ 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):