## Abstract

The concentration determination of optically active species in moderately turbid suspensions is studied both experimentally and with a validated three-dimensional polarization-sensitive Monte Carlo model. It is shown that the orientation of the polarization of the scattered light exhibits a strong dependence on exit position in the side or backscattered directions, but not in the forward direction. In addition, it is shown that the increased path length of photons due to multiple scattering in a 1 cm cuvette filled with forward-peaked scatterers (anisotropy around 0.93) increases the optical rotation by up to 15%, but only for scattering coefficients under 30 cm^{-1}, after which it decreases again. It is concluded that in order to avoid systematic errors in concentration determination of optically-active molecular species in turbid samples, the scattered light in the forward direction should be used.

©2005 Optical Society of America

## 1. Introduction

Recently, polarimetric measurements in tissues have gained importance by their ability to provide information about ballistic and near-ballistic photons with relative simplicity. For instance, polarization imaging has been used to enhance spatial resolution in turbid samples[1] and to increase contrast for surface imaging[2] by removing the depolarized intensity that consists of multiply scattered photons. By looking at the linear polarization orientation of the remaining polarized photons, one can obtain further information about tissues. Since chiral molecules rotate the polarization of light and tissue conserves a fraction of the light polarization, various experimental schemes have been proposed to relate the rotation of the linear polarization of an incident light beam to the concentration of optically active molecular species in a sample by using scattered light from the forward-, side- or backward directions[3, 4, 5]. A case of potential interest here is the determination of glucose concentration in turbid solutions and its application to blood glucose monitoring[6]. In a *clear* solution, the concentration *C* of an optically active molecular species can be obtained from the measurement of the optical rotation with:

where *α* is the optical rotation of the linear polarization, [*α*${]}_{\lambda}^{T}$
is the known rotatory power of the molecular species at a given wavelength *λ* and temperature *T*, and 〈*L*〉 the photon path length in the solution[7]. The typical rotation due to 1 mM of optically active molecules through 1 cm of clear solution is on the order of 1 millidegree, which is readily measurable[5]. Since parts of the eye are transparent and since glucose levels in the eye can be related to that in blood[8], polarimetric measurements of the aqueous humour have been proposed for monitoring blood glucose levels[9, 10]. Although promising, difficulties remain due to layered structure of the cornea that causes both birefringence and multiple reflections. It has therefore been suggested that optically thick multiply scattering tissue could be sampled since it partially preserves the polarization of light and its polarization rotation can be measured[4, 5]. However, difficulties arise in turbid tissue when trying to extract the concentration of optically active molecules from optical rotation measurements. One is that the polarization properties of the beam may be modified upon scattering in ways that are unrelated to glucose concentration. Another problem is that the paths photons travel in tissue are not unique and are obviously dependent on the detection geometry (measurements in reflection or transmission will detect photons that have travelled very different distances). Since the requirements for reflection or transmission experiments are very different, a theoretical study is first needed in order to determine an optimal geometry for measuring this rotation.

The problem of light propagation in tissue can be formulated using transport theory for polarized particles. To make the transport equations analytically solvable, polarization is usually neglected and certain simplifying assumptions are made to obtain the diffusion equation for highly scattering media[11, 12]. Although elegant, the approximations are not always appropriate and a numerical approach is often required. The Monte Carlo technique can be used to numerically solve the transport equations, whereby random photon paths are launched in tissue and, upon completion, the macroscopic properties of interest are computed from the statistical average of the photon properties. This technique is the most general approach currently available in tissue optics (albeit at the price of long computation time). Initial implementations of the Monte Carlo algorithm were only for intensity calculations, and the most common is the Monte Carlo Multi Layer (MCML) implementation of Wang *et al*.[13] for layered geometries, which has been validated extensively with intensity measurements (reflectance, absorbance, transmittance, point spread function). However, both the diffusion equation and MCML neglect the polarization of light. New Monte Carlo implementations and modifications of existing ones have been presented in the literature recently and incorporate the polarization modelling. Most of these do not study the rotation of linear polarization and instead focus on the surviving polarization fraction of light after propagating through phantoms and tissues[14, 15, 16, 17, 18]. The works of Mehrübe [19] and Wang[3] look at optical rotation, but only in the backscattering geometry of layered structures. These polarization-sensitive implementations are not available for others to use and, most importantly, have had very limited validation. It must not be underestimated that, although the formalism itself is not complicated, implementation is very challenging and prone to errors because of the multiple three-dimensional transformations of reference frames required when tracking the polarization of the photon. In the present article, the shortcomings of other implementations are addressed, and a polarization-sensitive Monte Carlo implementation (referred to as Pol-MC) that considers any three-dimensional geometry is presented, validated quantitatively with experimental data, and made publicly available[20].

Using this validated model, the dependence of the optical rotation on the details of the detection geometry and its relation to concentration of chiral molecules in turbid media is investigated.

The outline of this paper is as follows: first, sections 2 and 3 briefly introduce the Monte Carlo algorithm and the experimental setup used for the measurements. Section 4 validates the Pol-MC simulation results with experiments and section 5 discusses the dependence of the optical rotation measurements on the detection geometry and average photon path length in the context of determination chiral molecule concentration. Section 6 contains a brief conclusion.

