Three Monte Carlo programs were developed which keep track of the status of polarization of light traveling through mono-disperse solutions of micro-spheres. These programs were described in detail in our previous article . This paper illustrates a series of Monte Carlo simulations that model common experiments of light transmission and reflection of scattering media. Furthermore the codes were expanded to model light propagating through poly-disperse solutions of micro-spheres of different radii distributions.
© 2005 Optical Society of America
Three Monte Carlo programs for polarized light transfer into scattering media have been described and verified in detail . In this paper we present various numerical and experimental tests that were conducted to validate all three programs.
Polarized light transport in mono-disperse solutions of micro-spheres has been extensively studied; the literature is rich with papers dealing on this subject both experimentally and theoretically. When conducting a standard Monte Carlo program simulation some typical tests are executed, such as total transmission through a slab, or total reflection from a semi-infinite media. This is also true for polarized light Monte Carlo programs. In addition to the classic numerical test that may be compared with adding doubling simulations, other polarization specific simulations may be conducted. One such test is the evaluation of the Mueller matrix backscattered from a solution of micro-spheres. Backscattered photons will exit the solution and form 16 unique intensity patterns. Bartel et al.  used this test to validate their Monte Carlo program with experimental results, as well as Cameron et al. , and Rakovic et al. . Several authors are now using the deviation from this standard pattern to study solution chirality and composition [5, 6]. In almost all the previous implementations, as well as ours, the polarization is represented in terms of the Stokes vector. For this reason the Monte Carlo programs are unable to simulate coherent phenomena and should be used only in cases where information on the phase is not necessary. For coherent phenomena the code by Min Xu  could be used.
In this paper we present experimental and numerical results for the backscattered Mueller matrix, with the addition of skewed illumination. This simulation is important for breaking the experiment geometrical symmetry, and to facilitate understanding of the numerical results.
Polarized light transfer through a slab of scattering solutions is also used as a test for the validity of Monte Carlo programs. We conducted experiments in transmission with mono-disperse solutions of micro-spheres that agree with the corresponding Monte Carlo simulations moreover we extended the three Monte Carlo programs to handle light scattering through poly-disperse solutions of micro-spheres. The Monte Carlo results were compared to Adding Doubling simulations , and the particle size distribution was modeled with a modified Gamma distribution. The results agree with the Monte Carlo simulations within the statistical error determined by the number of photons launched (N), where error = 1/√N. The effect of the average anisotropy of a poly-disperse solution on a polarized beam traveling through the media was also studied.
Unless otherwise specified, the number of photons launched in all simulations was 106. The simulations were conducted on a 1.7 GHz Pentium 4 computer, using the Linux Red Hat 7.1 operating system.
2. Mueller Matrix of backscattered light from micro-spheres solutions
Monte Carlo simulations with polarization are often presented in Mueller Matrix form. Bartel et al.  and Cameron et al. , showed results of runs conducted with mono-disperse solutions of micro-spheres of diameter 2.02 μm, index of refraction equal to 1.59, and 1.33 for the index of refraction of water.
Cameron’s experimental set-up is shown in Fig. 1. The incident laser beam (λ= 0.543 μm) was normal to the sample, and the incident beam was polarized by polarizing optics, such as a linear polarizer and quarter-wave plate. A front surface mirror tilted at a ~45° angle to the normal was positioned in front of the sample. The incident light could reach the sample through a small hole in the mirror. The light back-scattered from the micro-sphere solution was reflected by the mirror in the direction of a CCD detector. Before reaching the detector the back-scattered light was filtered by polarizing optics and focused on the CCD by a lens.
The scattering coefficient for the micro-sphere solution was 11.88 cm-1. A more detailed explanation of the experimental set up is available in .
Cameron’s data is shown in Fig. 2 . To avoid saturation of the CCD a small beam block was introduced in front of the camera lens; this caused the small discontinuity visible in the center of all images.
We modeled Cameron’s experiment with the Monte Carlo programs, and the results are also shown in Fig. 2. Both experimental and simulation figures are composed of the 16 element Mueller matrix images. We normalized all the images by an off-center-value of the m11 images, respectively of Cameron’s and our data; this choice was dictated by the presence of the beam block at the center of Cameron’s images. The Monte Carlo simulation was constructed with the addition of the mirror reflection that flips the image from the left to the right. The Mueller matrix of a mirror was used to adjust the exiting beam polarization. All of the back-reflected shapes are successfully reconstructed by our simulations Figure 2 shows that many of the images (m 11, m 12, m 21, m 2, etc.) are axially symmetric. This symmetry could hide some issues in the Monte Carlo program, and is worth studying in more detail.
A particularly interesting illumination pattern, one that breaks symmetry, is shown in Fig. 3. This problem is also of particular interest to us, as our clinical work on the detection of skin cancer margins  uses a similar asymmetric polarized illumination to avoid glare.
A mono-disperse solution of micro-spheres was illuminated at an angle of approximately 30° by a laser beam of wavelength 632.8 nm and beam diameter 1.2 mm. The incident beam was polarized with a linear polarizer (Ealing, Rocklin, CA), whose orientation was parallel to the optical table (H horizontal), perpendicular to that plane (V vertical), and at +45 degrees to that plane (P). The P orientation was achieved rotating the linear polarizer counterclockwise 45 degrees looking into the incoming beam. The sample was a solution of 2.0μm diameter micro-spheres with index of refraction 1.59 and media index of refraction 1.33 (Duke Scientific Corporation, Palo Alto, CA). Once diluted the micro-sphere solution had a scattering coefficient equal to 11 cm-1.
A set of 9 images was collected with a 12-bit digital camera (Roper Scientific, Trenton, NJ). For every source polarization, the detector was oriented to H, V, and P. The 9 collected pictures HsourceHdetector, Hsource Vdetector, HsourcePdetector, VsourceHdetector, Vsource Vdetector, VsourcePdetector, PsourceHdetector, Psource Vdetector, PsourcePdetector, were combined as in Yao and Wang , to yield the first 9 terms of the backscattered Mueller matrix (Fig. 4 left-hand-side).
The results for the meridian plane Monte Carlo simulations are shown in Fig. 4 on the right-hand-side; identical results were achieved with all three methods. All values of the Mueller matrix were normalized by a maximum value of the m 11 figures. In these images we show only the positive values for all the Mueller matrix terms, the m1 1 values are in between 0 and 0.1, all other values were in between -0.01 and 0.01. The direction cosines at launch were changed to obtain an asymmetrical incident angle of 30°; a non-zero incident angle did not affect the status of polarization of the incident light since the beam reference plane did not change. The programs reproduced all of the shapes correctly. Since we used matched boundaries in the Monte Carlo simulations the glare due to the air/glass/liquid interface visible in the experimental figure on the left-hand-side was not reproduced in the Monte Carlo results.
This simulation set-up can highlight mistakes in the handling of the reference plane rotations, direction cosines, and exiting routines.
3. Imaging of laterally scattered polarized light
Imaging of polarized light scattering in a liquid media is done, almost exclusively, by analyzing the light transmitted or backscattered by the scattering solution. We wanted to observe what happens to the light in the media before it exits the solution in the two aforementioned directions. The images obtained with this experiment show how deep polarized light can penetrate before loosing its original state of polarization and, at the same time, give information on the photons lateral spread. For our particular sphere solution, we observed that light scattering in the lateral direction maintains the original polarization steps for several optical lengths.
A water-tight glass container (7x7x4 cm) was built with wall thickness of 0.25 cm. The container was filled with a solution of 0.5 μm micro-spheres (Duke Scientific Corporation, Palo Alto, CA), Fig.5. The scattering coefficient of the micro-sphere solution was determined with a collimated transmission experiment and was equal to 10.5 cm-1. The index of refraction of the micro-spheres was 1.59, and the index of refraction assumed for water was 1.33. A Helium Neon Laser (Melles-Griot, Carlsbad, CA), at 632.8 nm wavelength and 1.2 mm diameter was oriented so that the beam was parallel to one of the tank sides, and at 1 cm from the wall. A polarizer (Ealing, Rocklin, CA) was used to polarize the incident beam parallel to the optical table.
A 12-bit digital camera (MicroMax, Roper Scientific, Trenton, NJ) was positioned perpendicular to the sidewall of the tank, and focused on the inner wall. An analyzer was located between the camera and the tank. The analyzer was sequentially rotated parallel and perpendicular to the incident beam to acquire a co-polarized and cross-polarized image. Monte Carlo simulations were conducted for this micro-sphere density, size, and beam/tank geometry. The geometry of the experiment was implemented in the code without Fresnel reflectance; every photon reaching the walls of the tank was killed. Results of both experimental and calculated images were normalized to the image maximum value. A background noise of 250 counts was subtracted from the experimental images.
Figure 6 shows the results for the co-polarized experiment and the corresponding Monte Carlo results.
The general shape of the top image (experiment) is successfully reproduced by the bottom Monte Carlo simulation. The air/water and air/glass/water interfaces are visible in the experiment image.
Figure 7 shows exittance along the beam axis (black line in the insert). These values are plotted versus the relative distance from the front wall.
The R2 value for this test was 0.95; better agreement could have been achieved with a Monte Carlo simulation with higher spatial resolution, since this would have avoided some of the sharp peaks deviating from the general behavior.
In Fig. 8 we show the cross-polarized image. When increasing or decreasing the concentration of micro-spheres, we saw different shapes forming. Higher concentration generated well defined lobes in the cross-polarized image; this effect could be helpful when characterizing a micro-sphere solution.
4. Transmission through micro-sphere solutions
In this section we validate the Monte Carlo programs with an experiment measuring polarized transmission through slabs of micro-spheres. This experiment is particularly important because it gives quantitative results. Moreover the experiment shows the experimental error due to the polarizing optics. We measured the polarized light transmission through micro-sphere solutions. The experimental apparatus is shown in Fig. 9, a helium-neon laser (Melles-Griot, Carlsbad, CA) with nominal wavelength 543-nm and beam diameter equal to 1.5 mm was polarized with a polarizer oriented parallel to the optical table (Ealing Electrooptics, plc., Watford, UK). The cross-polarized extinction was equal to 10-5.
An analyzing linear polarizer was placed in front of the detector. An aperture with diameter d1 was located a distance R1 from the rear of the slab, and a detector with diameter d2 was located a distance R2 from the slab. The choices of d1, R1 and R2 were such that the solid angle of collection by the detector equals the solid angle of collection by the aperture (for a point source of diffuse exittance on axis at the slab’s rear surface). For d1 = 0.214 cm, d2 = 1 cm, R1 = 1.5 cm, and R2 = 5.5 cm, the solid angles of collection by the aperture and the detector match and equal 0.016 sr.
The samples were mono-disperse solutions of micro-spheres in water. The considered diameters were 482 nm and 308 nm (Pella Inc, Redding, CA). The length of the cuvette was adjustable allowing measurements of increasing optical thicknesses. The experiment was repeated three times.
The system of Fig. 9 was modeled with our Monte Carlo programs. Each simulation was repeated three times, each time changing the random number seed. All the Monte Carlo runs were conducted with 106 photons. The index of refraction for the spheres in the Monte Carlo program was 1.59 and the assumed index of refraction for water was 1.33. The stock solution scattering coefficient μs was 8 cm-1 for the 482 nm micro-spheres solution and 54 cm-1 for the 308 nm micro-spheres solution. The absorption coefficient was kept equal to 0 in the simulations. Results for experiments (red and blue symbols) and simulations (black symbols) for 482 nm spheres are shown in Fig. 10, and for the 308 nm spheres the results are shown in Fig. 11. All of the data was normalized by a co-polarized transmission measurement through water. The back line corresponds to a heuristic model proposed by Jacques et al. , based on a parameter χ indicative of the depolarizing ability of the solution. χ = 0.98 for the 308 nm solution and χ = 0.8 for the 482 nm solution.
As expected for a parallel polarizer and analyzer, the signal decreases with increasing optical thickness. The opposite is true in the cross-polarized case. The polarized beam randomizes more and more, traveling deeper and deeper in the micro-sphere solution, where some perpendicularly polarized signal is detectable. At an optical thickness τ ~ 10 the beam is completely randomized both parallel and perpendicular signals are equal.
This experiment emphasizes the difficulty implicit in transmission mode experiments with polarized light. It is evident that for certain scatterers and concentrations high precision polarizers are necessary to obtain reliable cross polarized results; an extinction ratio of at least 10-4 is desirable. Small errors in the determination of the solution optical thickness can also generate problems; for this reason we preferred changing the physical length of the slab to changing the solution concentration, since we had a better control of the former.
5. Scattering from poly-disperse solutions of micro-spheres
The scattering from mono-disperse solution of micro-spheres is often used as a gold standard to calibrate an experimental set up, but it is not a very common event in nature. Scattering from poly-disperse particle solutions is more frequent and has received great attention in a variety of fields. Polarized light transport through the atmosphere is a clear example of light scattering from poly-disperse media. The atmospheric optics community has dealt with this issue for a number of years. Allen and Plat  used Lidar measurements to evaluate the depolarization and multiple scattering through cirrus clouds. Zaccanti et al.  wrote a Semi Monte Carlo program to simulate polarized Lidar returns from clouds and poly-disperse media in general, and Bruscaglioni et al. simulated the transmission of a light pulse through a turbid media .
In medicine light scattering from biological media has been used diagnostically in many applications. One example is light transfer through blood; this is of great interest for diagnosis of diabetes and other blood related diseases. Red blood cells (RBC) dominate the scattering and absorption in blood and can be considered a poly-disperse solution of particles. The modeling of this problem was attempted by many. He et al.  used a Finite Difference Time Domain Method, while Tsinopoulos et al. studied aggregated RBC in a rouleau geometry . Others used unpolarized Monte Carlo models  to model light scattering through diluted solutions of blood.
All three Monte Carlo programs can be modified to treat the polarized light scattering from poly-disperse solutions. A photon traveling in a poly-disperse solution of micro-spheres will encounter particles of different radii depending on the solution particle size distribution and scattering coefficient μs, Fig. 12.
In our implementation the scattering coefficient μs is fixed at initialization; consequently, the mean photon step size does not change during the simulation. The programs are structured as described in the previous paper with Launch, Move, Drop, and Scatter sections. The only difference is that the particle size parameter changes for every scattering event in accordance with a predetermined particle size distribution. Consequently the corresponding terms of the scattering matrix s11, s12, s33, and s34 will also change at every scattering event.
Several distributions were considered, starting with simple uniform distributions to more complex. If the particle size distribution was particularly complex the rejection method  was used to select the radii. Monte Carlo runs were compared with the Adding Doubling model of Evans et al.  for upward and downward radiation of polarized light through a slab. The characteristics of Evans’ code were described in section 7 of the previous paper  and will not be restated here for brevity.
In the poly-disperse case Evans’ program calculates the first (I) and second (Q) term of the total backscattered Stokes vector for a modified Gamma distribution of particles N(r).
where r is the particle radius that is allowed to vary between two values, rmin and rmax.
Evans’ results are compared to the Monte Carlo results for different particle size distributions. Table 1 shows the results for transmission through a slab of thickness 4/μs.
Table 2 shows the results for total reflectance from the slab. In the simulations the spheres index of refraction was 1.59 and the media index of refraction was 1.0.
The wavelength was 0.6328 μm, with light delivered as a pencil beam normal to the slab. The incident light was unpolarized. All Monte Carlo simulations were conducted with 1 million photons.
The agreement between Monte Carlo and Adding Doubling results is good (error<1%) for small sphere solutions, and becomes progressively worse as the particle size increases as in the last case in table 1 and 2. This is possibly due to the fact that the number of integration points in the particle distribution influences the results both in Adding Doubling and Monte Carlo, and that subdivision process is different in the two cases.
We investigated the simple case of light traveling through poly-disperse solutions of micro-spheres with uniform size distribution. A slab-geometry was considered Fig. 13. The optical thickness (τ= Lμs) of the slab was varied to study the effect of scattering on the transmitted Stokes vector.
In all the simulations the detector was infinite and the polarization reference plane was the plane xy perpendicular to the slab. The simulations were conducted on four different uniform distributions where the upper and lower radius (r) limits were respectively 0.01 μm < r < 0.01 μm, 0.1 μm < r < 0.5 μm, 0.1 μm < r < 1.0 μm, and 0.25 μm < r < 0.5 μm. In all cases the index of refraction of the spheres was 1.59 and the index of refraction of the surrounding medium was 1.33. The incident light was a pencil beam polarized parallel to the reference plane with wavelength equal to 0.543 μm. The results of the simulations are shown in Fig. 14 and represent the degree of linear polarization (DLP) of light transmitted through the slab as a function of slab-thickness. Simulations were also conducted with mono-disperse solutions of micro-spheres whose anisotropy parameter g was equal to the average anisotropy gav of the poly-disperse simulations. The results of Fig. 14 show that in the simple case of uniformly distributed micro-spheres the degree of linear polarization for a uniformly distributed poly-disperse solution of anisotropy gav behaves as a mono-disperse solution of micro-spheres of identical anisotropy g. For example the DLP for a mono-disperse solution of micro-spheres of radius 0.063 μm and anisotropy g equal to 0.16 behaved just as the DLP for poly-disperse solutions, where the particle radii were uniformly distributed between 0.01 μm and 0.1 μm, and the average anisotropy was also 0.16.
In reality not many materials are composed of uniformly distributed micro-spheres. Moreover, as previously shown in Table 1 and 2, exponential distributions behave differently from uniform ones. These results are valid only for the geometry of Fig. 13 and the aforementioned experimental parameters, but are nevertheless interesting in relating the influence of anisotropy in the case of total transmission and incident linear polarization.
The polarized light Monte Carlo programs described in our previous paper  were applied to a variety of experimental geometries. Frequently used set-ups such as backscattering from micro-sphere solutions and transmission through slabs were successfully reproduced.
We were able to reproduce the Mueller matrix of light backscattered from a solution of micro-spheres obtaining identical patterns to the work published by Cameron et al. A skewed illumination-geometry was also implemented, both experimentally, and in the simulations. This geometry is particularly interesting since some groups are using skewed illumination angles and normal collection to avoid glare. The skewed illumination also guarantees the breaking of symmetry for some of the Monte Carlo images; this can help in thoroughly test the programs.
Another experimental geometry of interest is the transmission of polarized light through slabs of scattering solutions. Our Monte Carlo programs can easily treat this case and we were able to successfully match experiments conducted with different sphere sizes.
Finally we modified the Monte Carlo programs to handle the scattering from poly-disperse solutions of micro-spheres. Comparison between the Monte Carlo runs and a polarized Adding Doubling code showed good agreement for small particle sizes.
The programs were used to study polarized light transmission through slabs of poly-disperse solutions of micro-spheres in the limited case of uniform particle size distribution. It was noticed that poly-disperse solutions of average anisotropy value gav, behave just like mono-disperse solutions of equal anisotropy. One consequence of this finding is that the time of execution of the mono-disperse Monte Carlo is much shorter than the corresponding poly-disperse case. We treated a very simple case where the media was composed of spherical particle of uniformly distributed sizes noting that most naturally occurring media do not have such a simplified distribution.
Jessica Ramella-Roman wishes to thank Brent D. Cameron for the use of the Mueller matrix experimental data and K. F. Evans for its feedback on its polarized Adding Doubling program.
References and Links
1 . J.C. Ramella-Roman , S.A. Prahl , and S.L. Jacques , “ Three Monte Carlo programs of polarized light transport into scattering media: part I ,” Opt. Express 13 , 4420 – 4438 , ( 2005 ), http://www.opticsexpress.org/abstract.cfm?URI=OPEX-13-12-4420 [CrossRef] [PubMed]
2 . 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]
3 . B. D. Cameron , M. J. Rakovic , M. Mehrubeoglu , G. W. Kattawar , S. Rastegar , L. V. Wang , and G. Cote , “ Measurement and calculation of the two-dimensional backscattering Mueller matrix of a turbid medium ,” Opt. Lett. 23 , 485 – 487 , 1998 . [CrossRef]
4 . M. J. Rakovic , G. W. Kattawar , M. Mehrubeoglu , B. D. Cameron , L. V. Wang , S. Rastegar , and G. L. Cote , “ Light backscattering polarization patterns from turbid media: theory and experiment ,” Appl. Opt. 38 , 3399 – 3408 , ( 1999 ). [CrossRef]
5 . Daniel CÔté and I. Alex 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 . http://www.opticsexpress.org/abstract.cfm?URI=OPEX-13-1-148 [CrossRef] [PubMed]
6 . M. Mehrubeoglu , 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]
7 . M. Xu , “ Electric field Monte Carlo simulation of polarized light propagation in turbid media ”, Opt. Express 26 , 6530 – 6539 , ( 2004 ). http://www.opticsexpress.org/abstract.cfm?URI=OPEX-12-26-6530 [CrossRef]
8 . K.F. Evans and G.L. Stephens , “ A new polarized atmospheric radiative transfer model ,” J.Quant.Spectrosc. Radiat. Transfer. 46 , 413 – 423 , ( 1991 ). [CrossRef]
9 . The experimental data gently provided by Dr. Cameron.
11 . G. Yao and L-H. Wang , “ Two-dimensional depth-resolved Muller matrix characterization of biological tissue by optical coherence tomography ,” Opt. Letters 24 , 537 – 539 , ( 1999 ). [CrossRef]
13 . G. Zaccanti , P. Bruscaglioni , M. Gurioli , and P. Sansoni , “ Laboratory simulations of lidar returns from clouds: experimental and numerical results ,” Appl. Opt. 32 , 1590 – 1597 , ( 1993 ). [CrossRef] [PubMed]
14 . P. Bruscaglioni , G. Zaccanti , and Q. Wei , “ Transmission of a pulsed polarized light beam through a thick turbid media: numerical results ,” Appl. Opt. 32 , 6142 – 6150 , ( 1993 ). [CrossRef] [PubMed]
15 . J. He , A Karlsson , J. Swartling , and S. Anderson-Engles , “ Light scattering by multiple red blood cells ,” J. Opt. Soc. Am. A. 21 , 1953 – 1961 , ( 2004 ). [CrossRef]
17 . L. Whang , S.L. Jacques , and L. Zheng , “ MCML - Monte Carlo modeling of light transport in multi-layered tissues ,” Comput. Methods Programs Biomed. 47 , 131 – 46 , ( 1995 ). [CrossRef]
18 . W. H. Press , S. A. Teukolsky , W. T. Vetterling , and B. P. Flannery , Numerical Recipes in C, the art of Scientific Computing , ( Cambridge University Press , 1992 ).