## Abstract

Photon pathlength distributions as a function of the number of scattering events in cylindrical turbid samples are studied using a polarization-sensitive Monte Carlo model with linearly polarized light input. Sample scattering causes extensive depolarization, yielding a photon field comprised of polarized and depolarized sub-populations. It is found that the pathlength of polarization-preserving photons is distributed within a defined spatial range with strong angular dependence. This pathlength, averaged over the range, is 2-3X smaller than the one averaged over the widely-spread range of all (polarized + depolarized) collected photons. It is also demonstrated that changes in optical properties of the media affect the pathlength distributions.

© 2007 Optical Society of America

## 1. Introduction

When light is used for medical diagnostics, it is crucial to ensure that it reaches the intended target within tissue, and then carries the information to the detector. Because biological tissue is highly scattering, the photons transport occurs in complicated zigzag fashion. This makes the unambiguous determination of penetration depth and sampling volume of injected photons difficult. A detailed description of pathlength distribution may help in estimating these important parameters.

The pathlength distribution within multiply scattering media has been of great interest [1–4]. There are two main theoretical methods to address this problem. One is an analytical approach through the diffusion equation [5–7]. Another is a numerical approach through Monte Carlo simulations (MC) [8–11]. The analytical approach is fast, but the diffusion results may not be accurate near light sources, or in regions of high optical gradients (near surfaces and boundaries) [12]. Conversely, Monte Carlo simulations provide a flexible and accurate, if somewhat slower, platform that traces a large number of photon histories in complex turbid media, thus yielding the pathlength distributions directly.

So far, Monte Carlo simulations for pathlength distribution studies have mostly examined slab and semi-infinite geometries. In the context of noninvasive blood glucose monitoring, for example, the cylindrical geometry is of special relevance because the measuring sites are often of quasi-round curvature like finger tips and lips. Compared with slab and semi-infinite geometries, the cylindrical geometry may actually be advantageous. Specifically, it permits multiple-direction detections (0° to 360°, compared with 0° and 180° detections for slab and 180°-only detection for semi-infinite set-ups). This flexibility may enable a more complete study of the scattering properties of turbid media. This multiple-direction detection capability may also be more practical and convenient in the clinical environment. It is thus useful to expand the studied geometries further to include cylindrical samples.

Furthermore, polarized light has been found useful in studying highly turbid media, because the properties of the media can be extracted through changes in the photon’s polarization state, such as optical rotation α and surviving linear polarization fraction β_{L} [13–21]. Therefore, in this paper, we investigate the pathlength distribution of polarized and depolarized photon populations in highly scattering media in a cylinder irradiated with polarized incident beam, by means of the polarization-sensitive Monte Carlo simulations developed in our group [13,14]. The effects of sample geometry and sample optical property changes on the pathlength distribution are studied in detail.

## 2. Theory: Monte Carlo simulations in cylindrical geometry

In the polarization-sensitive Monte Carlo simulation, the photons are propagated between scattering events, as determined by a pseudo-random sampling of the scattering mean free path. In our validated approach [13,14], which expands upon the earlier models [15,16], the polarization information is tracked in the form of individual Stokes vectors [14]. These are summed over a large number of tracked photon histories, to yield the experimentally observable macroscopic polarization and intensity properties of interest.

Specifically, the position, propagation direction and polarization (described by Stokes vector formulism) of a photon are initialized at the entrance of a cylindrical sample (40 mm height and 8 mm diameter in this paper) which is characterized by a set of surface elements: rectangular on the sides and triangular on the bottom and top (48 on the sides, 48 on the top and 48 on the bottom in this study). The curvature of the side element can be ignored. So the rectangular element is treated as a plane (0.52 mm × 40 mm). The photon is moved by a scattering event within the sample if it does not cross an interface. The Stokes vector of the photon is transformed by the application of the Mueller matrix of the scattering event, which is calculated from Mie theory [22]. If the photon encounters an interface, its Stokes vector is transformed into the Fresnel reference frame. The photon’s Stokes vector and propagation direction are modified according to reflection or transmission. If the photon is reflected, the propagation inside the sample continues. If the photon is transmitted and outside the sample, the measured Stokes vector is computed from the sum and differences of intensities through linear and circular polarizers. The macroscopic values of the Stokes coefficients are obtained from the sum of the Stokes vector components of incoming photons, after transformation through polarizers.