## 2. Methods: simulations

The numerical work is performed using the Monte Carlo technique, a technique that has been discussed extensively in the past by several authors with [14, 16, 18, 19] and without[13] polarization. The present implementation builds upon the work of Kaplan[18] and Jaillon[16], except for the interface crossing algorithm (which here is applicable to any three-dimensional structure) and for the computation of the detected Stokes vector components (which here is computed from the transmitted intensities through linear polarizers). The reader is referred to Appendix A for mathematical details of the various transformations. The computer implementation (source code in ANSI C++) can be obtained from the Internet[20].

The samples into which photons propagate are characterized by a set of surface elements which are unambiguously determined by two vectors (**a** and **b**) that span the plane of the surface element and a normal *n̂* pointing outwards (see Fig. 6 in the appendix). The photon position, propagation direction and polarization (Stokes vector with coefficients *I*,*Q*,*U*,*V*, see appendix A.1) are initialized at the entrance of the sample. An interaction event occurs at a distance *d* (based on the probability distribution exp(-*µ*_{t}*d*), where the extinction coefficient *µ*_{t}
is the sum of the scattering and absorption coefficients *µ*_{s}
+*µ*_{a}
), by which the photon is moved if it does not cross an interface[13]. According to the Mueller matrix of the scatterer and the polarization state of the photon at that point, a scattering plane and a scattering angle are statistically sampled[16, 18]. The Stokes vector of the photon is first expressed in this new scattering plane, and then transformed by the application of the Mueller matrix of the scattering event (see appendices A.2 and A.3). In the present work, the Mueller matrix is calculated from Mie theory[21]. If the photon encounters an interface, the photon is moved to the interface and its Stokes vector is transformed into the Fresnel *ŝ* and *p̂* reference frame (see appendix A.2). The intersection point is calculated from three-dimensional linear vector algebra (see appendix A.4). The probability of transmission or reflection is calculated from the polarization-dependent Fresnel coefficient and the polarization state of the photon at that point[16]. The photon’s Stokes vector and propagation direction are modified according to reflection or transmission (see appendix A.5). If the photon is reflected, the propagation inside the object continues. If the photon is transmitted through an interface and outside of the object, then the measured Stokes vector is computed from the sum and differences of intensities through linear and circular polarizers (see appendix A.6). Since no interference effects are considered, the macroscopic values of the Stokes coefficients are obtained from the sum of the Stokes vector components of incoming photons, after transformation through linear polarizers. The macroscopic averages are tabulated according to the exit position and the cosine of the exit angle that the propagation vector of the photon makes with the normal to the interface. One quantity of importance is the orientation *γ* of the linear polarization of the beam, which is obtained by:

where *Q* and *U* are two coefficients of the Stokes vector[22]. The rotation *α* of the polarization is simply the difference between the final and initial orientations (*α*=*γ*-*γ*
_{∘})

To obtain convergence for intensity Monte Carlo calculations in optically thick tissue (*e*.*g*.within 1% of the asymptotic values), one typically requires millions of photons since the statistical noise can be assumed to scale with 1/√*N*, where *N* is the number of *detected* photons. The number of detected photons can be estimated from the effective extinction coefficient of the medium and the thickness of the sample when interested in transmitted photons (in another implementation[18], the number of detected photons is monitored dynamically to determine the number of photons required for the computation). In polarization-sensitive calculations, there are two properties in addition to the number of detected photons to consider. First, only a fraction *β*_{L}
(*β*_{C}
) of the photons (typically 5% to 95%) remain linearly (circularly) polarized and therefore the Stokes parameters *Q*,*U* (*V*) are smaller than the intensity by a factor *β*_{L}
(*β*_{C}
). Second, polarization effects appear as a small amplitude modulation of the Stokes vector coefficients. For instance, a polarization rotation of *α* radians appears as a sin2*α* amplitude modulation of *Q* and/or *U*. For *α*=1°, this is a 3% amplitude modulation, but for physiologically-relevant glucose concentrations — where the rotation is expected to be of the order of one millidegree — this is decreased by a factor of 1000. Hence, the required number of photons for good polarization-sensitive calculations is orders of magnitude larger than for intensity-only calculations. Typically, 10^{9} to as much as 10^{11} photons are used, with the largest numbers of photons required for high scattering (*µ*_{s}
>50 cm^{-1}). Although excessive for a single workstation, this can readily be split in smaller calculations and dispatched to a cluster of computers in order to complete the calculation in a reasonable amount of time (from a few hours to as much as three days on a 50 node cluster). With the addition of other variance reduction techniques in the computer implementation, the time required for computation could be reduced significantly (although the validity of such approximations would have to be experimentally verified).

## 3. Methods: experiments

