1.

Introduction

The use of polarimetric techniques has received considerable recent attention for biomedical diagnosis.¹ This is because the polarization properties of scattered light contain additional diagnostic information on tissue, which cannot be obtained from the use of unpolarized light. For example, collagen (a structural protein) and other fibrous tissues exhibit linear birefringence due to their oriented structures. In various kinds of tissue abnormalities, the structural (and functional) properties of collagen and other fibrous structures show distinct changes. Measurement of the linear birefringence of tissue may thus find useful applications in noninvasive diagnosis of various pathological diseases.^{2, 3, 4, 5} Glucose is another important tissue constituent that possesses circular birefringence due to its asymmetric chiral structure. Its presence in tissue leads to rotation of the plane of linearly polarized light about the axis of propagation (known as optical rotation). Measurements of optical rotation may offer an attractive approach for noninvasive monitoring of tissue glucose levels, and several studies have investigated the polarization properties of scattered light from chiral turbid medium. ^{6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16} In addition to these applications, polarimetric measurements are also finding use for contrast enhancement in biomedical imaging.^{17, 18, 19, 20, 21}

The use of the polarimetric techniques for biological tissue characterization is, however, severely hampered by multiple scattering effects. In a turbid medium such as tissue, the incident polarized light becomes strongly depolarized due to multiple scattering, the magnitude and rate of depolarization being determined by a large number of parameters such as the concentration, size, shape, refractive index of the scatterers, detection geometry, and by the incident light’s state of polarization.^{22, 23, 24, 25} Polarization modulation and synchronous detection methods, therefore, have been developed to extract the weak fraction of the polarization preserving component of scattered light from a turbid medium and analyze it for determination of the diagnostically relevant polarization properties such as linear birefringence and optical activity.^{7, 9, 10, 26} In particular, this approach is being investigated to determine of the concentration of optically active glucose in human tissue. However, studies have shown that even when depolarization due to multiple scattering is taken into account, the relatively small rotation of the polarization vector arising from optical activity in the medium is hidden by the much larger change in the orientation angle of the polarization vector due to scattering.^{7, 27} In addition to the scattering effects, interpretation of optical rotation using the conventional Stokes polarimetry-based approach is further confounded by the effects of linear birefringence and its complex orientation in real tissues.

Because the Mueller matrix (which transforms the Stokes vector describing the polarization properties of incident light to the Stokes vector of scattered light) description contains complete information about all the polarization properties of a medium,^{28, 29} measurements of full Mueller matrix elements may be used to address these issues. However, when many optical polarization effects are simultaneously occurring in the sample (as is the case for tissue), the resulting elements of the Mueller matrix reflect several “lumped” effects, thus hindering their unique interpretation. A methodology that allows the extraction of the individual “intrinsic” material polarimetry characteristics from the lumped Mueller matrix description is thus needed. The independent constituent polarization parameters—linear retardance ( $\delta$ , difference in phase between two orthogonal linear polarization, and its orientation angle $\theta$ ), circular retardance or optical rotation ( $\psi$ , difference in phase between right and left circularly polarized light), diattenuation ( $d$ , differential attenuation of orthogonal polarizations for both linear and circular polarization states), and depolarization coefficient ( $\Delta$ , linear and circular)—can be efficiently extracted by decomposing the measured Mueller matrix into the “basis” Mueller matrices of a retarder, diattenuator, and depolarizer.³⁰ The individual polarization parameters obtained through this process should provide more intrinsic information on the underlying morphological features and physiological state of tissue and may thus offer a more fundamental (less empirical) way to understand its structure and behavior.

Indeed, preliminary studies have shown that Mueller matrix decomposition can be used to extract the pure optical rotation component corresponding to the concentration of chiral substances by decoupling the contribution of scattering-induced rotation of the polarization vector that primarily arises due to the contribution of singly (or weakly) scattered photons in a scattering medium.³¹ In a highly scattering medium such as thick tissue, however, the contribution of singly scattered photons is expected to be small and the majority of the detected photons will be multiply scattered. Further, unlike the ideal situation, in actual tissue, all the polarization effects are exhibited simultaneously and not in an ordered, sequential manner. To realize the potential of the Mueller matrix decomposition approach for quantification and interpretation of the useful polarization parameters (linear retardance and optical rotation in particular) in real tissues, it is necessary to test the validity of the decomposition process on multiply scattering media exhibiting simultaneous linear birefringence, optical activity, and depolarization. A quantitative understanding of the confounding effects of multiple scattering on the values for the polarization parameters estimated through the decomposition process is also essential.

We therefore investigate these aspects both theoretically and experimentally. The effects of multiple scattering, the propagation path of scattered photons, and detection geometry on the values for $\delta$ and $\psi$ were studied systematically using polarization-sensitive Monte Carlo simulations of the propagation of polarized light in birefringent, chiral, multiply scattering media having varying optical properties. In the simulations, the simultaneous effects of linear birefringence and optical activity on the photons between successive scattering events were incorporated using the Jones N-matrix formalism.³² In the corresponding experimental studies, Mueller matrices were recorded from realistic tissue phantoms exhibiting controllable scattering, linear birefringence, and optical activity. The values for the polarization parameters estimated through the decomposition process were compared with the known optical properties of the phantoms. The results of these theoretical investigations and the experimental studies are reported in this communication.

