## Abstract

Tissue polarimetry has demonstrated its great potential in biomedical field presently. In this study, the polarization characteristics of red blood cell (RBC) suspensions in a back-detection geometry have been investigated with experimental measurements and Monte Carlo (MC) simulation based on Mueller matrix decomposition. It is revealed that the simulated dependence of degree of polarization (DOP) and diattenuation on the distance away from incident point is qualitatively consistent with experimental result. DOP and diattenuation decay with increasing radial distance except in the region adjacent to the incident point. Further analysis shows that the number of scattering events and the scattering angle simultaneously influence the trends of DOP and diattenuation curves in the central region.

© 2012 OSA

## 1. Introduction

Polarimetry has been shown to be effective for analyzing suspensions in fields as varied as analysis of material, biology, medicine, etc [1]. As the optical fingerprint, Mueller matrix provides complete information about the polarization properties of an optical system [2]. Mueller matrix of turbid media has been studied by many researchers [3–8]. However, the meaning of 16-element Mueller matrix is not so straightforward. Recently, increasing attentions have been devoted to polar decomposition to extract these polarization characteristics from the full Mueller matrix measurements [1, 9–11]. Tissue polarimetry, based on Mueller matrix decomposition, has demonstrated great potential in biomedical field presently. For instance, the method of Mueller matrix decomposition has been applied to image precancer in tissue [9, 10]. Another potential application is noninvasive glucose monitoring in diabetic patients. Polarimetry, based on the chiral nature of the glucose molecules and their associated optical activity, is promising as it is potentially specific to glucose. Multiple scattering is a universal phenomenon in tissue. It not only causes extensive depolarization, but also alters the polarization state of the residual polarization preserving signal. Depolarization caused by multiple scattering is the most prominent polarimetry effect in biological tissues. Birefringence is the other important polarimetry property stemming from the anisotropic organized nature of many tissues. Optical rotation in tissue arises from the presence of asymmetric optically active chiral molecules like glucose, proteins, and lipids. A diattenuator has an intensity transmission that depends on the incident polarization state of light [12]. Detailed investigation of these polarization characteristics is necessary to achieve some clinical goals via tissue polarimetry.

Analysis of polarization characteristics in backscattering plane of diffusely scattering sample has not been reported. This study is based on Mueller matrix decomposition to extract polarization characteristics of RBC suspensions. Study of polarization characteristics on backscattering surface of turbid media is also theoretically valuable in the field of tissue polarimetry.

In this study, we choose RBC suspensions as sample because it is clinically meaningful. Mueller matrices of RBC suspensions with different concentrations at backscattering surface are experimentally measured. To validate the experimental results, backscattering Mueller matrices of RBC suspensions are simulated by polarization-sensitive Monte Carlo method. To match experiments, the results for an infinitely narrow photon beam are converted to responses for Gaussian beams using Green’s convolution. Polar decomposition is applied to the backscattering Mueller matrices of RBC suspensions to extract intrinsic polarizing properties [13]. In this paper, we discussed degree of polarization (DOP) and diattenuation in the backscattering plane of RBC suspensions. And some specific properties about DOP and diattenuation are found. The findings will be helpful for diagnosis of some blood diseases with the abnormal size, number and refractive index of RBCs.

## 2. Experiment

The experimental setup is depicted in Fig. 1 . The light source was a He-Ne laser with wavelength of 632.8 nm and beam diameter of 2 mm. The light beam with linear polarization was incident on the sample with oblique angle 10 deg and collected by a CCD in normal direction. The arrangement can effectively avoid glare and does not introduce more optical elements.

In order to get 16 elements of Mueller matrix, at least 16 kinds of intensity measurements corresponding to 16 different combinations of polarizer and analyzer have to be performed. Here, however, we performed 36 measurements to remove the influence of background light for greater accuracy [14]. It was achieved by rotating polarizing block (composed of $H{W}_{1}$ and $Q{W}_{1}$) and analyzing block (composed of $Q{W}_{2}$ and $A$) at proper angles.