The polarimetric measurements are performed with the experimental setup presented previously[5] as shown schematically in Fig 1. The unpolarized laser beam of approximately 10 mW of power at 633 nm goes through a mechanical chopper, operating at a frequency *f*_{c}
=(2*π*)^{-1}
*ω*
_{c}=200 Hz, then through a polarizer *P*
_{1}(*θ*
_{1}) oriented with its axis making an angle *π*/8 (22.5°) with the horizontal, and a photo-elastic modulator (PEM, Hinds, PEM-90) modulating its retardation *δ*(*t*) at frequency *f*_{p}
=(2*π*)^{-1}
*ω*
_{p}=50 kHz according to a sinusoidal function *δ*(*t*)=*δ*
_{∘} sin*ω*
_{p}
*t* with its axis horizontal. The amplitude of the retardation modulation *δ*
_{∘} is user-specified. After the sample, the two orthogonal linear polarizations are separated with a polarizing beam splitter and sent to a pair of balanced detectors that give the real time difference of the two intensities (essentially the *Q* coefficient of the Stokes vector). Based on the same analysis as in previous work[5], it has been shown that from the measurements with a lock-in amplifier at frequencies *f*_{c}
and ${f}_{2{f}_{p}}$, the rotation *α* can be obtained with:

where *Q*_{fc}
and ${Q}_{{f}_{2{f}_{p}}}$ are the measured signals at modulation frequencies *f*_{c}
and ${f}_{2{f}_{p}}$, *ℱ*(2*f*
_{p})/*ℱ*(*f*
_{c}) is the ratio of the frequency responses of the detector at those two frequencies (0.56 in our case), and *J*
_{2}(*δ*
_{∘}) is the second-order Bessel function of the first kind (at *δ*
_{∘}=3, *J*
_{2}(*δ*
_{∘}) is maximum and has a value of 0.48).

The turbid samples are suspensions of polystyrene spheres (*n*_{s}
=1.59, density *ρ*=1.05 g/cm^{3}) in distilled water (*n*=1.33) with weight fraction *f*_{w}
(weight of microspheres to weight of water) ranging from 0 to 0.1%. There is negligible absorption at the wavelength used. The spheres have a diameter of 2*r*=1.0 *µ*m or 1.4 *µ*m, and Mie theory allows the calculation of the scattering efficiency *Q*
_{scat} (*i*.*e*., the scattering cross-section is *Q*
_{scat}
*πr*
^{2}) and anisotropy *g* (the average cosine of the scattering angle)[21, 23]. For non-absorbing polystyrene microspheres of diameter 1.0 and 1.4 *µ*m in pure water (no glucose), the scattering efficiency *Q*
_{scat} is 2.59 and 3.56, and the anisotropy is 0.917 and 0.930 respectively. The solutions with glucose have a nominal concentration of 1 M, which is used in order to induce a rotation easily measurable. The actual concentration is determined at the time the solutions are made by measuring the mass of glucose being added to the known water volume. The addition of glucose changes the index of refraction of water by 0.027 per molar of glucose[24], which in turn changes the anisotropy and scattering efficiency. With 1.2Mof glucose and 1.4 *µ*m spheres, the anisotropy is 0.943 and the scattering efficiency is 3.206. The scattering coefficient is obtained from *µ*_{s}
=3*Q*
_{scat}
*f*_{w}
*ρ*
_{∘}/4*rρ*(with *ρ*
_{∘} the density of water 1 g cm^{-3}), and the reduced scattering coefficient *µ′*_{s}
=(1-*g*)*µ*_{s}
. For very dense suspensions, this linear relation between *µ*_{s}
and *f*_{w}
does not hold[25, 26], but the suspensions in the present work are dilute enough for the relation to hold. Sets of suspensions with varying scattering coefficients are prepared via serial dilution to ensure well defined scattering coefficients. The results here are also applicable to physiological concentrations of glucose since the optical rotation is proportional to the concentration of glucose even at low concentrations[5].

## 4. Validation and results

The validation of the polarization-sensitive Monte Carlo implementation is presented here, before it is used in conjunction with experiments to examine the dependence of rotation on detection geometry and scattering properties of the sample. First, the transmittance, reflectance and absorbance for layer geometries are computed and compared with the MCML implementation fromWang *et al*.Table 1 shows that good agreement with MCML is obtained for an infinite layer of thickness 1 cm and index 1.33 in air with moderate scattering (*µ*_{s}
=10 cm^{-1}, *µ*_{a}
=0.1 cm^{-1}) and forward peaked scattering (Henyey-Greenstein phase function with anisotropy *g*=0.93, from 1.4 *µ*m-diameter polystyrene scatterers in a medium with index 1.33). The number of photons launched is 10^{6}. In the case where the Henyey-Greenstein phase function is used in Pol-MC virtually identical results are obtained (within uncertainties for number of photons detected) as expected since the algorithms are identical in this particular case (polarization is not considered because the Henyey-Greenstein distribution is for intensity only). When the complete Mie scattering phase function is used in Pol-MC, the results are slightly different. The slight discrepancies are due to the polarization-sensitive Fresnel coefficients used at interfaces and to the different distributions of scattering angles (since the Mie phase function is used in Pol-MC for sampling the scattering angles and computing the polarization changes instead of

the Henyey-Greenstein phase function used in MCML). The exact details of the phase function are known to produce slightly different results for transmittance and reflectance calculations even with the same anisotropy[27]. Table 2 demonstrates the effects of finite geometry on diffuse

reflectance and transmittance measurements. As the width of the sample is reduced, the transmittance and reflectance decrease due to the loss of light through the sides of the sample. This stresses the importance of finite sample size in calculations since, as can be seen from the table, even when the sample is 4 times as wide as it is deep, there is still a small discrepancy with the infinite geometry for intensity measurements. In the case where 1 cm path length cuvettes are used, as is the case in most experiments, one must consider the finite size of the cuvette.

Next, polarimetric calculations of the surviving polarization fraction and optical rotation are validated against experimental measurements. Calculations and measurements of the surviving polarization fraction *β*_{L}
≠(*I*
_{90°}-*I*
_{0°})/(*I*
_{90°}+*I*
_{0°}) for turbid suspensions of polystyrene microspheres (without glucose) are performed as a function of scattering coefficient *µ*_{s}
in a 1 cm cuvette, where *I*
_{0°} and *I*
_{90°} are the intensities measured with analyzers at 0° and 90° just before the detector. Figure 2a shows the surviving polarization fraction *β*_{L}
of a beam initially linearly polarized along 90° as a function of the scattering coefficient for 1 *µ*m microspheres with no absorption (phase function with anisotropy *g* of 0.917). The imaged area is approximately 5 mm in diameter with an acceptance angle of 15° for the apertures and lenses used. Very good agreement is obtained between the experiments and calculations, without any adjustable parameters in the calculations.

Additional validation has been done with data presented by other researchers, but an extensive survey is not presented here. In general, good qualitative and moderate quantitative agreement is reached. For example, Fig. 2b shows the measurements of Ghosh *et al*.[28] and calculations using the Monte Carlo model presented here (varying concentrations of polystyrene microspheres of 1.072 *µ*m diameter in water, *g* 0.923, cuvette thickness of 0.5 and 1.0 cm, no absorption). The imaged area (not mentionned in Ref. [28]) is taken to be 1 cm in diameter with an acceptance angle of 20°, as inferred from the text. The quantitative discrepancy observed on Fig. 2b can potentially be explained as follows. First, small uncertainties in scattering coefficients due to aggregation of microspheres or evaporation of water from the solutions over time can be important since the surviving polarization fraction follows an exponential-like (*i*.*e*., rapidly-varying) behavior with the scattering coefficient. Second, the area and the acceptance angle over which the signals are integrated when calculating *β*_{L}
influence the result signficantly. Often the acceptance angle of the detection system is adjusted with an aperture, but this also affects the imaged area. This makes comparison with results in the literature difficult, but nevertheless, when reasonable assumptions are made, as was done here, general agreement is usually obtained.

Experiments were performed to validate optical rotation predictions. Calculations and measurements of the polarization rotation in the forward-scattering direction as a function of the scattering coefficient are shown on Fig. 3 for water suspensions of 1.4 *µ*m diameter microspheres with 1.2 M of glucose (anisotropy from Mie scattering phase function is 0.943). The agreement between calculations and experiments is very good, except at larger scattering coefficients where uncertainties in the experimental measurements are quite large due to the low intensities reaching the detector. The rotation calculated from the average photon path length, as obtained from the Monte Carlo calculation and Eq. (1) if it held, is also shown on the graph.

Initially, the increase of photon path length leads to an increase in optical rotation (by approximately 10–15%). However, at high turbidity (*µ*_{s}
≳30 cm^{-1}), although the path length 〈*L*〉 is still increasing, no increase in optical rotation is observed, and it decreases back towards its original value, as further addressed in the Discussion.

Using the validated model, Monte Carlo calculations have been performed to investigate the dependence of the rotation on the detection geometry in forward-, side- and back-scattering geometry for a turbid, cubic, aqueous suspension of polystyrene microspheres containing 1 M of glucose. The curves on Fig. 4 show the orientation (as calculated from Eq. (2)) of a 45°-polarized beam on the forward, side and back face of a cube of 1 cm filled with polystyrene microspheres of 1.0, 1.4 and 1.8 *µ*m diameter (scattering coefficient *µ*_{s}
=20 cm^{-1}, *µ*_{a}
=0 cm^{-1}), with a corresponding anisotropy calculated from their Mie phase function of 0.917, 0.930 and 0.922 respectively (background index of 1.33 assumed for solution). The rotation is obtained by subtracting 45° from the orientation. In the forward direction, the optical rotation is not strongly dependent on the position where the measurement is taken, as can be seen from the uniformity of the rotation across the front face of the sample. For a 1 cm path length, the expected rotation is 0.80° in clear media, close to the 0.88° calculated by Pol-MC in turbid media. This approximately 10% turbidity-induced increase in polarization rotation is in agreement with the trend shown on Fig. 3 and will be addressed further in the Discussion below. The orientation of the beam in the side- and back-scattering directions varies greatly as a function of the exit position and scatterer size. It varies between from 47° to 45° in the lateral direction, with a sharp increase when approaching the backscattering direction, at which point it is 135°.

## 5. Discussion

The experimental determination of the concentration of optically active molecules from rotation measurements is only possible if one can measure the rotation accurately and then relate it to the concentration. If the measurement strongly depends on the exact position where the measurement is done, concentrations are difficult to extract reliably. Figure 4 shows that the actual orientation of the beam polarization changes significantly with the position where the detection is performed. This is in agreement with what has been observed in other work, where a non-zero optical rotation of the polarization is observed even in the absence of glucose in solution when the measurements are done in the side- or back-scattering geometry[4, 19, 30, 31]. This rotation due to scattering can be understood with the help of Fig. 5, which shows the orientation of the polarization of the beam after a single scattering event (polystyrene spherical scatterer of 1.4 µm diameter in water, *g*=0.93). Upon scattering by optically inactive spheres, the orientation of the electric field of the photon cannot be modified in the forward or backward direction because there is no break of symmetry. In the forward direction, the field remains oriented in the same direction. In the back scattering direction however, because the direction of the photon changes, the actual orientation of the polarization is rotated by 90° in the frame of reference of the photon, as can be seen on the inset of Fig. 5. Since the orientation of the electric field will not change discontinuously, there must be a continuous change from the original 45° orientation to 135°, and this occurs along the side face of the cube. This is similar in spirit to the helicity flip of circularly polarized light discussed by Schmitt[1], but applied to linearly polarized light. As can be seen from the other calculations with slightly different parameters also shown on Fig. 4, the exact shape of the curve depends on the wavelength, anisotropy factor, scattering coefficient, and diameter of scatterers (since they have different scattering phase functions). Hence, any scheme for concentration determination based on the measurement in the side or back-scattered light will suffer from very high sensitivity to experimental conditions (position of detector) in addition to having a strong dependence on the details of the scattering parameters. The use of a reference image[19] with zero concentration to extract the rotation is difficult for in vivo applications. For such reasons, only in the forward scattering geometry will the rotation be directly related to glucose concentration, and be minimally affected by scattering artifacts.

The amount of polarization rotation a beam undergoes is related to the concentration of optically active molecules, but its relation to average path length is not linear, as illustrated by Fig. 3. In the forward geometry, the results shown on Fig. 3 indicate that the photons travelling a much longer distance due to scattering do not contribute significantly to optical rotation because they become part of the depolarized background. The linear dependence of optical rotation on the optical path length predicted from Eq. (1) for the modelled turbid chiral media therefore applies to average photon path length values no more than ≈15% larger than the physical thickness of the sample (in this case where *µ*_{s}
≈30 cm^{-1} or µ*′*_{s}
=2.1 cm^{-1} in a 1 cm thick sample with scattering anisotropy *g*=0.943). As typical scattering values for biological tissues in the near infrared are 3–5 times larger, this suggests that 1 cm-thick samples may be beyond that linear regime. Thus, although at moderate turbidity one must take into consideration the slight increase in photon path length when extracting the concentration from the optical rotation as per Eq. (1), at higher turbidity (as in tissue) the optical rotation returns to its original smaller value since only photons propagating in the straightforward direction contribute to rotation. Hence, a good estimate of the concentration from measurements of the optical rotation may be obtained from Eq. (1) with the physical thickness of the sample as the path length, with an additional small correction that can be determined from Fig. 3 if the scattering coefficient can be estimated.

## 6. Conclusion

In conclusion, we have introduced a three-dimensional polarization-sensitive Monte Carlo implementation that is available publicly[20]. The model was validated experimentally with surviving polarization fraction and optical rotation measurements in tissue phantoms of forward-peaked scatterers (g around 0.93), and was also compared to another Monte Carlo implementation and polarization measurements. The model was used to show the effect detection geometry on optical rotation measurements in the context of concentration determination of optically active molecules in turbid samples. It has been confirmed that in the side- or back-scattering direction, the orientation of the linear polarization of the light is changing rapidly, but this is related not only to the presence or absence of optically active molecules, but rather is mainly due to the scattering itself. On the other hand, the optical rotation of the light polarization on the illumination axis in the forward scattering direction is not strongly dependent on the position of the detectors nor on the scatterer size, and thus its measurement can be used to infer the concentration of known molecular species in the turbid media. Contrary to what might have been expected from Eq. (1), the effect of increased average path length on optical rotation in multiply scattering regime is not linear. It has been shown to increase the optical rotation of the beam by no more than 15% (in a 1 cm thick sample with scattering anisotropy *g*=0.943), and can be accounted for with an estimate of the scattering coefficient. Future work will investigate the exact dependence of this increase in rotation on absorption, scattering coefficient and anisotropy, but preliminary calculations for a variety of parameters indicate that the increase is always within 10 to 20%.

The authors would like to acknowledge Josh Grimes for technical help and Cesar Rendon for comments on early drafts of this manuscript. The computing facilities were provided by The Shared Hierarchical Academic Research Computing Network (SHARCNET). Financial support was provided by Natural Sciences and Engineering Research Council (NSERC) of Canada.

A. Implementation details

Below for completeness are listed the implementation details of the polarization-sensitive Monte Carlo algorithm described in section 2. For complete details on the Monte Carlo algorithm itself, the reader is referred to the work of Wang[13] (without polarization) and others[14, 15, 16, 17, 18] (with polarization). The work here differs from previous work in the way it treats interface crossing (described in section A.4) and detection (described in section A.6). The source code which implements the polarization-sensitive, three-dimensional Monte Carlo algorithm can be downloaded from the Internet[20].

*A.1. Definitions*

Polarized light and its interaction with matter is best described using Stokes parameters and Mueller calculus. The state of polarization of the light, with respect to a chosen set of orthonormal axes *ê*
_{‖} and *ê*
_{⊥}, is given by a Stokes vector **S** of the form:

where the same notation and definitions as in Bohren and Huffman[32] are used. We have:

with *E*
_{‖} and *E*
_{⊥} the complex electric field components (with their complex conjugate ${E}_{\Vert}^{*}$ and ${E}_{\perp}^{*}$, and *i*≠√-1). *I* represents the intensity of the beam, *Q* and *U* represent the linear polarization (respectively in the frame of reference made of *ê*
_{‖} and *ê*
_{⊥}, and in another frame rotated by 45° with respect to the former), and *V* represents the circular polarization. We define detector operators: when one applies a detector operator onto a Stokes vector S, the scalar value of that particular Stokes parameter is obtained. Mathematically, detector operators for the first three Stokes parameters are defined as:

*A.2. Reference frame manipulation*

Upon scattering or crossing of an interface, the reference frame of the Stokes vector must be modified to be expressed in either the scattering plane or the Fresnel plane. For a photon propagating in the *ê*
_{3} direction, the polarization state is described with a Stokes vector **S** whose components are given with respect to a set of arbitrary orthonormal axes *ê*
_{⊥} and *ê*
_{‖}, with *ê*
_{3}=*ê*
_{⊥}×*ê*
_{‖}[5, 22, 32]. To rotate the reference frame of the Stokes vector by an angle *ϕ* around *ê*
_{3} without changing the polarization state, one must transform both the reference frame and the Stokes vector components according to:

The rotation matrix *R*
_{3}(*ϕ*) that rotates a real vector in the basis {*ê*
_{⊥},*ê*
_{‖},*ê*
_{3}} by an angle *ϕ* around *ê*
_{3} is :

the matrix *R*_{S}
(*ϕ*) that transforms the Stokes vector coordinates is:

and, because all vectors are ultimately expressed in Cartesian coordinates {*x̂*,*ŷ*,*ẑ*} in the lab frame, the basis change matrix *B* is:

These combined operations do not change the orientation of the polarization: they merely describe the same polarization state in a different basis. When a scattering plane is selected, this determines the reference frame in which the Stokes vector of the photon must be described before it is transformed.

*A.3. Scattering*

Upon scattering, the direction and polarization state of the light are modified. When a photon scatters by an angle *θ* in the plane normal to *ê*
_{⊥}, its polarization state and its reference frame are transformed according to:

The rotation matrix *R*
_{⊥}(*θ*) is given by:

the basis change *B* is given by Eq. (17), and *M*_{S}
(*θ*) is the Mueller matrix of the scattering event, as obtained in the present work from Mie theory[21] although any formalism may be used. The probability density distribution *P*(*θ*,*ϕ*) of scattering events in a small solid angle around *θ* and *ϕ* is *P*(*θ*,*ϕ*)≠**I**
^{†}·**S′** sin *θ* (note the sin*θ* factor). Various techniques to sample the scattering angles have been presented by other authors. The best two methods are from Jaillon[16] and Kaplan[18], give correct results and are the fastest since they can be implemented with lookup tables that do not depend on the polarization of the incident photon.

*A.4. Interface crossing*

The sample is delimited by a collection of surface elements. The intersection of a photon path with any surface element can be found with vector algebra. The parametric equation for a photon path $\overline{\mathit{SD}}$ starting at a point *S* and ending at *D* (with $\mathbf{d}\equiv \overline{\mathit{SD}}$) is:

where *C* is any point on the line $\overline{\mathit{SD}}$ if *t*≥0 and *t*≤1. The parametric equation for an infinite plane is:

where the point *O* is on the plane and *n̂* is the normal to the plane and *C* is any point on the plane. Using linear algebra, one can show that the intersect between the photon path and the plane is obtained with :

When 0≤*t*≤1, then $\overline{\mathit{SD}}$ crosses the plane and the point *S*+*t*
**d** is on the plane. If the vectors **a** and **b** correspond to the sides of the triangle, the intersection point has coordinates in the basis (*â*,*b̂*) given by :

It is a simple matter to check if the point (*u*,*v*) is on the triangle or not. When the photon propagates inside an object, one must compute the distance to all of the surface elements using Eq. (24) and keep the surface element which has the closest intersection point.

*A.5. Fresnel reflection and transmission*

Upon reaching an interface, the direction of propagation is modified according to Snell’s law and the Stokes vector of the photon is transformed using Fresnel Mueller matrices. The reference frame of the Stokes vector is first rotated such that *ê*
_{⊥} and *ê*
_{‖} are parallel to Fresnel *ŝ* and *p̂* vectors using Eqns (12–14), with *ŝ*=*ê*
_{3}×*n̂′* and *ŝ*×*p̂*=*ê*
_{3}. The surface normal *n̂′* used in Fresnel calculations always points towards the medium that the photon is entering. The photon is reflected if a random number between zero and one is smaller than (${r}_{p}^{2}$
+${r}_{s}^{2}$
)/2+(${r}_{p}^{2}$
-${r}_{s}^{2}$
)**Q**
^{†}
**S**/2**I**
^{†}
**S**, the normalized intensity of the reflected beam for that polarization state (see below)[16]. If the positive angle of incidence is *θ*_{i}
, the reference frame and the polarization state are transformed upon transmission according to:

with *θ*_{t}
the angle of refraction (transmission) with the normal *n̂′*, and the Mueller matrix for Fresnel transmission:

with *t*_{p}
and *t*_{s}
the Fresnel transmission coefficients for the fields. If the relative index of refraction is *m* (index of the destination medium divided by that of the source medium, *i*.*e*., air to glass is 1.5), *t*_{p}
=2cos*θ*_{i}
(cos*θ*_{t}
+mcos*θ*_{i}
)-1 and *t*_{s}
=2cos*θ*_{i}
(cos*θ*_{i}
+*m*cos*θ*_{t}
)^{-1}. In reflection, the polarization state and the reference frame are transformed according to:

with *𝓡*(*θ*_{i}
) the Mueller matrix for Fresnel reflection and is similar to *𝓣*(*θ*_{i}
) but with reflection coefficients *r*_{p}
and *r*_{s}
instead of transmission *t*_{p}
and *t*_{s}
(*r*_{p}
=*m*
*t*_{p}
-1 and *r*_{s}
=*t*_{s}
-1).

*A.6. Detection*

The Stokes vector of the detected photon is calculated in a way that simulates the measurement process in the laboratory. The intensities *I*
_{‖}, *I*
_{⊥}, *I*
_{45}, *I*
_{-45}, *I*
_{R} and *I*
_{L} are calculated explicitly and the Stokes vector is reconstructed. For instance, to calculate the intensity *I*
_{‖} through a polarizer with normal *n̂* that lets through photons polarized along *d̂*
_{‖}, the following is done (see Fig. 7): using Eqns (12), (13) and (14), the reference frame of the Stokes vector is rotated by an angle *ϕ*
_{‖} such that *ê′*
_{‖} is in the plane spanned by *n̂* and *d̂*
_{‖} (which is obtained by looking for a vector perpendicular to the plan spanned by *d̂*
_{⊥} and *ê*
_{3}), then the transmitted intensity is obtained by *I*‖=(**I**
^{†}
**S′**+**Q**
^{†}
**S′**)/2. The process is repeated with axes *d̂*
_{⊥}, *d̂*
_{45}, *d̂*
_{-45}, for intensities *I*
_{⊥}, *I*
_{45} and *I*
_{-45} respectively. For a right (left) circular analyzer, the reference frame is rotated such that *ê′*
_{‖} is in the plane spanned by *n̂* and *d̂*
_{‖}, a *π*/2 (-*π*/2) waveplate transformation is then applied, and finally the detection through linear polarizer along +45° gives *I*_{R}
(*I*_{L}
).