The organization of this paper is as follows. Section 2 provides a description of the experimental systems and methods as well as the optical phantoms used in this study. The forward polarization-sensitive Monte Carlo simulations of propagation of polarized light in multiply scattering media and the inverse process for polar decomposition of the resulting Mueller matrix are described in Sec. 3. The results of decomposition of the experimental and the Monte Carlo–generated Mueller matrices are presented in Sec. 4 to validate the method and to elucidate the observed trends. Section 5 concludes the paper with a discussion on the implications of these results and potential use of the Mueller matrix decomposition methodology in diagnostic photomedicine.

2.

Experimental Methods and Materials

The Mueller matrices were recorded from the tissue phantoms using a photoelastic modulation (PEM)–based polarimeter.^{7, 32} A schematic of the experimental setup is shown in Fig. 1 . Unpolarized light at $632.8\mathrm{nm}$ from a He-Ne laser (LHRR-1200M, Research Electro-Optics, Boulder, Colorado) was first passed through a mechanical chopper (C) operating at a frequency $f_c=500\mathrm{Hz}$ and then through a linear polarizer (P1) and a removable quarter wave plate (WP1). The linear polarizer and quarter wave plate combination enabled generation of the four required input polarization states, $0\mathrm{deg}$ (Stokes vector ${\left[1100\right]}^T$ ), $45\mathrm{deg}$ (Stokes vector ${\left[1010\right]}^T$ ), $90\mathrm{deg}$ (Stokes vector ${[1-100]}^T$ ) linear polarization, and circular polarization (Stokes vector ${\left[1001\right]}^T$ ). The four set of Stokes vectors $\left[IQUV\right]$ of light exiting the sample in selected detection directions were measured for each of the four input polarization states. For this study, we chose forward detection geometry; other detection geometries (including backscattering) are ongoing. The detection optics consists of another removable quarter wave plate (WP2) with its fast axis oriented at $-45\mathrm{deg}$ . When WP2 was in place, we could measure the Stokes parameters $Q$ and $U$ (linear polarization),⁷ and when it was removed, we could measure the Stokes parameter $V$ (circular polarization).³² The forward scattered light from the sample then passed through a photoelastic modulator (PEM, IS-90, Hinds Instruments, Hillsboro, Oregon, operating at a frequency $f_p=50\mathrm{kHz}$ ) with its fast axis oriented at $0\mathrm{deg}$ . The light finally passed through an analyzer (P2) oriented at $45\mathrm{deg}$ . The resulting modulated intensity of the scattered light was collected using a pair of lenses (detection area of $1{\mathrm{mm}}^2$ and acceptance angle $\sim 18\mathrm{deg}$ ) and was relayed to an avalanche photodiode (APD) detector (C5460, Hamamatsu Photonics, Hamamatsu City, Japan). The detected signal was sent to a lock-in amplifier (SR 830, Stanford Research Systems, Sunnyvale, California) with reference input of the amplifier toggled between the chopper $\left(500\mathrm{Hz}\right)$ and the PEM controller ( $50\mathrm{kHz}$ and harmonics).

Fig. 1

Schematic of the experimental setup. C, mechanical chopper; P1 and P2, polarizers; WP1 and WP2, removable quarter wave plates, PEM, photoelastic modulator; APD, avalanche photodiode; ${\mathrm{f}}_{\mathrm{c}}$ and ${\mathrm{f}}_{\mathrm{p}}$ modulation frequencies of the chopper and PEM, respectively.

The Stokes vector of the output light for each of the four input polarization states were measured using the detected signal at the first harmonic of the chopper frequency and the first and the second harmonics of the PEM frequencies.^{7, 32} The Mueller matrix was generated using standard relationships between its 16 elements and the measured output Stokes parameters $\left[IQUV\right]$ for each of the four input polarization states.^{28, 29}.

The optical phantoms used in this study were polyacrylamide gels having strain-induced linear birefringence, sucrose-induced optical activity (magnitude of optical activity was $\chi =1.96\mathrm{deg}{\mathrm{cm}}^{-1}$ , corresponding to $1\mathrm{M}$ concentration of sucrose) and polystyrene microspheres–induced scattering (mean diameter $D=1.40\mu \mathrm{m}$ , refractive index $n_s=1.59$ ).^{32, 33} The refractive index of polyacrylamide was measured to be $n_m\sim 1.393$ at the wavelength $\lambda =632.8\mathrm{nm}$ . The phantom preparation details are described in Ref. 32. The values for scattering coefficient $\left({\mu }_s\right)$ of the phantoms at $632.8\mathrm{nm}$ were calculated using the known concentration of the scatterers and Mie theory.³⁴ The value for anisotropy parameter ( $g$ , average cosine of scattering angle) of the scattering phantom at $632.8\mathrm{nm}$ was calculated³⁴ to be $g=0.95$ .

To apply controllable strain to produce linear birefringence, one end of the polyacrylamide phantoms (dimension of $1\times 1\times 4\mathrm{cm}$ ) was clamped to a mount and the other end to a linear translational stage. The phantoms were stretched along the vertical direction (along the long axis of the sample) to introduce varying linear birefringence with its axis along the direction of strain. Mueller matrices were recorded in the forward detection geometry from the phantoms with varying extension (varying degree of strain-induced birefringence).

3.

Theory

4.

Results and Discussion