The transmitted photons are binned spatially into the surface elements, which is divided into smaller surface detection element, based on where they exit the sample and binned angularly within each surface detection element based on exit angle that the propagation vector of the photon makes with the normal to the interface, which is defined as acceptance angle ψ (from 0° to 90°). The size of surface detection element and acceptance angle affects the detected photon intensity, degree of polarization and pathlength distribution. The extent of the effects varies with detection geometry. In this paper, we chose the sizes of surface detection element (~ 0.6 mm^{2}) and acceptance angle (48°) similar to our experimental system. The detected photons can be further binned based on the number of scattering events they undergo within the sample, which is the key investigation method employed in this study.

The cylindrical geometry, as shown in Fig. 1, mimics the experimental conditions used in our research [21]. A 632.8 nm horizontally polarized beam of 1 mm diameter is incident at O on the center of a vertically orientated cylindrical sample of 8 mm diameter and 40 mm height. The scattered photons at point P(z, θ), within acceptance angle ψ(48°), are summed over a 0.6 mm^{2} detection area on the surface of the cylinder. The detection angle θ, which is the angle between forward direction and the normal of the surface detection element, varies from 0° to 180°. The vertical position of the surface detection element z ranges from -20 mm to +20 mm, with the signs indicating the relative position with respect to the horizontal incident plane. The samples simulate water suspensions of 4.1 μm diameter polystyrene microspheres. For most simulations (except Fig. 8), the scattering coefficient of the medium μ_{s} is set to 100 cm^{-1}. The scattering coefficient range is chosen to approximate typical turbidity of biological tissue. The refractive index of the scattering particles is 1.59. Calculated from Mie theory [22], the scattering anisotropy g (a measure of the amount of forward direction retained after a single scattering event, described by mean value of cos γ, where γ is scattering angle) is 0.88 and the scattering efficiency Q_{sca} is 2.71. The absorption coefficient is set to 0.00326 cm^{-1} for water. To reduce the statistical uncertainty of the values obtained from the Monte Carlo simulations, a large number (10^{9}) of horizontally polarized photons are launched in the simulations.

In highly turbid media, only a small fraction of input photons preserve their polarization. The study of this photon subpopulation is of great importance, since they could reveal some unique information not available through intensity only (depolarized) signal channels (e.g., optical rotation of the incident polarization vector, relative depolarization rates, etc.). However the properties of the surviving polarized photons can be easily masked by vast depolarized population. In order to get a clearer picture of characteristics of the surviving polarized photons such as their intensity and pathlength distributions as functions of number of scattering events, the collected photons are binned based on the number of scattering events they experienced within the sample and compared with the total (N-unresolved) averaged ones. The bin number, referring to the number of scattering events, runs from N=1 to N~70. In our preliminary simulations, after 70 multiple scattering events, the surviving linear polarization fraction dropped to the noise level, that is, no polarized photons were observable for N≥70. The N-unresolved cumulative bin is referred to as Total, containing all the collected photons (N-unresolved, cumulative) without discrimination, that is, N from 1 to infinity. The parameters obtained from the simulations are thus indexed ones (I_{N}, β_{LN} and L_{N}) and cumulative ones (I_{total}, β_{Ltotal}, L_{total} and N_{total}). I_{N} is the summed number of photons in bin N, including both polarized and depolarized ones. β_{LN} is the surviving linear polarization fraction of those photons, and L_{N} is their average pathlength. I_{total}, β_{Ltotal} and L_{total} are the corresponding parameters for the bin Total. N_{total} is the average number of scattering events of the total collected photons,