## References and links

**1. **J. M. Schmitt, A. H. Gandjbakhche, and R. F. Bonner, “Use of polarized light to discriminate short-path photons in a multiply scattering medium,” Appl. Opt. **31**, 6535–6546 (1992). [CrossRef] [PubMed]

**2. **S. L. Jacques, J. C. Ramella-Roman, and K. Lee, “Imaging skin pathology with polarized light,” J. Biomed. Opt. **7**, 329–340 (2002). [CrossRef] [PubMed]

**3. **X. Wang, G. Yao, and L. V. Wang, “Monte Carlo model and single-scattering approximation of the propagation of polarized light in turbid media containing glucose,” Appl. Opt. **41**, 792–801 (2002). [CrossRef] [PubMed]

**4. **I. A. Vitkin and E. Hoskinson, “Polarization studies in multiply scattering chiral media,” Opt. Eng. **39**, 353–362 (2000). [CrossRef]

**5. **D. Côté and I. A. Vitkin, “Balanced detection for low-noise precision polarimetric measurements of optically active, multiply scattering tissue phantoms,” J. Biomed. Opt. **9**, 213–220 (2004). [CrossRef] [PubMed]

**6. **R. J. McNichols and G. L. Coté, “Optical glucose sensing in biological fluids: an overview,” J. Biomed. Opt. **5**, 5–16 (2000). [CrossRef] [PubMed]