Mueller matrices were experimentally recorded in the forward scattering geometry from the polyacrylamide phantoms having varying degrees of strain-induced birefringence ( $\text{extension}=0\text{to}4\mathrm{mm}$ applied along the vertical direction), varying scattering coefficients ( ${\mu }_s=0$ , 3, 4, and $6{\mathrm{mm}}^{-1}$ ), and sucrose-induced optical activity (concentration of $1\mathrm{M}$ corresponding to magnitude of optical activity $\chi =1.96\mathrm{deg}{\mathrm{cm}}^{-1}$ ). Measurements were also performed on a clear (nonscattering) polyacrylamide phantom at different extensions to obtain the calibration curve for the strain-induced birefringence with extension. Measurements were performed three times from each of the experimental phantoms and the mean values were taken to construct the Mueller matrices. Standard deviations in the values for the normalized Mueller matrix elements were found to be in the range 0.0015 to 0.0050.

Table 1 shows the experimental Mueller matrix and the decomposed (following the procedure described in Sec. 3.2) depolarization $\left({\mathrm{M}}_{\Delta }\right)$ , retardance $\left({\mathrm{M}}_{\mathrm{R}}\right)$ , and diattenuation $\left({\mathrm{M}}_{\mathrm{D}}\right)$ matrices for the nonscattering phantom ( ${\mu }_s=0{\mathrm{mm}}^{-1}$ , thickness $t=10\mathrm{mm}$ ) with extension of $4\mathrm{mm}$ and having a concentration of sucrose of $1\mathrm{M}$ . The values for all the retrieved polarization parameters [ $R,d$ ( $\delta$ and $\theta$ ), $\psi$ , and $\Delta$ , determined from the equations in Secs. 3.2.1, 3.2.2, 3.2.3] are also listed in the table. As expected, the nonscattering phantom does not show any appreciable depolarization $(\Delta \sim 1)$ . The value for $\delta$ is $\sim 1.345\mathrm{rad}$ for an extension of $4\mathrm{mm}$ and the estimated orientation angle of the axis of linear retardance $(\theta =89.05\mathrm{deg})$ is reasonably close to the actual orientation of the linear birefringence axis ( $\theta =90\mathrm{deg}$ for strain applied along the vertical direction). The value for $\delta$ increased from $0.088\text{to}1.345\mathrm{rad}$ for varying extension of $0\text{to}4\mathrm{mm}$ . Due to presence of optical rotation, the estimated value for the total retardance parameter $R$ [that represents the combined effect of linear retardance and optical rotation, see Eq. 10] is marginally larger $(R=1.346\mathrm{rad})$ than the corresponding value for $\delta$ . The estimate for $\psi$ $\left(1.77\mathrm{deg}\right)$ , however, is, lower than that expected for concentration of sucrose of $1\mathrm{M}$ ( $\psi =1.96\mathrm{deg}$ for a path length of $10\mathrm{mm}$ ). The contraction of the phantom due to longitudinal stretching reduces the effective path length of transmitted light, contributing to the observed lower value of $\psi$ . The reduction in path length at $4\text{-}\mathrm{mm}$ extension is $\sim 0.33\mathrm{mm}$ (path length of $9.67\mathrm{mm}$ as compared with $10\mathrm{mm}$ for the phantom with no strain) using the value for Poisson ratio $(\sim 0.33)$ of polyacrylamide.^{32, 40} However, even when accounting for this thickness reduction, the new expected $\psi$ $\left(1.90\mathrm{deg}\right)$ is still higher than the one derived from polar decomposition. This difference is possibly due to an uncertainty in the concentration of sucrose during the fabrication of the phantom. It is also seen from Table 1 that decomposition of Mueller matrix yields a small but nonzero value for the diattenuation parameter $d$ . The presence of a small amount of dichroic absorption (at $\lambda =632.8\mathrm{nm}$ ) due to anisotropic alignment of the polymer molecules in the polyacrylamide phantom may contribute to this slight nonzero value of diattenuation.

Table 1

(a) The experimentally recorded Mueller matrix and the decomposed basis matrices for a nonscattering phantom (μs=0mm−1) with strain-induced birefringence (extension of 4mm ) and a concentration of sucrose of 1M (χ=1.96degcm−1) . (b) The values for the different polarization parameters extracted from the decomposed matrices. M(1.0000−0.02290.00270.0058−0.01860.9956−0.03610.0318−0.01290.03920.2207−0.96560.00140.02800.97060.2231) MΔ(1.00000000.00410.996900−0.007000.99150−0.0019000.9966)MR(1.000000000.9988−0.03620.028100.03930.2227−0.974100.03200.97420.2239)MD(1.0000−0.02290.00270.0058−0.02291.0000−0.0000−0.00010.0027−0.00000.99970.00000.0058−0.00010.00000.9997)