The components of the Mueller matrix can be calculated by the combination of six polarization state measurements as $H,V,P,Q,R$ and $L$.$H$ is linear horizontal polarization, $V$ is linear vertical polarization,$P$ is linear 45 deg polarization,$Q$ is linear −45 deg polarization,$R$ and $L$ are circular right-handed and left-handed polarizations, respectively. The Mueller matrix can be computed by the following formula:

In Eq. (1), each element of Mueller matrix is calculated from four measurements at different incident and acceptance polarization states. For example, $HV$ means that the linear vertical polarization component is analyzed with the incidence of linear horizontal polarized light.

In order to measure the spatially distributed Mueller matrices at backscattering surface, $Q{W}_{2}$ and $A$ with 25 mm optical apertures were used. The distance from upper surface of the sample to the CCD was adjusted to be 25 cm to guarantee the object-image relationship between the backscattering surface and the CCD. The processed area was about 4$\times $4 cm^{2}. The receiving angle was about 3 deg in order to make $Q{W}_{2}$ and $A$ work as expected.

Blood of New Zealand rabbit was adopted as sample, which was supplied by Jiangsu Academy of Agricultural Sciences in China. Because the shape, average size and number of RBCs in unit volume of blood of New Zealand rabbit are close to the counterparts of human being blood. RBCs were centrifugated from blood plasma. And the RBCs were mixed with physiological saline with different weight ratios to make a range of RBC suspensions.

## 3. Theory

#### 3.1 Polarization-sensitive Monte Carlo simulation

The method of polarization-sensitive Monte Carlo (PSMC) simulation has been discussed extensively by several authors [15–19]. The geometry of multiple scattering events in a polystyrene-microsphere solution is shown in Fig. 2 . A narrow pencil beam propagates downward along the z-axis into a plane-parallel slab of the scattering medium with a thickness of $L$. The incident point is (0, 0, 0) in the laboratory coordinate system (x, y, z). Photon packets are scattered in the medium by microspheres before exiting the upper or lower surface of the medium. At each scattering event, the scattering angles of the photon packets are statistically selected. During propagation, the polarization evolutions of photon packets are traced through the Stokes-Mueller formalism. The scattering histories of a large number of such photon packets are tracked as they propagate through the medium and are summed to yield the output Stokes vector. The polarization state of the incident light is changed sequentially to horizontal (Stokes vector $[\begin{array}{cccc}1& 1& 0& 0{]}^{T}\end{array}$), vertical ($[\begin{array}{cccc}1& -1& 0& 0{]}^{T}\end{array}$), + 45°($[\begin{array}{cccc}1& 0& 1& 0{]}^{T}\end{array}$) linear polarization, and circular polarization ($[\begin{array}{cccc}1& 0& 0& 1{]}^{T}\end{array}$). The Stokes vector $[{\begin{array}{cccc}I& Q& U& V]\end{array}}^{T}$ of emergent light are recorded for each given incident polarization state. The Mueller matrix is generated by performing simple algebraic manipulations using the recorded output Stokes parameters for each input state [20].

In traditional PSMC, the scattering angle and the azimuthal angle are chosen by a “rejection method”. But this method is extraordinarily time-consuming especially for RBC suspensions with high scattering coefficients (for RBC suspension with concentration of 10%, the scattering coefficient is up to 800 cm^{−1}). In our implementation, the Henyey-Greenstein (HG) phase function was used for sampling the scattering angle, and azimuthal angle was chosen uniformly between 0 and $2\pi $ as conventional Monte Carlo algorithm. It is reasonable because HG phase function is more suitable than Mie phase function for RBC due to the non-spherical shape of RBC [21]. The optical parameters were computed using Mie theory for a set of input parameters (diameter of RBC, refractive index of physiological saline, refractive index of the surrounding medium and wavelength of light). In the MC simulation, the receiving area, the area grid and the receiving angle in the backscattering plane were the same as those in the experiment. To have a satisfactory statistics, all the simulations were run with 5 × 10^{7} input photon packets.