## 3. Results and discussion

#### 3.1 Photon intensity distributions

Figure 2 shows the calculated photons intensity distribution as functions of number of scattering events N and detection angle θ in the incident plane (z=0). The data from the forward hemisphere are not shown due to large statistical uncertainty resulting from inadequate photon number, as the detected intensity is mostly in the backwards hemisphere. Fig. 2(a) presents the indexed intensity I_{N} distribution. The photons (comprised of both polarized and depolarized sub-populations) are not evenly distributed in N, and the shape of the distribution is detection-direction dependent. In backward direction, θ=180°, the photon intensity peaks at N=1 (note the logarithmic scale on vertical axis), and decreases rapidly with N value. It indicates that single-scattering dominates the detected signal in the exact back-scattered geometry. For smaller angles, the intensity distributions are not monotonically decreasing, but instead show broad peaks. Both the minimum number of scattering events needed for photons to escape from a certain geometry and the most likely number of scattering events contributing to the detected signal can be estimated from the presented data. Taking θ=151° curve for example, photons have to be scattered at least 8 times before they can exit, while the photons which undergo ~25 scattering events have more chances to re-emit in this geometry than photons with other scattering histories.

The intensity distribution of total collected photons I_{total} as a function of detection angle is plotted in Fig. 2(b). It shows general increase of photon intensity with detection angle (i.e., detected intensity is higher in the backwards hemisphere), but lacks detailed description of photon scattering interactions, represented by I_{N}.

To examine in detail the polarized and the depolarized subpopulations, plotted in Fig. 3(a) are the indexed surviving linear polarization fractions as functions of N and θ. It is seen that the surviving linear polarization fraction generally decreases with N, which is understandable because scattering randomizes polarization. The simulation results show that the total degree of polarization DOP, defined as (Q^{2}+U^{2}+V^{2})^{1/2}/I (I, Q, U and V are elements of the Stokes vector of detected photons) [22], has similar value to β_{L}, for example, for θ=180° DOP_{N=1}=0.22492493, β_{L1}=0.22491533 and for θ=135° DOP_{N=22}=0.55322760, β_{L22}=0.55322686. It implies that the low β_{LN} value (< 100%) is the result of true depolarization, and not due to transformation of the polarization form (linearly polarized to elliptically/circularly polarized light). However, the low β_{LN} value at 180°, N=1 (~ 23%) is surprising in that single scattering by identical particles with spherical symmetry should not decrease the degree of polarization of 100% polarized incident light, although the nature of polarization may be changed [2,22]. It is also found that the β_{L1} value for θ=180° is dependent on sample (anisotropy factor g and scatterer radius r) and detection system (acceptance angle ψ). Different g values were obtained by varying r while keeping the scattering coefficient constant in the simulations. Note that g is not unique with scatterer size r. Therefore, samples with the same g may be simulated differently, as our Monte Carlo simulations use the full calculated Mie scattering phase function. For example, for the same g=0.88 as in our current simulations, but with smaller scatterer size [r=0.33 μm, versus 2.05 μm in Fig. 2(a)], β_{L1} value reaches ~65% for ψ=48° [in contrast to ~23% in Fig. 3(a)]. For θ=180° and ψ approaching zero, β_{L1} could be as high as ~94% for g=0.74 (r=0.2 μm). So the degree of polarization after a single scattering interaction appears to vary widely in magnitude, and to depend on medium properties and detection geometry. Further studies with circularly polarized light input and birefringent media properties are under way to check if these simulation results are seen under varying conditions, and to help interpret this phenomenon.

It is found from the simulations that same amount of scattering events may cause greater depolarization at higher angles than at smaller angles, as seen in Fig 3(b). Taking N=34 curve for example, only 4.7% of the photons escaping from θ=165° preserves their polarization after being scattered 34 times, in contrast with 46% when escaping from θ=121° following a similar scattering history. Consequently, both the number of scattering events and detection direction (geometry) determine resultant depolarization. 44 scattering events are large enough to effectively randomize the polarization of the photons at θ≥158°, but the light field is still ~ 40% polarized at this same N values detected at θ~100°. For θ>165°, the angular dependence of depolarization is weakened and somewhat reversed (slight change in β_{L}).