Table 2 displays the experimental Mueller matrix, the decomposed basis matrices, and the retrieved polarization parameters for a nonbirefringent $(\text{extension}=0\mathrm{mm})$ , achiral (concentration of $\text{sucrose}=0\mathrm{M}$ , $\chi =0\mathrm{deg}{\mathrm{cm}}^{-1}$ ), turbid phantom ( ${\mu }_s=3{\mathrm{mm}}^{-1}$ , $t=10\mathrm{mm}$ , optical thickness $\tau ={\mu }_s\times t=30$ , $g=0.95$ ). Ideally, the Mueller matrix of this phantom should model an isotropic depolarizer. The values for $\delta$ , $d$ , and $\psi$ obtained following decomposition are indeed found to be low for this nonbirefringent, achiral sample. The small nonzero value for $\delta$ (0.087) may arise due to the small amount of strain present during clamping of the sample to the mount and translational stage. Another possible cause, the scattering-induced linear retardance,³¹ which arises primarily from singly (or weakly) scattered photons, is not expected to contribute to the nonzero value for $\delta$ , because for this scattering sample having a value for optical thickness $\tau =30$ , multiply scattered photons are the dominant contributor to the detected photons in the forward detection geometry. Further, in agreement with previous reports,^{23, 24, 25} for this anisotropic scattering sample ( $g=0.95$ , prepared using large-sized scatterer with diameter $D⩾\lambda$ ), depolarization of circularly polarized light is weaker than depolarization of linearly polarized light [ ${\mathrm{M}}_{\Delta }(4,4)>{\mathrm{M}}_{\Delta }(2,2)$ and ${\mathrm{M}}_{\Delta }(3,3)$ ]. This was previously shown to arise due to the fact that in scattering media composed of large scatterers $D⩾\lambda ,g⩾0.7$ , the randomization of helicity (which is the major cause of depolarization of circularly polarized light) is considerably reduced due to predominantly forward scattering, leading to weaker depolarization of circularly polarized light.^{23, 24, 25} Note that the depolarization of linearly polarized light is almost independent of the orientation angle (horizontal, vertical, and $45\mathrm{deg}$ ) of the incident linear polarization vector [depolarization coefficients ${\mathrm{M}}_{\Delta }(2,2)\sim {\mathrm{M}}_{\Delta }(3,3)$ ]. This makes sense because in an isotropic, multiply scattering medium, randomization of the field vector’s direction is the primary cause of depolarization of incident linearly polarized light. The linear depolarization coefficients of the resultant depolarizing Mueller matrix in such a situation will be the same for any orientation of the incident linear polarization vector.⁴¹

Table 2

(a) The experimentally recorded Mueller matrix and the decomposed basis matrices for a nonbirefringent (extension=0mm) , achiral (concentration of sucrose=0M , χ=degcm−1 ), turbid phantom (optical thickness τ=30 , μs=3mm−1 , thickness t=10mm , g=0.95 ). (b) The values for the different polarization parameters extracted from the decomposed matrices. M(1.0000−0.01900.00300.0114−0.01830.7454−0.00250.00190.00130.00850.7378−0.06480.00340.00370.07470.8530) MΔ(1.0000000−0.00420.7468000.000000.73940−0.0065000.8566)MR(1.000000001.0000−0.0076−0.000500.00750.9961−0.087400.00120.08740.9962)MD(1.0000−0.01900.00300.0114−0.01900.9999−0.0000−0.00010.0030−0.00000.99980.00000.0114−0.00010.00000.9998)

The experimental Mueller matrix, the decomposed basis matrices, and the retrieved polarization parameters for the birefringent $(\text{extension}=4\mathrm{mm})$ , chiral ( $1\mathrm{M}$ sucrose, $\chi =1.96\mathrm{deg}{\mathrm{cm}}^{-1}$ ), and turbid phantom ( ${\mu }_s=3{\mathrm{mm}}^{-1}$ , $g=0.95$ ) are shown in Table 3 . Note the complicated nature of the Mueller matrix of this complex phantom and relatively unequivocal nature of the three basis matrices derived from the decomposition process. The elements of the depolarization matrix $\left({\mathrm{M}}_{\Delta }\right)$ and the resulting value for the depolarization coefficient ( $\Delta$ ) are similar to the corresponding values for the pure depolarizing phantom having the same scattering coefficient ${\mu }_s$ (Table 2, difference in the value for $\Delta <2\%$ ). The decomposition process thus successfully decouples the depolarization effects (due to multiple scattering) and retrieves the values for the parameters $\delta$ (and $\theta$ ) and $\psi$ . However, the estimate for $\psi$ is found to be larger than that for the clear phantom having the same concentration of sucrose ( $\psi =2.04\mathrm{deg}$ for the scattering phantom as compared with $1.77\mathrm{deg}$ for the clear phantom). This is due to the increase in the propagation path of light in a multiply scattering medium. The average photon path length of light exiting the scattering sample in the forward direction was calculated to be $11.7\mathrm{mm}$ from Monte Carlo simulation (optical parameters of the medium: ${\mu }_s=3{\mathrm{mm}}^{-1}$ , diameter of scatterer $D=1.40\mu \mathrm{m}$ , $n_s=1.59$ , $n_m=1.393$ , $g=0.95$ ). For the calculation of average photon path length, the thickness of the scattering medium was taken to be $9.67\mathrm{mm}$ as per the Poisson ratio discussion earlier. The value for $\psi$ calculated using the optical rotation of the nonscattering phantom ( $\psi =1.77\mathrm{deg}$ from Table 1) and the Monte Carlo–generated average photon path length comes out to be reasonably close to that estimated through decomposition of the Mueller matrix recorded from the birefringent, chiral, turbid phantom (calculated value of $\psi =2.07\mathrm{deg}$ as compared with $\psi =2.04\mathrm{deg}$ estimated for the turbid phantom).

Table 3

(a) The experimentally recorded Mueller matrix and the decomposed basis matrices for a birefringent (extension=4mm) , chiral (concentration of sucrose=1M , χ=1.96degcm−1 ), turbid ( μs=3mm−1 , g=0.95 ) phantom. (b) The values for the different polarization parameters extracted from the decomposed matrices. M(1.0000−0.03120.0029−0.0066−0.02140.7678−0.03700.0204−0.00550.02300.1043−0.77350.00140.03900.79720.1920) MΔ(1.00000000.00280.755200−0.010200.768900.0016000.8454)MR(1.000000000.9984−0.05130.024200.03330.1844−0.982300.04600.98150.1858)MD(1.0000−0.03120.0029−0.0066−0.03121.0000−0.00000.00010.0029−0.00000.9995−0.0000−0.00660.0001−0.00000.9995)