In general Monte Carlo simulation, the response to an infinitely narrow photon beam is processed. However, all photon beams have finite size in reality. Theoretically, we can compute the response to a finite size photon beam directly by distributing the initial positions of the launched photon packets, but this method is not efficient. Fortunately, the system we are dealing with is linear and spatially invariant. Therefore, if we assume the photon beam of finite size is collimated, the response of an infinitely narrow photon beam should be a Green’s function of the tissue system. So the response of the finite size photon beam can be computed from the convolution of the Green’s function according to the profile of the finite size photon beam. One should note that the response mentioned above is Stokes vector $[{\begin{array}{cccc}I& Q& U& V]\end{array}}^{T}$. For each element of the Stokes vector, we denote the response as $C(x,y)$ and the Green’s function as $G(x,y)$. If the photon beam has the intensity profile $S(x,y)$, the response can be obtained with the use of Eq. (2)

#### 3.2 Polar decomposition of the Mueller matrices

Polar decomposition of Mueller matrix is a robust mathematical tool for interpretation of the polarization characteristics of any medium. The method decomposes an arbitrary Mueller matrix $M$ into the product of three elementary matrices representing a depolarizer, a retarder and a diattenuator.

${M}_{\Delta}$ accounts for the depolarizing effects of the medium, ${M}_{R}$ describes the effects of linear birefringence and optical activity and ${M}_{D}$ indicates the effects of linear and circular diattenuations.The process and validity of polar decomposition were first demonstrated by Lu and Chipman in clear media [13]. It is important to note that because the multiplication of the Mueller matrices is not commutative, the results of the decomposition depend upon the order in which the three elementary matrices are multiplied. It has been shown previously that the six possible decompositions can be classified into two families, depending upon the order in which the depolarization and the diattenuator matrices are multiplied. Because the family, the diattenuator matrix coming ahead of the depolarizer matrix, always leads to a physically realizable Mueller matrix [22], this order of decomposition is adopted in the present study.

The decomposition process can be referred to the Ref [13]. If ${M}_{\Delta}$, ${M}_{R}$ and ${M}_{D}$ are known, degree of polarization is defined as

${M}_{\Delta}(2,2)$ and ${M}_{\Delta}(3,3)$ are depolarization coefficients for incident horizontal (or vertical) and 45 deg (or −45 deg) linearly polarized light, and ${M}_{\Delta}(4,4)$ is the depolarization coefficients for circularly polarized incident light.${M}_{R}$ can be expressed as a combination of a matrix for a linear retarder and a circular retarder [20,23]. The values for linear retardance $\delta $ and optical rotation $\psi $ can be determined from the elements of the matrix ${M}_{R}$ as

The magnitude of diattenuation can be obtained as

## 4. Results and discussion

In our previous work [24], backscattering Mueller matrices were measured for RBC suspensions with different concentrations. Generally, it is a multiple scattering process for light propagating in a highly scattering sample (the scattering coefficient of RBC suspensions varies from 320 cm^{−1} to 800 cm^{−1}). The phase relationship by different scattering events will be randomized. Hence, the coherent backscattering effects are very weak and can be ignored. After polar decomposition for Mueller matrices, distribution of DOP and diattenuation in the backscattering plane can be obtained. As illustrated in Fig. 3
, the distribution of DOP and diattenuation is azimuthally independent. In order to discuss the dependence of DOP and diattenuation on radial distance, Fig. 4
shows the curves along central horizontal axis L_{A}.

As can be seen from Fig. 4, DOP and diattenuation decay with increasing concentrations except the central area. It is reasonable that when suspension concentration increases, scattering coefficient will increase simultaneously, therefore DOP and diattenuation decreases with increased scattering events. Also, because the light emergent at larger radial distance from the incident point experiences more scattering events, DOP and diattenuation decay with increasing distance. However, the trends of the curves in the central area seem abnormal. Number of scattering events increases with increasing radial distance. Hence, DOP should decrease with increasing radial distance according to the previous study [15]. To validate the experimental results, numerical computation based on Monte Carlo simulation was performed.

In the simulation, the shape of RBC was regarded as spherical approximately and the diameter was set as 5.56$\mu m$. The refractive indices of RBC and physiological saline were set as 1.402 and 1.33 according to Ref [25]. Scattering coefficient${\mu}_{s}$and anisotropic parameter$g$were obtained by Mie theory. In our implementation, scattering number for each photon packet exiting from the backscattering plane of RBC suspensions was recorded in order to analyze the abnormal trends in the central region of experimental curves.