It is also important to know the fraction of the indexed surviving polarized photon subpopulation (β_{L}I_{N}) in the total collected photon populations I_{total}. Presented in Fig. 4(a) are plots of the ratio R_{N} of indexed surviving polarized photon intensity to the total collected photon intensity as a function of number of scattering events. R_{N} is calculated from

It is seen that the polarized photons are distributed in defined N ranges, with peak position decreasing with θ, while the magnitude of the peak generally increasing with θ (note that there are no lower-N bounds for θ=180° and 173° curves). For example, at θ=135°, the polarized photons are distributed between N=17 and 52, peaking at N~32. The numerical integration of the ratio over the non-zero N range, ∑^{52}
_{N=17}R_{N}(135°,z=0), is equal to the cumulative surviving linear polarization fraction β_{Ltotal} [7.5% for this case, see numerical labels on Fig. 3(a)], which means all the surviving polarized photons detected at that angle are contained in bins N=17 to N=52. For θ=180°, polarized photons are distributed between N=1 and N=20. The peak intensity (at N=1) is about four times higher than the one at θ=135° (N=32), implying the potential advantage of backwards detection geometries for enhanced polarization preservation.

Illustrated in Fig. 4(b) are the cumulative surviving linear polarization fractions β_{Ltotal} at different detection angles. The differential β_{LN} and cumulative β_{Ltotal} show opposite angular trends: β_{LN} was shown [Fig. 3(b)] to decrease with detection angle θ, while β_{Ltotal} increases with θ. The increase of β_{Ltotal} with detection angle in Fig. 4(b) is the result of the larger portion of depolarized photons of high N values at smaller angles. For example, at θ=135° the total average number of scattering events N_{total}=102, compared with N_{total}=17 at θ=180°.

In Fig. 4(b), the data point at θ=165° (g=0.88) appears to deviate from the main trend. It has been found that the behavior at θ=165° is related to the anisotropy factor g as shown in the inset: when g decreases from 0.93 to 0.5, the “anomalous” θ=165° behavior disappears. Further studies are currently addressing this phenomenon.

#### 3.2 Photon pathlength distributions

Figure 5(a) illustrates the indexed average pathlength distribution of photons (both polarized and depolarized). For θ=180° and 173°, average pathlength increases linearly with N, with a slope of 0.0108 cm per scattering event, close to the mean free path (mfp), calculated from 1/μ_{s} (0.01 cm). However, the average pathlengths for smaller detection angles deviate from the main trend, especially at smaller N. For instance, the 121° curve starts converging at N~45, whereas the 158° merges at N~22. Sample geometry effects can explain this, as shown in Fig. 5(b). The flat parts of the curves in Fig. 5(b) correspond to the linear regions in Fig.5(a). The reference line is the direct distance between entrance and exit of the sample (OP in the top view of the sample cylinder) at all detection angles [see inset in Fig. 5(b)]. This reference line sets lower pathlength limit, in that only photons that are able to travel at least that pathlength can escape in a given direction. For example, N~20 scattering events are sufficient for photons to exit at θ=180°, 173° and 165° [flat region of N=20 curve in Fig. 5(b)], but are insufficient for the photons to escape at θ=135°. Only those that travel a longer pathlength after being scattered N~20 times can exit at that angle. However, these photons only account for a small fraction of total photon population at θ=135° [see 135° curve in Fig. 2(a)]. This geometry effect on pathlength is also seen when simulations are performed off the incident plane (z≠0, results not shown).