Although the estimated value for $\delta$ of the turbid phantom is larger than that for the clear phantom ( $\delta =1.384\mathrm{rad}$ for the turbid phantom as compared with $1.345\mathrm{rad}$ for the clear phantom having a path length of $10\mathrm{mm}$ ), the value is lower than that one would expect for average photon path length of $11.7\mathrm{mm}$ $(\delta =1.345\times 1.17=1.574\mathrm{rad})$ . This is illustrated in Fig. 2 where the estimates for $\delta$ and $\psi$ of the chiral $(\chi =1.96\mathrm{deg}{\mathrm{cm}}^{-1})$ phantoms having varying degrees of strain-induced birefringence (extension of $0\text{to}4\mathrm{mm}$ ) are displayed. Results are shown for both the clear $({\mu }_s=0{\mathrm{mm}}^{-1})$ and the turbid $({\mu }_s=3{\mathrm{mm}}^{-1})$ phantoms. The increase in the value for $\psi$ as a result of increased average photon path length in the turbid phantom as compared with the clear phantom can be seen quite clearly. The gradual decrease in the value for $\psi$ with increasing longitudinal stretching is consistent with resulting lateral contraction of the phantom, reducing the effective path length. In contrast to $\psi$ , the expected increase in $\delta$ as a result of increased average photon path length in the turbid phantom is not seen. The results of studies on phantoms having other values of scattering coefficients ( ${\mu }_s=4$ and $6{\mathrm{mm}}^{-1}$ , data not shown) also yielded similar trends to that observed for the phantom with ${\mu }_s=3{\mathrm{mm}}^{-1}$ .

Fig. 2

The values for linear retardance $\left(\delta \right)$ and optical rotation $\left(\psi \right)$ estimated from the decomposition of experimentally recorded Mueller matrices from the chiral $(\chi =1.96\mathrm{deg}{\mathrm{cm}}^{-1})$ phantoms having varying degrees of strain-induced birefringence (extension of $0\text{to}4\mathrm{mm}$ ). Results are shown for both clear $({\mu }_s=0{\mathrm{mm}}^{-1})$ and turbid $({\mu }_s=3{\mathrm{mm}}^{-1})$ phantoms.

To understand these trends further and to investigate the effect of multiple scattering, propagation path, and detection geometry on $\delta$ and $\psi$ parameters, Monte Carlo simulations were carried out.

The Monte Carlo simulations were run with input optical parameters chosen to mimic the optical properties of the phantoms used in the experimental studies. The photon collection geometry was chosen to have a detection area of $1{\mathrm{mm}}^2$ and an acceptance angle of $20\mathrm{deg}$ around the forward directed ballistic beam. Standard deviations (originating from random statistical errors) in the different elements of the Monte Carlo–generated Mueller matrices were found to be in the range 0.0004 to 0.0025. Table 4 shows the Mueller matrix, the decomposed basis matrices, and the retrieved polarization parameters for a nonbirefringent, achiral, turbid medium (linear birefringence $\Delta n=0$ , $\chi =0°{\mathrm{cm}}^{-1}$ , $\tau =30$ , ${\mu }_s=3{\mathrm{mm}}^{-1}$ , $t=10\mathrm{mm}$ ) composed of monodisperse spherical scatterers ( $D=1.40\mu \mathrm{m}$ , $n_s=1.59$ , $n_m=1.393$ , $g=0.95$ ). In agreement with the experimental results (Table 2), polarization preservation of circularly polarized light is found to be larger than linearly polarized light [ $M_{\Delta }(4,4)>M_{\Delta }(2,2)$ and $M_{\Delta }(3,3)$ ] for this forward scattering $(g=0.95)$ turbid medium. The elements of the decomposed depolarization matrix ${\mathrm{M}}_{\Delta }$ are also close to the experimental Mueller matrix results (difference in the value for $\Delta \sim 4\%$ ) for the phantom having similar scattering coefficient (small differences may arise due to slight mismatch in the collection geometry employed in the experiments and in the simulation, and possibly due to subtle errors in concentration of scatterers in the phantom). Further, as expected, the values for $\delta$ , $d$ , and $\psi$ of this nonbirefringent, achiral, turbid medium are also small. The small value of $\delta$ $\left(0.0023\mathrm{rad}\right)$ for this nonbirefringent turbid medium confirms that scattering-induced linear retardance should not contribute to the experimentally observed slightly larger value for $\delta$ $\left(0.087\mathrm{rad}\right)$ for the phantom having the same scattering coefficient of ${\mu }_s=6{\mathrm{mm}}^{-1}$ .

Table 4

(a) The Monte Carlo simulation–generated Mueller matrix and the decomposed basis matrices for a nonbirefringent, achiral, turbid medium ( Δn=0 , χ=0degcm−1 , τ=30 , μs=3mm−1 , t=10mm , g=0.95 ). (b) The values for the different polarization parameters extracted from the decomposed matrices. M(1.0000−0.0025−0.0026−0.0024−0.00080.77870.00010.0009−0.00200.00030.77670.00340.00020.0004−0.00040.8809) MΔ(1.00000000.00120.7788000.000000.776600.0024000.8810)MR(1.000000001.0000−0.00010.000300.00011.00000.00230−0.0003−0.00231.0000)MD(1.0000−0.0025−0.0026−0.0024−0.00251.00000.00000.0000−0.00260.00001.00000.0000−0.00240.00000.00001.0000)