The simulation curves along central horizontal axis are shown in Fig. 5 . The curves are not very smooth, especially when the radial distance gets larger. In this case, more simulation photon packets should be simulated. One can find the abnormality of DOP and diattenuation in central region of Fig. 4 also shows up in Fig. 5. In fact, there have been reports to show the similar trends in the central region by PSMC simulation, but the authors failed to explain the mechanism of the trends [20, 23]. It is noted that the abnormal range and the relation between DOP/diattenuation and position in central area in Fig. 5 do not accord with the experimental results. The inconsistency can be due to the infinitely narrow beam used in our Monte Carlo simulation. Considering the actual incident light is Gaussian beam, Green’s convolution was introduced to produce Mueller matrices for Gaussian beam from the results for infinitely narrow beam.

To match experiment, the radius of Gaussian beam was set as 2 mm. As shown in Fig. 6 , the DOP/diattenuation curves are smoother than the corresponding curves in Fig. 5.

It is evident that the radius of Gaussian beam should be considered based on the following facts: Firstly, the central abnormal range is extended for both DOP and diattenuation when radius of Gaussian beam is considered; the second, the trends of curves in the central area are greatly affected by the radius of beam especially for DOP; the third, the value of diattenuation at incident point is influenced significantly by the radius of beam. Obviously, both the abnormal ranges and the shapes of curves based on Gaussian beam are more consistent with the experimental results than the ones based on infinitely narrow beam.

The experimental results and the simulation results agree with each other qualitatively rather than quantitatively. For DOP, experimental values are evidently larger than numerical ones; as to diattenuation, the experimental value at the incident point is far larger than the simulation value. The reason causing the difference is that discrepancy exists in experimental samples and simulation samples. Experimentally, RBCs in physiological saline are not uniformly distributed on account of RBCs subsiding in the process of rotating wave-plates. Although we shook up the RBC suspensions several times during the experimental process, non-uniform distribution of RBCs in physiological saline was unavoidable. So the concentrations of experimental samples especially for supernatant liquid are lower than simulation samples, which accounts for the higher DOP in experimental results. At the other hand, the heterogeneity of RBCs in experimental samples maybe accounts for larger diattenuation at the incident point in experimental results.

To explore what factors contribute to the abnormality in central region of the curves, scattering-number-resolved analysis of polarization characteristics on the backscattering surface was performed. Figure 7 shows the DOP dependence on radial distance for photon packets experiencing different number of scattering events. It is seen that DOP is not monotonously dependent on the scattering number at P2, P3, P4 and P5 (marked in Fig. 7). At P1 and P6, DOP decreases with increasing scattering number. However, most previous studies show that DOP decays exponentially with an increasing number of scattering events in the forward-scattering plane of turbid media [15]. We suppose most large scattering events will take place in the backscattering of turbid media, when scattering number is relatively small, the events of large scattering events will influence DOP greatly. The influence gradually decreases with increasing scattering number. In fact, the similar phenomenon has been reported in Ref [26]. We argue it is reasonable because the underlying mechanism of depolarization due to multiple scattering is the scrambling of photon’s reference frame as a consequence of random sequence of scattering events in a variety of scattering directions. Hence, the abnormality of DOP in the central area should be caused by the combined effects of scattering angle and scattering number. Affected by the mechanism, it is interesting to find some specific properties from the concentration dependence in Fig. 5(a) and Fig. 6(a). DOP decreases with increasing concentration at the incident point and at the positions outside the central area. Within the transition region, the dependence relationship is not monotonous. And the transition range is sensitive to the concentration. The finding will be helpful for diagnosis of some blood diseases. Because many blood diseases usually accompany with the abnormal size, number and refractive index of RBCs. And DOP curves are also sensitive to the size and refractive index of RBCs besides the concentration of RBC suspensions. For instance, hemolysis is a common blood disease. And it will take place in plasma when osmolarity is lower than normal physiological status. Plasma osmolarity influences the size and shape of RBCs; meanwhile, the refractive index within the cell will also be altered [27]. Hence, it is expected that DOP curves will be sensitive to different hemolytic extent. This work is being investigated by our group.