Illustrated in Fig. 6 is the indexed average pathlength divided by the number of scattering events (L_{N}/N), which we shall term “the unit pathlength”. The ratio is not generally equal to the calculated mean free path 1/μ_{s}=0.01 cm (see the reference line), but asymptotically approaches it at large N values. This behavior can be interpreted as follows. Mean free path of photons is a statistical estimation of the distance which photons travel between successive scattering events, equal to L_{N}/(N+1). It can be approximated to L_{N}/N when N is large enough so that the “1” in the denominator can be neglected, which can help explain some deviation at low (N≤10) values. With the increase of N, the effect is diminishing and L_{N}/N gradually approaches 1/μ_{s}. For the smaller detection angles (θ<173°), another more important mechanism is involved: the geometry induced longer pathlength as discussed above [see Fig. 5(b)]. It makes the difference between L_{N}/N and 1/μ_{s} larger than that at the high θ values. The smaller θ value curves approach 180° and 173° results when the geometry effect decreases at high N values.

So far, the pathlength discussion has dealt with the average pathlength distribution of all photons in bins without polarization differentiation. The pathlength distribution of the polarized photon fraction is also of interest. The MC simulations with both polarized and unpolarized light inputs demonstrate that the indexed pathlength L_{N} is polarization-independent. For a given number of scattering events N, polarized and unpolarized photons, escaping from same detection direction traverse similar pathlengths (provided that N is small enough such that the polarized subpopulation still exits). Therefore, the pathlength distribution as a function of number of scattering events, discussed above in Fig. 5(a), can be used to represent pathlength distribution of polarized photons. Combining Fig. 4(a) with Fig. 5(a), we get distribution of polarized photons as a function of pathlength, shown in Fig. 7(a). It reveals that the polarized photons have a defined pathlength range for a given detection geometry. For example, at θ=135°, the pathlength of the polarized subpopulation ranges from ~ 0.33 cm to ~ 0.57 cm. Here, the photon population threshold for the upper and lower pathlength limits is set to 10% of the peak value. The range can be made smaller by raising the threshold. We define the average of this pathlength range as the characteristic polarization pathlength, L_{CP}. That is,

where i and j are minimum and maximum number of scattering events the surviving polarized photons experienced for a given detection geometry. For example, i=8, j=40 for θ=158° [see Fig. 4(a)]. It better characterizes the pathlength of polarized photons than L_{total},

whose calculation includes longer pathlengths at N>j values where no polarized photons are present. This is illustrated in Fig. 7(b), which plots the described characteristic polarization pathlength L_{CP} (using 10% of peak value threshold) as a function of θ, and compares it with the total average pathlength L_{total}. Both curves show similar trends, with pathlengths decreasing with detection angle θ. However, L_{total} is 2-3X longer than L_{CP}. If pathlength is used to estimate the penetration depth of injected photons, L_{CP} will be more appropriate than L_{total} to describe the interaction extent of polarization-preserving subpopulation.

#### 3.3 Effects of optical properties on pathlength

The Monte Carlo simulations demonstrate that the change in optical properties of the sample affects the pathlength distribution, as shown in Fig. 8. For a fixed number of scattering events (N=30 in this example), decreasing the scattering coefficient μ_{s} results in longer pathlength because of increasing the mean free path. The anisotropy factor g doesn’t have as much influence on the pathlength as μ_{s}, especially at backward directions (θ=180° and 173°). When the lumped optical property reduced scattering coefficient μ_{s}’, which incorporates scattering coefficient μ_{s} and anisotropy g by μ_{s}’=μ_{s}(1-g), increases from 6 cm^{-1} to 25 cm^{-1}, the pathlength varies in a non-monotonic way. In general, for any optical property set, the pathlength (and hence the sampling volume) decreases as the detection direction approaches the exact backward geometry (θ=180°).

## 4. Conclusion