Figure 3 shows the variation of the parameters $d$ , $\delta$ , $\Delta$ , and $\psi$ of transmitted light as a function of optical thickness $\tau$ (varying ${\mu }_s$ and fixed path length of $t=10\mathrm{mm}$ ) for the nonbirefringent, achiral, turbid medium. As expected, depolarization decreases with increasing value of $\tau$ . In the absence of chiral molecules, the values for $\psi$ are low at all the optical thickness values of the achiral, turbid medium. Further, both the parameters $\delta$ and $d$ are $<0.02$ for the range of $\tau =10$ to 80. As discussed previously, the scattering-induced linear retardance and diattenuation arises primarily from singly (or weakly) scattered photons and the magnitude is greater at large scattering angles. Because for detection around the forward direction, the contribution of such photons is negligible, the scattering-induced retardance and diattenuation are weak in this geometry.

Fig. 3

Variation of the parameters $\delta$ , $d$ , $\psi$ , and $\Delta$ (estimated from the decomposition of Monte Carlo–generated Mueller matrices) as a function of optical thickness $\tau$ (varying ${\mu }_s$ and fixed thickness of $t=10\mathrm{mm}$ ) for nonbirefringent $(\Delta n=0)$ , achiral $(\chi =0\mathrm{deg}{\mathrm{cm}}^{-1})$ , turbid medium composed of monodisperse spherical scatterers ( $D=1.40\mu \mathrm{m}$ , $n_s=1.59$ , $n_m=1.393$ , $g=0.95$ ). The results are shown for transmitted light collected from the ballistic beam position at the exit face of the medium.

Table 5 presents the Monte Carlo–generated Mueller matrix, the decomposed basis matrices, and the retrieved polarization parameters for a birefringent, chiral, turbid medium ( $\Delta n=1.36\times {10}^{-5}$ , corresponding to $\delta =1.35\mathrm{rad}$ for a path length of $t=10\mathrm{mm}$ , optical activity $\chi =1.96\mathrm{deg}{\mathrm{cm}}^{-1}$ , $\tau =30$ , ${\mu }_s=3{\mathrm{mm}}^{-1}$ , $g=0.95$ ). The axis of linear birefringence was along the vertical $(y,\theta =90\mathrm{deg})$ direction in the simulation. The values for $\delta$ and $\psi$ estimated following the decomposition process are reasonably close to those of the experimental phantom having similar properties (Table 3). Slightly larger disagreement in $\psi$ $(\sim 14\%)$ possibly arises due to an uncertainty in the concentration of sucrose during the fabrication of the phantom (as noted previously). Incorporation of the path length reduction due to longitudinal stretching ( $t=9.67\mathrm{mm}$ as compared with $t=10\mathrm{mm}$ ) in the simulations did not result in significant changes in the estimated values for $\delta$ and $\psi$ (reduction in $\delta$ and $\psi$ were $⩽3\%$ ). As for the nonbirefringent, achiral, turbid medium (Table 4), this birefringent, chiral, turbid medium prepared using large-sized scatterers (diameter $D⩾\lambda$ ) also exhibits weaker depolarization of circularly polarized light than linearly polarized light $[M_{\Delta }(4,4)>M_{\Delta }(2,2),M_{\Delta }(3,3)]$ . The elements of ${\mathrm{M}}_{\Delta }$ and the resulting value for $\Delta$ are also similar to the corresponding values for the nonbirefringent, achiral, turbid medium with the same scattering properties.

Table 5

(a) The Monte Carlo simulation–generated Mueller matrix and the decomposed basis matrices for a birefringent, chiral, turbid medium ( Δn=1.36×10−5 , χ=1.96degcm−1 , τ=30 , μs=3mm−1 , t=10mm , g=0.95 ). (b) The values for the different polarization parameters extracted from the decomposed matrices. M(1.0000−0.00390.0022−0.00970.00120.77420.0310−0.0393−0.0061−0.04290.1234−0.7949−0.0009−0.02510.81460.1983) MΔ(1.00000000.00380.775800−0.014200.78360−0.0009000.8595)MR(1.000000000.99790.0412−0.04900−0.05610.1948−0.97920−0.03080.98000.1967)MD(1.0000−0.00390.0022−0.0097−0.00391.0000−0.00000.00000.0022−0.00000.9999−0.0000−0.00970.0000−0.00001.0000)

Note that at an even larger value of the input linear birefringence $(\Delta n>1.36\times {10}^{-5})$ , the decomposed circular depolarization coefficient $\left[M_{\Delta }(4,4)\right]$ was found to be lower than the corresponding value for the nonbirefringent, turbid medium (the trend can be seen from Tables 4, 5 also). Because for this forward scattering medium $(g\sim 0.95)$ , depolarization of circularly polarized light is weaker than depolarization of linearly polarized light,^{23, 24, 25} while propagating through such a birefringent (with axis of birefringence along the vertical direction), turbid medium, a greater fraction of incident circularly polarized light gets transferred to $45\text{-}\mathrm{deg}$ linear polarization as compared with transfer of $45\text{-}\mathrm{deg}$ linear polarization to circular polarization. This also led to a difference in the value for the depolarization coefficient for $45\text{-}\mathrm{deg}$ linearly polarized light and horizontally (or vertically) polarized light $[M_{\Delta }(3,3)\ne M_{\Delta }(2,2)]$ . This suggests that even in the presence of strong multiple scattering, for a turbid medium that exhibits linear birefringence, the depolarization of linearly polarized light may depend on the orientation angle of the incident linear polarization vector.

