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(-µ_td), 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 (α=γ-γ _∘)

 

Fig. 1. Schematic diagram of the experimental setup used for the measurements. Chop.-mechanical chopper; P₁: polarizer; PEM- photo-elastic modulator; A_1,2: apertures; PC: polarizing beam splitter cube; L_x,y- lenses; D_x,y- photodetectors. f_c and f_p are the modulation frequencies for the mechanical chopper and PEM respectively.

Download Full Size | PDF

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⁹ to as much as 10¹¹ 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 ₁(θ ₁) 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_{2f_p}$, the rotation α can be obtained with:

where Q_fc and $Q_{f_{2f_p}}$ are the measured signals at modulation frequencies f_c and $f_{2f_p}$, ℱ(2f _p)/ℱ(f _c) is the ratio of the frequency responses of the detector at those two frequencies (0.56 in our case), and J ₂(δ _∘) is the second-order Bessel function of the first kind (at δ _∘=3, J ₂(δ _∘) 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³) 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 2r=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 ²) 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 =3Q _scat f_w ρ _∘/4rρ(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⁶. 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

Table 1. Table comparing Pol-MC results with MCML. The MCML simulation is done with three layers (air, tissue, air) of indices of (1.0,1.33,1.0), the Henyey-Greenstein phase function with anisotropy parameter g=0.93, and µ_s =10 cm^-1, µ_a =0.1 cm^-1 in the tissue layer. The Pol-MC calculation is done using two different phase functions: the Henyey-Greenstein phase function with anisotropy parameter g=0.93 (no polarization considered in this case), and the phase function from Mie theory for spherical scatterers of 1.4 µm-diameter (with an anisotropy g of 0.93). In all three cases, there is a 2% specular reflection at the first interface because of the index mismatch between air and tissue.

View Table | View all tables in this article

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

Table 2. Table showing the effect of finite size on diffuse reflectance (back scattering) and transmittance (forward scattering) for calculations done with three layers (air, tissue, air) of indices of (1.0,1.33,1.0) with µ_s =10 cm^-1 and µ_a =0.1 cm^-1 with spherical polystyrene scatterers of 1.4 µm-diameter (which have a Mie scattering phase function with an anisotropy g of 0.93). As the lateral extent of the sample goes from infinite to finite, the transmittance and reflectance go down since a larger fraction of the light is lost through the sides of the sample.

View Table | View all tables in this article

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.

 

Fig. 2. a) Experimental measurements (squares) and calculation (line) of surviving linear polarization fraction for the forward-scattered beam in 1 µm diameter microspheres suspensions of various scattering coefficients, with no absorption (phase function with anisotropy g of 0.917). The error bars (shown on graph) are smaller than the symbols for most measurements. The cuvette has a thickness of 1 cm, the imaged area is approximately 5 mm in diameter with an acceptance angle of 15° for the apertures and lenses used. No adjustable parameters were used. b) Experimental measurements (symbols, no error bars provided) from Fig. 2a and 2c in Ghosh et al.[28] for surviving linear polarization fraction of the forward-scattered beam in cuvettes for thicknesses 0.5 cm (blue line) and 1 cm (red line) with 1.072 µm diameter (g=0.923) microspheres suspensions of various scattering coefficients, and corresponding calculations (lines, same colors as symbols) with the model presented here.

Download Full Size | PDF

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.

 

Fig. 3. Experimental measurements (squares) and Monte Carlo predictions of the optical rotation (solid black line) as a function of scattering coefficient µ_s , in the transmission direction through a turbid 1 cm cuvette for solutions containing 1.2 M of glucose and 1.4 µm diameter microspheres (phase function with anisotropy of 0.943, no adjustable parameters were used in the calculation). The absorption coefficient µ_a is negligible. The rotation as predicted from the average path length if Eq. (1) held is also shown on the graph (blue line), and can be seen to increase with scattering coefficient, unlike the rotation which reaches a plateau.

Download Full Size | PDF

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°.

 

Fig. 4. Monte Carlo calculations of the beam orientation as a function of exit position for the (a) forward, (b) side and (c) back scattering geometry. The scattering coefficient is 20 cm^-1 and the scatterer diameters are 1.0 (blue line), 1.4 (red) and 1.8 (green) µm, absorption is negligible. The scatterers have different phase functions and different anisotropies (0.917, 0.930 and 0.922 respectively) because of the different sizes. The beam is incident with a 45° polarization and all three solutions contain 1 M of glucose (effect of glucose on background index not considered). Regardless of the scatterer size, the rotation on axis in the forward direction is the same. Notice the sensitivity to exit position in directions other than the forward direction, with a sharp orientation change in the back direction. Another solution with no glucose and 1.8 µm scatterer diameter (black line) demonstrates that this change of orientation is solely due to scattering since the sharp change in orientation is also observed in that case.

Download Full Size | PDF

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.

 

Fig. 5. Orientation of the polarization of a 45° incident beam after single scattering from a spherical scatterer (polystyrene spherical scatterer of 1.4 µm diameter in water, g=0.93). The polarization changes from its original 45° in the forward direction to 135° in the back scattering direction. Inset: Diagram illustrating the orientation of polarization in two limiting cases (forward and backward scattering).

Download Full Size | PDF

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 ê ₃ 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 ê ₃=ê _⊥×ê _‖[5, 22, 32]. To rotate the reference frame of the Stokes vector by an angle ϕ around ê ₃ without changing the polarization state, one must transform both the reference frame and the Stokes vector components according to:

The rotation matrix R ₃(ϕ) that rotates a real vector in the basis {ê _⊥,ê _‖,ê ₃} by an angle ϕ around ê ₃ 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.

 

Fig. 6. Schematic diagram of the important points and vectors used in calculating the intersection between the photon path and a triangular surface element.

Download Full Size | PDF

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 ŝ=ê ₃×n̂′ and ŝ×p̂=ê ₃. 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/2I ^† 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 +mcosθ_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).

 

Fig. 7. Polarized photon incident on detecting interface, where d̂ _‖ and d̂ _⊥ represent the orientation of linear polarizers in the laboratory. The reference frame of the Stokes vector is rotated such that ê _‖ is parallel to d̂ _⊥×ê ₃. The intensity that is detected is obtained by I _‖ in that reference frame.

Download Full Size | PDF

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 ₄₅, 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 ê ₃), then the transmitted intensity is obtained by I‖=(I ^† S′+Q ^† S′)/2. The process is repeated with axes d̂ _⊥, d̂ ₄₅, d̂ _-45, for intensities I _⊥, I ₄₅ 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 ).