Utilizing a polarization-sensitive Monte Carlo simulation model, we studied the pathlength distribution of photons as a function of number of scattering events in cylindrical turbid samples (μ_{s}~100 cm^{-1}). The pathlength of the escaping photons spreads out due to multiple scattering. However, the pathlength of the surviving polarized photons is distributed in a more defined range, the average over which is 2-3X smaller than the average pathlength of all collected photons, which is dominated by the longer pathlengths of depolarized photon populations. Therefore, the average over this range, defined as characteristic polarization pathlength, is more descriptive as the pathlength of the polarization-preserving photon subpopulation. The strong angular dependence of the characteristic polarization pathlength suggests that penetration depth of injected photons, as estimated by pathlength distribution, can be controlled by adjusting the detection angle. The simulation results demonstrate that the change in scattering coefficient μ_{s} and anisotropy factor g of turbid media impacts the photon pathlength distribution.

## Acknowledgments

The authors thank the Natural Sciences and Engineering Research Council of Canada (NSERC) for financial support of this research.

## References and links

**1. **G. H. Weiss, R. Nossal, and R. F. Bonner, “Statistics of penetration depth of photons re-emitted from irradiated tissue,” J. Mod. Opt. **36**,349–359 (1989). [CrossRef]

**2. **M. Dogariu and T. Asakura, “Photon pathlength distribution from polarized backscattering in random media,” Opt. Eng. **35**,2234–2239 (1996). [CrossRef]

**3. **Y. Tsuchiya, “Photon path distribution and optical responses of turbid media: theoretical analysis based on the microscopic Beer-Lambert law,” Phys. Med. Biol. **46**,2067–2084 (2001). [CrossRef]

**4. **Y. Liu, Y. L. Kim, X. Li, and V. Backman, “Investigation of depth selectivity of polarization gating for tissue characterization,” Opt. Express **13**,601–611 (2005). [CrossRef]

**5. **M. S. Patterson, B. Chance, and B. C. Wilson, “Time resolved reflectance and transmittance for the non-invasive measurement of tissue optical properties,” App. Opt. **28**,2331–2336 (1989). [CrossRef]

**6. **S. R. Arridge, M. Cope, and D. T. Delpy, “The theoretical basis for the determination of optical pathlengths in tissue: temporal and frequency analysis,” Phys. Med. Biol. **37**,1531–1560 (1992). [CrossRef]

**7. **G. Popescu and A. Dogariu, “Dynamic light scattering in subdiffusive regimes,” Appl. Opt. **40**,4215–4221 (2001). [CrossRef]

**8. **B. C. Wilson and G. Adam, “A Monte Carlo model for the absorption and flux distributions of light in tissue,” Med. Phys. **10**,824-830 (1983). [CrossRef]

**9. **S. T. Flock, M. S. Patterson, B. C. Wilson, and D. R. Wyman, “Monte Carlo modeling of light propagation in highly scattering tissues-I: model predictions and comparison with diffusion theory,” IEEE Trans. Biomed. Eng. **36**,1162–1168 (1989). [CrossRef]

**10. **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]

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

**12. **L. Wang and S. L. Jacques, “Hybrid model of Monte Carlo simulation and diffusion theory for light reflectance by turbid media,” J. Opt. Soc. Am. A **10**,1746–1752 (1993). [CrossRef]

**13. **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]

**14. **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**,148–163 (2005). [CrossRef]

**15. **B. Kaplan, G. Ledanois, and B. Drévillon, “Muller Matrix of dense polystyrene latex sphere suspensions: measurements and Monte Carlo simulation,” Appl. Opt. **40**,2769–2777 (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]

**17. **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]

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

**19. **M. Mehrübeoǧlu, 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). [CrossRef]

**20. **G. L. Coté, M. D. Fox, and R. B. Northrop, “Noninvasive optical polarimetric glucose sensing using a true phase measurement technique,” IEEE Trans. Biomed. Eng. **39**,752–756 (1992). [CrossRef]

**21. **X. Guo, M. F. G. Wood, and I. A. Vitkin, “Angular measurements of light scattered by turbid chiral media using linear Stokes polarimeter,” J. Biomed. Opt. **11**,041105 (2006). [CrossRef]

**22. **C. F. Bohren and D. R. Huffman, *Absorption and Scattering of Light by Small Particles* (Wiley, New York, 1983).