To gain quantitative understanding of the dependence of the $\delta$ and $\psi$ parameters on the propagation path of multiply scattered photons, Fig. 4 shows the Monte Carlo results for the variation of $\delta$ and $\psi$ as a function of the thickness of the same birefringent, chiral, scattering medium. The top axis shows the calculated average photon path length of light exiting the scattering medium. The value for $\psi$ increases with increasing average photon path length. However, the absolute values for $\psi$ were marginally lower compared to calculations using the linear relationship ( $\psi =\chi \times \text{average}$ photon path length). It should be noted here that the average path length has contributions from both the polarization-preserving and the depolarized photons. The fact that the propagation path of the polarization-preserving photons (which would show experimentally detectable optical rotation) is shorter than the average photon path length⁴² should account for the observed lower value for $\psi$ .

Fig. 4

Variation of $\delta$ and $\psi$ (estimated from the decomposition of Monte Carlo–generated Mueller matrices) as a function of the thickness of a birefringent ( $\Delta n=1.36\times {10}^{-5}$ , corresponding to a value of $\delta =1.35$ for path length of $10\mathrm{mm}$ ), chiral $(\chi =1.96\mathrm{deg}{\mathrm{cm}}^{-1})$ , turbid medium ( ${\mu }_s=3{\mathrm{mm}}^{-1}$ , $\tau =30$ , $g=0.95$ ). The Monte Carlo–calculated average photon path length of light exiting the scattering medium is shown on the top axis.

It is interesting that although the value for $\delta$ increases with increasing average photon path length, the increase is much slower as compared with $\psi$ . A similar trend was also observed in the experimental studies. This is because in a turbid medium, the forward scattered light does not travel in a straight line but rather along a curved path, the curvature being controlled by the values for ${\mu }_s$ and $g$ . Though such paths do not have any influence on $\psi$ (as long as the average photon path length is the same), the net value for $\delta$ would be significantly affected because the linear birefringence has a specific orientation axis ( $\theta =90\mathrm{deg}$ in the simulations and in the experiments).

To understand this, Mueller matrices were generated for transmitted light collected at the exit face of the birefringent ( $\Delta n=1.36\times {10}^{-5}$ and $0.68\times {10}^{-5}$ ), chiral $(\chi =1.96\mathrm{deg}{\mathrm{cm}}^{-1})$ , turbid medium ( $\tau =30$ , ${\mu }_s=3{\mathrm{mm}}^{-1}$ , $g=0.95$ ) at different spatial positions away from the position of the ballistic beam along the horizontal ( $x$ axis, perpendicular to the direction of linear birefringence) and vertical ( $y$ axis, parallel to the direction of linear birefringence) directions. The variations for the parameter $\delta$ as a function of the distance from the ballistic beam position (along both $x$ and $y$ axis) at the exit face are displayed in Fig. 5 . At spatial positions along the $x$ axis, $\delta$ increases with increasing distance from the center of the ballistic beam (i.e., with increasing average photon path length). In contrast, for detection at spatial positions along the $y$ axis, the value for $\delta$ shows gradual decrease with increasing distance from the center of the ballistic beam. This is because a larger component of the photon propagation path is along the axis of birefringence leading to a reduction in net linear retardance $\delta$ (because propagation along the direction of the birefringence axis does not yield any retardance). Because such differences in the photon propagation path for the two different detection geometries should have no influence on $\psi$ , the estimates for $\psi$ were found to be identical for similar detection positions either along the $x$ or the $y$ axis at the exit face of the medium (data not shown). It thus appears that the value for $\delta$ is influenced by the detection position and the orientation angle of the axis of linear birefringence in a birefringent, turbid medium.

Fig. 5

Variation of $\delta$ (estimated from the decomposition of Monte Carlo–generated Mueller matrices) of transmitted light as a function of distance from ballistic beam position at the exit face of a birefringent, chiral $(\chi =1.96\mathrm{deg}{\mathrm{cm}}^{-1})$ , turbid medium ( $\tau =30$ , ${\mu }_s=3{\mathrm{mm}}^{-1}$ , $g=0.95$ ). The results are shown for two different values of birefringence ( $\Delta n=1.36\times {10}^{-5}$ and $0.68\times {10}^{-5}$ ). The axis of linear birefringence was kept along the vertical ( $y$ axis) direction and the results are shown for transmitted light collected at different spatial positions along the horizontal ( $x$ axis, perpendicular to the direction of the axis of birefringence) and vertical ( $y$ axis, parallel to the direction of the axis of birefringence) direction at the exit face of the medium. The Monte Carlo–calculated average photon path length of light exiting the scattering medium is shown on the top axis. The inset shows the detection geometry.