**7. **D. R. Lide, (ed.), *CRC Handbook of Chemistry and Physics* (CRC Press LLC, Boca Raton, Florida, 1998), pp. 3–12,8–64., 79th edn.

**8. **W. F. March, B. Rabinovitch, and R. L. Adams, “Noninvasive glucose monitoring of the aqueous humor of the eye: Part II. Animal studies and the scleral lens,” Diabetes Care **5**, 259–265 (1982). [CrossRef] [PubMed]

**9. **B. D. Cameron and G. L. Coté, “Noninvasive glucose sensing utilizing a digital closed-loop polarimetric approach,” IEEE Trans. Biomed. Eng. **44**, 1221–1227 (1997). [CrossRef] [PubMed]

**10. **R. R. Ansari, S. Bockle, and L. Rovati, “New optical scheme for a polarimetric-based glucose sensor,” J. Biomed. Opt. **9**, 103–115 (2004). [CrossRef] [PubMed]

**11. **A. J. Welch, G. Yoon, and M. J. van Gemert, “Practical models for light distribution in laser-irradiated tissue,” Lasers Surg Med **6**, 488–493 (1987). [CrossRef] [PubMed]

**12. **M. S. Patterson, B. C. Wilson, and D. R. Wyman, “The propagation of optical radiation in tissue I. Models of radiation transport and their application,” Lasers in Medical Science **6**, 155–166 (1990). [CrossRef]