A diattenuator has an intensity transmission that depends on the incident polarization state. Figure 8 shows the diattenuation dependence on the radial distance for photon packets experiencing different number of scattering events. It is seen that diattenuation increases with radial distance for the photon packets encountering the same number of scattering events; diattenuation decay monotonously with number of scattering events for the photon packets emergent from the same distance. At certain position, diattenuation is affected by many photon packets experiencing different number of scattering events. The preliminary conclusion comes out that diattenuation dependence on radial distance is also affected by scattering number and scattering angle. But the influence mechanism of scattering number and scattering angle on diattenuation is different from that on DOP.

## 5. Conclusion

We have studied the backscattering polarization characteristics of RBC suspensions experimentally and theoretically. According to the dependence of DOP and diattenuation on radial distance, we find that the abnormal variations in the central area are not caused by experimental errors, since simulation results exhibit the similar abnormality. It is found that the measured polarization characteristics and the simulation ones are in qualitative agreement even in the central region. Further analysis shows the abnormality of DOP in the central region in backscattering plane is caused by the combined effects of scattering angle and scattering number. Diattenuation is also affected by scattering angle and scattering number, but the influence mechanism is different from that on DOP. More studies should be conducted to interpret the detailed mechanism and the reasons to induce the quantitative discrepancy between experimental results and simulation results.

## References and links

**1. **F. Boulvert, G. Le Brun, B. Le Jeune, J. Cariou, and L. Martin, “Decomposition algorithm of an experimental Mueller matrix,” Opt. Commun. **282**(5), 692–704 (2009). [CrossRef]

**2. **M. H. Smith, “Interpreting Mueller matrix images of tissues,” Proc. SPIE **4257**, 82–89 (2001). [CrossRef]

**3. **J. L. Pezzaniti and R. A. Chipman, “Mueller matrix imaging polarimetry,” Opt. Eng. **34**(6), 1558–1568 (1995). [CrossRef]

**4. **A. H. Hielscher, A. A. Eick, J. R. Mourant, D. Shen, J. P. Freyer, and I. J. Bigio, “Diffuse backscattering Mueller matricesof highly scattering media,” Opt. Express **1**(13), 441–453 (1997). [CrossRef] [PubMed]

**5. **B. D. Cameron, M. J. Raković, M. Mehrübeoğlu, G. W. Kattawar, S. Rastegar, L. V. Wang, and G. L. Coté, “Measurement and calculation of the two-dimensional backscattering Mueller matrix of a turbid medium,” Opt. Lett. **23**(7), 485–487 (1998). [CrossRef] [PubMed]

**6. **M. J. Raković, G. W. Kattawar, M. B. Mehrübeoğlu, B. D. Cameron, L. V. Wang, S. Rastegar, and G. L. Coté, “Light backscattering polarization patterns from turbid media: theory and experiment,” Appl. Opt. **38**(15), 3399–3408 (1999). [CrossRef] [PubMed]

**7. **P. Yang, H. Wei, G. W. Kattawar, Y. X. Hu, D. M. Winker, C. A. Hostetler, and B. A. Baum, “Sensitivity of the backscattering Mueller matrix to particle shape and thermodynamic phase,” Appl. Opt. **42**(21), 4389–4395 (2003). [CrossRef] [PubMed]

**8. **M. R. Antonelli, A. Pierangelo, T. Novikova, P. Validire, A. Benali, B. Gayet, and A. De Martino, “Mueller matrix imaging of human colon tissue for cancer diagnostics: how Monte Carlo modeling can help in the interpretation of experimental data,” Opt. Express **18**(10), 10200–10208 (2010). [CrossRef] [PubMed]

**9. **X. Li and G. Yao, “Mueller matrix decomposition of diffuse reflectance imaging in skeletal muscle,” Appl. Opt. **48**(14), 2625–2631 (2009). [CrossRef] [PubMed]

**10. **J. Chung, W. Jung, M. J. Hammer-Wilson, P. Wilder-Smith, and Z. Chen, “Use of polar decomposition for the diagnosis of oral precancer,” Appl. Opt. **46**(15), 3038–3045 (2007). [CrossRef] [PubMed]