To examine this further, in Fig. 6 , we show the variation of $\delta$ as a function of the orientation angle $\left(\theta \right)$ of the axis of linear birefringence for the birefringent, chiral, turbid medium ( $\Delta n=1.36\times {10}^{-5}$ , $\chi =1.96\mathrm{deg}{\mathrm{cm}}^{-1}$ , $\tau =30$ , ${\mu }_s=3{\mathrm{mm}}^{-1}$ , $g=0.95$ ). Results are shown for three different detection positions at the exit face of the medium, around the position of the ballistic beam [coordinate (0,0)], at $3\mathrm{mm}$ away from the position of the ballistic beam along the horizontal $\left(x\right)$ [coordinate (3,0)] and vertical $\left(y\right)$ [coordinate (0,3)] axes, respectively. The value for $\delta$ is observed to vary considerably with a change in the birefringence orientation angle for off-axis detection [at spatial positions (3,0) or (0,3)]. As noted previously, the minimum values for $\delta$ are found at spatial positions along the direction of the birefringence axis. Due to the symmetry of position coordinates with respect to the direction of birefringence axis, the two values of $\delta$ [for spatial positions (3,0) and (0,3)] are also observed to converge for the value of $\theta =45\mathrm{deg}$ . In contrast, for detection around the position of the ballistic beam [coordinate (0, 0)], the estimates for $\delta$ do not show any significant variation with a change in the orientation angle of the axis of birefringence. Thus for simultaneous determination of the intrinsic values for the parameters $\delta$ and $\psi$ of a birefringent, chiral, turbid medium in the forward scattering geometry, detection around the direction of propagation of the ballistic beam may be preferable.

Fig. 6

Variation of $\delta$ (estimated from the decomposition of Monte Carlo–generated Mueller matrices) as a function of the orientation angle $\theta$ (with respect to the horizontal axis) of the axis of linear birefringence for a birefringent, chiral, turbid medium ( $\Delta n=1.36\times {10}^{-5}$ , $\chi =1.96\mathrm{deg}{\mathrm{cm}}^{-1}$ , $\tau =30$ , ${\mu }_s=3{\mathrm{mm}}^{-1}$ , $g=0.95$ ). Results are shown for three different detection positions at the exit face of the medium, detection around the position of the ballistic beam [position coordinate (0,0)], detection at spatial positions $3\mathrm{mm}$ away from the ballistic beam position along the horizontal $\left(x\right)$ [coordinate (3,0)] and vertical $\left(y\right)$ [position coordinate (0,3)] axis, respectively.

The presented results demonstrate that decomposition of the Mueller matrix can be used for simultaneous determination of the intrinsic parameters $\delta$ (and $\theta$ ) and $\psi$ of a birefringent, chiral, turbid medium. For conceptual and practical reasons, the extension of the developed methodology to backward detection geometry is warranted. However, this will be more complex, because in the backward detection geometry, the contribution of backscattered (large angle scattered) photons would be greater and thus the interplay of the scattering-induced linear retardance and diattenuation with the intrinsic values for $\delta$ and $\psi$ would be more coupled. Future studies will investigate this aspect in detail and will concentrate on developing an appropriate strategy to minimize the scattering-induced artifacts for simultaneous determination of $\delta$ and $\psi$ in the backscattering geometry using Mueller matrix decomposition.

5.

Conclusions

We have investigated the validity of Mueller matrix decomposition methodology to extract the individual intrinsic polarimetry characteristics from a multiply scattering medium exhibiting simultaneous linear birefringence and optical activity. Decomposition of experimental Mueller matrices recorded from tissue phantoms exhibiting controllable scattering, linear birefringence, and optical activity yielded satisfactory estimates for linear retardance $\delta$ and optical rotation $\psi$ . Reliable estimates for linear and circular depolarization coefficients were also obtained. The effects of multiple scattering, the propagation path of scattered photons, and detection geometry on $\delta$ and $\psi$ were further investigated using polarization-sensitive Monte Carlo simulations. The simulation results were found to corroborate the experimental findings and demonstrate that the Mueller matrix decomposition methodology can successfully extract the individual polarization characteristics of a medium that exhibits simultaneous linear birefringence, optical activity, and multiple scattering effects.

The ability of this technique to determine the constituent contributions of individual processes should prove valuable in diagnostic photomedicine. The value for linear retardance estimated using this approach would represent its average value over the volume of the turbid medium (tissue) probed by the polarization-preserving photons⁴³ and could monitor changes in tissue birefringence as a signature of tissue abnormality and/or response to therapy. Note that in actual tissues due to larger value of ${\mu }_s$ and due to presence of absorption, the sampling depth of the polarization-preserving photons will be considerably reduced (typically a few millimeters) as compared with that of the phantoms used in this study. As a representative application, we are currently investigating the use of the decomposition-derived linear retardance for monitoring myocardial tissue regeneration following stem cell treatments of myocardial infarction. Measurement (in the forward detection geometry) through $\sim 1\text{-}\mathrm{mm}$ thick ex vivo myocardial tissue samples are yielding encouraging results and show promise for the use of this approach for monitoring of stem-cell–based treatments of myocardial infarction. These results will be reported in a separate publication. Further, the decomposition method may also find useful applications in quantifying the optical rotations due to blood glucose in diabetic patients, potentially impacting the (currently unsolved) noninvasive glucose-monitoring problem. The results presented here demonstrate the ability to quantify optical activity (albeit at higher than physiological levels) despite depolarization and other effects in multiply scattering media, using the decomposition methodology. Ordinarily, the small chirality-caused optical rotation is swamped by much larger scattering-induced rotation, and it is the decomposition method that enables the glucose-related signal to be derived and quantified from the Mueller matrix results. The application of this approach to real tissues with lower blood glucose levels ( $\sim 3\text{to}30\mathrm{mM}$ ), however, will, require further refinements of the highly sensitive Mueller matrix measurements to detect small changes in the matrix elements corresponding to the physiological glucose levels. This and the combined use of the Mueller matrix decomposition approach and spectroscopic polarimetry (coupled with chemometric data analysis to isolate the rotation due to glucose from that caused by other chiral constituent confounders in tissue) is currently under investigation in our laboratory.