**13. **L. Wang, S. L. Jacques, and L. Zheng, “MCML-Monte Carlo modeling of light transport in multi-layered tissues,” Comput. Meth. Prog. Biomed. **47**, 131–146 (1995). [CrossRef]

**14. **M. Moscoso, J. B. Keller, and G. Papanicolaou, “Depolarization and blurring of optical images by biological tissue,” J. Opt. Soc. Am. A **18**, 948–960 (2001). [CrossRef]

**15. **J. R. Mourant, T. M. Johnson, and J. P. Freyer, “Characterizing mammalian cells and cell phantoms by polarized backscattering fiber-optic measurements,” Appl. Opt. **40**, 5114–5123 (2001). [CrossRef]

**16. **F. Jaillon and H. Saint-Jalmes, “Description and time reduction of a Monte Carlo code to simulate propagation of polarized light through scattering media,” Appl. Opt. **42**, 3290–3296 (2003). [CrossRef] [PubMed]

**17. **S. Bartel and A. H. Hielscher, “Monte Carlo simulations of the diffuse backscattering Mueller matrix for highly scattering media,” Appl. Opt. **39**, 1580–1588 (2000). [CrossRef]

**18. **B. Kaplan, G. Ledanois, and B. Drévillon, “Mueller Matrix of dense polystyrene latex sphere supsensions: measurements and Monte Carlo simulations,” Appl. Opt. **40**, 2769–2777 (2001). [CrossRef]