**11. **N. Ghosh, M. F. G. Wood, S. H. Li, R. D. Weisel, B. C. Wilson, R. K. Li, and I. A. Vitkin, “Mueller matrix decomposition for polarized light assessment of biological tissues,” J Biophotonics **2**(3), 145–156 (2009). [CrossRef] [PubMed]

**12. **N. Ghosh and I. A. Vitkin, “Tissue polarimetry: concepts, challenges, applications, and outlook,” J. Biomed. Opt. **16**(11), 110801 (2011). [CrossRef] [PubMed]

**13. **S.-Y. Lu and R. A. Chipman, “Interpretation of Mueller matrices based on polar decomposition,” J. Opt. Soc. Am. A **13**(5), 1106–1113 (1996). [CrossRef]

**14. **M. Itoh, M. Yamanari, Y. Yasuno, and T. Yatagai, “Polarization characteristics of multiple backscattering in human blood cell suspensions,” Opt. Quantum Electron. **37**(13-15), 1277–1285 (2005). [CrossRef]

**15. **X. Wang, L. V. Wang, C. W. Sun, and C. C. Yang, “Polarized light propagation through scattering media: time-resolved Monte Carlo simulations and experiments,” J. Biomed. Opt. **8**(4), 608–617 (2003). [CrossRef] [PubMed]

**16. **D. Côté and I. A. Vitkin, “Robust concentration determination of optically active molecules in turbid media with validated three-dimensional polarization sensitive Monte Carlo calculations,” Opt. Express **13**(1), 148–163 (2005). [CrossRef] [PubMed]

**17. **J. C. Ramella-Roman, S. A. Prahl, and S. L. Jacques, “Three Monte Carlo programs of polarized light transport into scattering media: part I,” Opt. Express **13**(12), 4420–4438 (2005). [CrossRef] [PubMed]

**18. **J. C. Ramella-Roman, S. A. Prahl, and S. L. Jacques, “Three Monte Carlo programs of polarized light transport into scattering media: part II,” Opt. Express **13**(25), 10392–10405 (2005). [CrossRef] [PubMed]

**19. **http://omlc.ogi.edu/software/polarization/

**20. **N. Ghosh, M. F. G. Wood, and I. A. Vitkin, “Polarimetry in turbid, birefringent, optically active media: A Monte Carlo study of Mueller matrix decomposition in the backscattering geometry,” J. Appl. Phys. **105**(10), 102023 (2009). [CrossRef]

**21. **J. Q. Lu, P. Yang, and X. H. Hu, “Simulations of light scattering from a biconcave red blood cell using the finite-difference time-domain method,” J. Biomed. Opt. **10**(2), 024022 (2005). [CrossRef] [PubMed]

**22. **J. Morio and F. Goudail, “Influence of the order of diattenuator, retarder, and polarizer in polar decomposition of Mueller matrices,” Opt. Lett. **29**(19), 2234–2236 (2004). [CrossRef] [PubMed]

**23. **S. Manhas, M. K. Swami, P. Buddhiwant, N. Ghosh, P. K. Gupta, and J. Singh, “Mueller matrix approach for determination of optical rotation in chiral turbid media in backscattering geometry,” Opt. Express **14**(1), 190–202 (2006). [CrossRef] [PubMed]

**24. **X. Wang, L. Yang, J. Lai, and Z. Li, “Polar decomposition applied to light back-scattering by erythrocyte suspensions,” Proc. SPIE **8192**, 81924T, 81924T-6 (2011). [CrossRef]

**25. **A. Roggan, M. Friebel, K. Dörschel, A. Hahn, and G. Müller, “Optical properties of circulating human blood in the wavelength range 400-2500 nm,” J. Biomed. Opt. **4**(1), 36–46 (1999). [CrossRef]

**26. **X. Guo, M. F. G. Wood, and A. Vitkin, “Monte Carlo study of pathlength distribution of polarized light in turbid media,” Opt. Express **15**(3), 1348–1360 (2007). [CrossRef] [PubMed]

**27. **M. Friebel, J. Helfmann, and M. C. Meinke, “Influence of osmolarity on the optical properties of human erythrocytes,” J. Biomed. Opt. **15**(5), 055005 (2010). [CrossRef] [PubMed]