**19. **M. Mehrübeoglu, N. Kehtarnavaz, S. Rastegar, and L. V. Wang, “Effect of molecular concentrations in tissue-simulating phantoms on images obtained using diffuse reflectance polarimetry,” Opt. Express **3**, 286–297 (1998), http://www.opticsexpress.org/abstract.cfm?URI=OPEX-3-7-286. [CrossRef] [PubMed]

**20. **D. Côté and I. A. Vitkin, “Pol-MC: a three-dimensional polarization-sensitive Monte Carlo implementation for light propagation in tissue,” Available at http://www.novajo.ca/ont-canc-inst-biophotonics/.

**21. **T. A. Germer, “SCATMECH: Polarized Light Scattering C++ Class Library,” Available at http://physics.nist.gov/scatmech.

**22. **H. C. van de Hulst, *Light scattering by small particles* (Dover, New York, 1981).

**23. **W. J. Wiscombe, “Improved Mie scattering algorithms,” Appl. Opt. **19**, 1505–1509 (1980). [CrossRef] [PubMed]

**24. **J. S. Maier, S. A. Walker, S. Fantini, M. A. Franceschini, and E. Gratton, “Possible correlation between blood glucose concentration and the reduced scattering coefficient of tissues in the near infrared,” Opt. Lett. **19**, 2062–2064 (1994). [CrossRef] [PubMed]

**25. **J. M. Steinke and A. P. Shepherd, “Diffusion model of the optical absorbance of whole blood,” J. Opt. Soc. Am. B **5**, 813–822 (1988). [CrossRef]

**26. **V. Sankaran, J. T. Walsh, and D. J. Maitland, “Polarized light propagation through tissue phantoms containing densely packed scatterers,” Opt. Lett. **25**, 239–241 (2000). [CrossRef]

**27. **A. N. Yaroslavsky, I. V. Yaroslavsky, T. Goldbach, and H.-J. Schwarsmaier, “Influence of the scattering phase function approximation on the optical properties of blood determined from the integrating sphere measurements,” J. Biomed. Opt. **4**, 47–53 (1999). [CrossRef]

**28. **N. Ghosh, P. K. Gupta, H. S. Patel, B. Jain, and B. N. Singh, “Depolarization of light in tissue phantoms -effect of collection geometry,” Opt. Comm. **222**, 93–100 (2003). [CrossRef]

**29. **V. Sankaran, K. Schönenberger, J. T. Walsh, and D. J. Maitland, “Polarization discrimination of coherently propagating light in turbid media,” Appl. Opt. **38**, 4252–4261 (1999). [CrossRef]

**30. **M. P. Silverman, W. Strange, J. Badoz, and I. A. Vitkin, “Enhanced optical rotation and diminished depolarization in diffusive scattering from a chiral liquid,” Opt. Comm. **132**, 410–416 (1996). [CrossRef]

**31. **K. C. Hadley and I. A. Vitkin, “Optical rotation and linear and circular depolarization rates in diffusively scattered light from chiral, racemic, and achiral turbid media,” J. Biomed. Opt. **7**, 291–299 (2002). [CrossRef] [PubMed]

**32. **C. F. Bohren and D. R. Huffman, *Absorption and scattering of light by small particles* (Wiley, New York, 1983), *chap. 2*, pp. 46–56.