Electric field Monte Carlo simulation of polarized light propagation in turbid media

A Monte Carlo method based on tracing the multiply scattered electric field is presented to simulate the propagation of polarized light in turbid media. Multiple scattering of light comprises a series of updates of the parallel and perpendicular components of the complex electric field with respect to the scattering plane by the amplitude scattering matrix and rotations of the local coordinate system spanned by the unit vectors in the directions of the parallel and perpendicular electric field components and the propagation direction of light. The backscattering speckle pattern and the backscattering Mueller matrix of an aqueous suspension of polystyrene spheres in a slab geometry are computed using this Electric Field Monte Carlo (EMC) method. An efficient algorithm computing the Mueller matrix in the pure backscattering direction is detailed in the paper. © 2004 Optical Society of America OCIS codes: (170.5280) Photon migration; (030.5620) Radiative transfer (290.4210) Multiple scattering; (290.1350) Backscattering (290.7050) Turbid media; (030.6140) Speckle References and links 1. A. Ishimaru, Wave propagation and scattering in random media, I and II (Academic, New York, 1978). 2. A. Yodh and B. Chance, “Spectroscopy and imaging with diffusing light,” Phys. Today 48(3), 38–40 (1995). 3. S. K. Gayen and R. R. Alfano, “Emerging optical biomedical imaging techniques,” Opt. Photon. News 7(3), 17–22 (1996). 4. S. R. Arridge, “Optical tomography in medical imaging,” Inverse Problems 15, R41–R93 (1999). 5. S. Chandrasekhar, Radiative transfer (Dover, New York, 1960). 6. K. F. Evans and G. L. Stephens, “A new polarized atmospheric radiative transfer model,” J. Quant. Spectrosc. Radiat. Transfer 46, 413–423 (1991). 7. A. D. Kim and M. Moscoso, “Chebyshev Spectral methods for radiative transfer,” SIAM J. Sci. Comput. 23, 2074–2094 (2002). 8. A. D. Kim and J. B. Keller, “Light propagation in biological tissue,” J. Opt. Soc. Am. A 20, 92–98 (2003). 9. G. W. Kattawar and G. N. Plass, “Radiance and polarization of multiple scattered light from haze and clouds,” Appl. Opt. 7, 1519–1527 (1968). 10. I. Lux and L. Koblinger, Monte Carlo Particle Transport Methods: Neutron and Photon Calculations (CRC Press, Boca Raton, Fla., 1991). 11. 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(30), 6535– (1992). 12. P. Bruscaglioni, G. Zaccanti, and Q. Wei, “Transmission of a pulsed polarized light beam through thick turbid media: numerical results,” Appl. Opt. 32(30), 6142–6150 (1993). 13. M. J. Rakovic, G. W. Kattawar, M. Mehrbeolu, 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(15), 3399–3408 (1999). 14. S. Bartel and A. H. Hielscher, “Monte Carlo Simulations of the Diffuse Backscattering Mueller Matrix for Highly Scattering Media,” Appl. Opt. 39(10), 1580–1588 (2000). (C) 2004 OSA 27 December 2004 / Vol. 12, No. 26 / OPTICS EXPRESS 6530 #5731 $15.00 US Received 15 November 2004; revised 13 December 2004; accepted 14 December 2004 15. M. Moscoso, J. B. Keller, and G. Papanicolaou, “Depolarization and blurring of optical images by biological tissue,” J. Opt. Soc. Am. A 18(4), 948–960 (2001). 16. H. H. Tynes, G. W. Kattawar, E. P. Zege, I. L. Katsev, A. S. Prikhach, and L. I. Chaikovskaya, “Monte Carlo and Multicomponent Approximation Methods for Vector Radiative Transfer by use of Effective Mueller Matrix Calculations,” Appl. Opt. 40(3), 400–412 (2001). 17. B. Kaplan, G. Ledanois, and B. villon, “Mueller Matrix of Dense Polystyrene Latex Sphere Suspensions: Measurements and Monte Carlo Simulation,” Appl. Opt. 40(16), 2769–2777 (2001). 18. X. Wang and L. V. Wang, “Propagation of polarized light in birefringent turbid media: A Monte Carlo study,” J. Biomed. Opt. 7, 279–290 (2002). 19. G. W. Kattawar, M. J. Raković, and B. D. Cameron, “Laser backscattering polarization patterns from turbid media: theory and experiments,” in Advances in Optical Imaging and Photon Migration, J. G. Fujimoto and M. S. Patterson, eds., vol. 21 of OSA TOPS, pp. 105–110 (1998). 20. J. C. Ramella-Roman, “Imaging skin pathologies with polarized light: empirical and theoretical studies,” Ph.D. thesis, OGI School of Science & Engineering at Oregon Health & Science University (2004). 21. H. C. van de Hulst, Light scattering by small particles (Dover, New York, 1981). 22. C. F. Bohren and D. R. Huffman, Absorption and scattering of light by small particles (John Wiley & Sons, 1983). 23. R. Y. Rubinstein, Simulation and the Monte Carlo method (John Wiley & Sons, 1981). 24. J. von Neumann, “Various techniques used in connection with random digits,” J. Res. Natl. Bur. Stand. 5, 36–38 (1951). 25. W. H. Press, S. A. Teukolsky, W. T. Vetterling, and B. P. Flannery, Numerical recipes in C (Cambridge university press, 1996). 26. P.-E. Wolf and G. Maret, “Weak Localization and Coherent Backscattering of Photons in Disordered Media,” Phys. Rev. Lett. 55(24), 2696–2699 (1985). 27. M. P. V. Albada and A. Lagendijk, “Observation of Weak Localization of Light in a Random Medium,” Phys. Rev. Lett. 55(24), 2692–2695 (1985). 28. E. Akkermans, P. E. Wolf, and R. Maynard, “Coherent backscattering of light by disordered media: analysis of the peak line shape,” Phys. Rev. Lett. 56(14), 1471–1474 (1986). 29. J. W. Goodman, “Statistical properties of laser speckle patterns,” in Laser speckle and related phenomena, J. C. Dainty, ed., pp. 9–75 (Springer-Verlag, Berlin, 1975). 30. D. S. Saxon, “Tensor Scattering Matrix for the Electromagnetic Field,” Phys. Rev. 100(6), 1771–1775 (1955). 31. B. D. Cameron, M. J. Rakovic, M. Mehrbeoglu, G. W. Kattawar, S. Rastegar, L. V. Wang, and G. L. Cot, “Measurement and calculationof the two-dimensional backscattering Mueller matrix of a turbid medium,” Opt. Lett. 23(7), 485–487 (1998). 32. I. Berezhnyy and A. Dogariu, “Time-resolved Mueller matrix imaging polarimetry,” Opt. Exp. 12(19), 4635– 4649 (2004). 33. EMC is available at http://www.sci.ccny.cuny.edu/ minxu.


Introduction
The propagation of polarized light in turbid media is fundamental to many practical applications of considerable interest including remote sensing of clouds and imaging of colloidal suspensions and biological materials [1][2][3][4].Due to the lack of analytical solutions to radiative transfer [5] of polarized light within a bounded medium, numerical solutions [6][7][8] of the radiative transfer equation and Monte Carlo simulations [9][10][11][12][13][14][15][16][17][18] of propagation of polarized light in turbid media have been pursued extensively and applied to characterization of clouds, particle suspensions, and biological materials.
Most Monte Carlo methods trace the Stokes vector I = (I, Q,U,V ) T in simulation where r , E l and E r are two orthogonal complex electric field components perpendicular to the propagation direction, the superscript "T " denotes transpose, and means ensemble average.Light scattering involves a rotation of the Stokes vector to a local scattering reference frame and the multiplication of the Stokes vector by the 4 × 4 Mueller matrix which prescribes how polarized light is scattered by an isolated particle in that reference frame.The Stokes vector is summed up at the detector assuming the detected light is incoherent.Many implementations of the above Monte Carlo approach to simulate polarized light propagation in turbid media have appeared in the literature [9][10][11][12][13][14][15][16][17][18].Some recent work includes [13,14,[16][17][18].Symmetry considerations [19], randomly polarized incident light [17], next event point estimators [17] and other techniques [14,16] have been investigated to improve the efficiency of the Monte Carlo computation.Various implementations of the above Monte Carlo approach have also been compared [20].
In this article, we present a Monte Carlo method based on tracing the multiply scattered electric field to simulate the propagation of polarized light in turbid media.Multiple scattering of light comprises a series of updates of the parallel and perpendicular components of the complex electric field with respect to the scattering plane by the amplitude scattering matrix and rotations of the local coordinate system spanned by the unit vectors in the directions of the parallel and perpendicular electric field components and the propagation direction of light.The phase of light accrues from phase delays due to the interaction of light with scatterers and propagation through the host medium.In contrast with the conventional Monte Carlo approach based on tracing the Stokes vector, the Electric Field Monte Carlo (EMC) method traces the electric field and phase of light and makes it possible to simulate also coherent properties of multiply scattered light.The algorithm of EMC is straightforward and can be easily adapted to simulate the propagation of polarized light in optically active medium.
After outlining the theoretical formalism of the Electric Field Monte Carlo method in section 2, we present EMC computations of the backscattering speckle pattern and the backscattering Mueller matrix of an aqueous suspension of polystyrene spheres in a slab geometry in section 3. Special considerations to improve the efficiency of computing the Mueller matrix in the pure backscattering direction are detailed in section 3. We conclude in section 4.

Theoretical formalism
Light scattering by small particles is succinctly described by the amplitude scattering matrix [21].For spherical or randomly-oriented aspheric scatterers, the amplitude scattering matrix is diagonal and depends only on the scattering angle θ due to the symmetry [22].The parallel and perpendicular components of the electric field with respect to the scattering plane are scattered according to S j (θ ) where j = 2, 1, respectively [21].
The scattering of photons takes a simple form in the local orthonormal coordinate system (m, n, s) where m and n are the unit vectors in the directions of the parallel and perpendicular components of the electric field with respect to the scattering plane of the previous scattering event and s is the photon's propagation direction prior to the current scattering [see Fig. 1].The propagation direction s of the photon after the current scattering is given by: s = m sin θ cos φ + n sin θ sin φ + s cos θ where θ and φ are the scattering and azimuthal angles of the current scattering respectively.The current scattering plane is the plane spanned by s and s .The unit vectors in the directions of the parallel and perpendicular electric fields with respect to the current scattering plane are given by e 1 = m cos φ + n sin φ (2) e 2 = −m sin φ + n cos φ prior to scattering and   We can now summarize the updating rules of the local coordinate system and the electric field in EMC.For each scattering with the scattering angle θ and the azimuthal angle φ , the local coordinate system is updated by with and the electric field by We have introduced an extra factor F(θ , φ ), the scattered wave intensity propagating along the (θ , φ ) direction, given by to normalize the light intensity.The intensity of the incident light 1 is assumed unity in Eq. ( 8).The scattered light intensity in the series of scattering events.Absorption of light, if any, is accounted for by adjusting the photon weight as in the conventional Monte Carlo simulations [23].
The orthogonal matrix A rotates the local coordinate system (m, n, s) each time the photon is scattered by a scatterer.The electric field is updated simultaneously by the matrix B. The free path between consecutive scattering events is sampled in the same fashion as that in a conventional Monte Carlo method.The state of a multiply scattered photon is specified by the final local coordinate system m (n) , n (n) , s (n) after consecutive rotations, the final complex electric field components E (n) 1,2 after consecutive updates, and the optical path l inside the host medium where n denotes the number of scattering events the photon have experienced before being detected.The detected electric field is given by where k = 2π/λ is the wave number, λ is the wavelength of light in the host medium, and we have assumed a convention that the temporal dependence of a plane wave of frequency ω is exp(−iωt).The phase of the detected photon accrues from phase delays due to the interaction of light with scatterers and propagation through the host medium.The photon weight w, unity at incidence, is multiplied by the albedo of the scatterer at each scattering event.Once the photon hits a detector, the electric field at the detector is increased by w 1/2 E d and the Stokes vector is increased by the Stokes vector computed from the electric field w 1/2 E d .
Propagation of multiply scattered light is describable by the radiative transport equation in which any interference effects are neglected [5].The conventional Monte Carlo methods based on tracing the Stokes vector of light solve the radiative transfer equation numerically and do not include any correlation effect of multiply scattered light.By tracing the electric field associated with a wave packet, EMC provides much more detailed information about the propagation of multiply scattered light than just the transfer of energy.EMC should be regarded as a method to sample the probability distribution of the electric field at selected spatial positions where the detectors are located.Detectors of finer resolution than the speckle size are required to resolve the interference pattern of light well.Detectors of a larger cell size will smear the interference pattern.It should be noted that EMC only probes one instantaneous picture of the disordered medium when accumulating the electric field coherently while it yields the detected light ensemble averaged over all pictures of the disordered medium when accumulating the Stokes vector incoherently.This point is addressed in more detail in section 3.1.
One key step in the Monte Carlo simulation of polarized light is the sampling of the scattering angles (θ , φ ) which distribute according to the normalized phase function p(θ , φ ) = F(θ , φ )/πQ sca x 2 where Q sca is the scattering efficiency, x = ka is the size parameter, and a is the radius of the particle.The rejection technique [24,25] has been widely used to sample such a distribution function.The procedure is to choose a doublet (µ ≡ cos θ , φ ) where µ and φ are uniformly distributed over [−1, 1] and [0, 2π] respectively and a random number f uni-formly distributed over [0, max θ ,φ p(θ , φ )], the doublet (µ, φ ) is accepted as the scattering and azimuthal angles if f < p(arccos µ, φ ), otherwise generate a new doublet (µ, φ ), a new f and start over.Each trial involves one Mie calculation of the amplitude scattering matrix at that trial scattering angle.The number of rejections per scattering event is large (on the order of one for Rayleigh particles and can be much more significant for larger particles).This limits the efficiency of the Monte Carlo simulation.Mie calculations of the amplitude scattering matrix are hence usually performed on a fixed set of scattering angles in advance to generate a table of scattering matrices and reduce the computation load.
We generate the scattering angle θ by drawing of one random number uniformly distributed over [0, 1] and looking up an inverse table of the marginal distribution function which does not depend on the electric field and is pre-calculated before simulation.The azimuthal angle φ is then obtained by the rejection method for the conditional probability p(φ |θ ) = p(θ , φ )/p(θ ).One Mie calculation for the looked up scattering angle is required if the table of scattering matrices is not generated in advance for each scattering event.
A different sampling strategy is to sample the scattering angle θ according to p(θ ) while to sample φ uniformly distributed over [0, 2π].At the same time, we modify F(θ , φ ) in the update rule of the electric field (6) to that the light intensity is no longer conserved in the series of scattering events.This strategy saves the rejection sampling of the angle φ but the simulation result is prone to be contaminated by hotspots and has a larger variance compared to the first strategy because it is unfavorable to the statistical variance when photons of varying weights, rather than equal weights, are accumulated by the detector.The second sampling strategy is not used in simulation.

Backscattering speckle pattern
Due to the large ratio between the velocities of light and of the scatterers, coherent light probes an instantaneous picture of the disordered medium and realizes speckle patterns when the multiple scattered light emerges from within the medium.The α = x, y, z component of the outgoing light intensity I α (θ , φ ) in direction (θ , φ ) on the boundary of the medium comprises a multitude of independently-phased additive complex electric fields in that direction (outside the coherent backscattering cone [26][27][28]) and follows a Rayleigh distribution [29] p(I α ) = 1 I α exp(−I α / I α ) where I α is the mean intensity.The normalized light intensity η = I α / I α follows the exponential distribution exp(−η).
We perform EMC to study the plane wave light backscattering from an aqueous solution of polystyrene spheres inside a slab.Incident light is polarized in the x direction and propagates in the z direction (normal to the surface of the slab).The size of the polystyrene sphere is 0.46µm.The wavelength of the incident wave is λ vac = 0.52µm in vaccum.The refractive indices of the polystyrene sphere and water is n sct = 1.59 and n bg = 1.33, respectively.The mean scattering length of light inside the solution is l s = 2.80µm with 2πn bg l s /λ vac = 45.The thickness of the slab is L z = 4l s .Total 5 × 10 8 photons are launched into the medium.Both backscattering electric field and Stokes vector are recorded versus the direction (θ , φ ) of the backscattered light where θ is the angle away from the surface normal (0 ≤ θ ≤ π/2) and φ is the azimuthal angle (0 ≤ φ ≤ 2π).The ranges of the angles θ and φ are uniformly divided into 1125 and 360 bins in simulation, respectively.EMC, like an experiment probing one static random medium, only probes one instantaneous picture of the disordered medium.The ensemble average is realized in experiments, for example, through the movement of the scatterers in the medium or sampling different regions of the random media.EMC, as a numerical experiment, can record both the electric field and the Stokes vectors in simulation.The recorded electric field added coherently yields the emergent light from the one instantaneous picture of the disordered medium.The recorded Stokes vector added incoherently, on the other hand, yields the emergent light ensemble averaged over all pictures of the disordered medium.
Figure 2(a) displays the x component I x / I x where I x = |E x | 2 and I x is estimated by the mean of the first and second Stokes vectors.The appearance of speckles is obvious in Fig. 2(a).The first order statistics about I x / I x can be computed from its histogram.This yields a negative exponential distribution exp(−η) as expected [see Fig. 2(b)].The other two y and z components of light behave similarly and not displayed.

Backscattering Mueller matrix
We then compute the backscattering Mueller matrix from a pencil beam incident normally on an aqueous solution of polystyrene spheres inside a slab.The size of the polystyrene sphere is 2µm.The wavelength of the incident wave is λ vac = 0.63µm in vacuum.The refractive indices of the polystyrene sphere and water is n sct = 1.59 and n bg = 1.33, respectively.The thickness of the slab is L z = 4l s .Total 3 × 10 8 photons are launched into the medium.
To improve its efficiency of the Monte Carlo simulation, we combined symmetry considerations [19], randomly polarized incident light [17] and the next event point estimator [17] in computing the backscattering Mueller matrix.The computation time is less than 2 hours for each 10 8 photons launched using one Pentium III 1GHz CPU.In computation of the Stokes vector in Monte Carlo simulations, one should note that the Stokes vector upon which Mueller matrix is defined depends on the reference system used.The backscattered Stokes vector is defined in the reference system whose x y z axes coincide with −x, y, −z axes of the reference system of the incident Stokes vector respectively.
In the backscattering geometry investigated here where the incident light is normal to the surface of the medium and the outgoing beam is exactly backscattering, the Mueller matrix M bs (ρ, φ ) of the medium shall satisfy the following relation [19,30] where ρ, φ is the polar coordinate of the position of the outgoing beam and R(φ ) is the rotation matrix The reduced Mueller matrix M bs 0 (ρ) ≡ M bs (ρ, φ = 0) relates where I i,o are the incident and outgoing Stokes vectors with respect to the φ = 0 plane.In our simulation, the incident electric field is randomly polarized as E i = cos αe −iβ x + sin αe iβ ŷ where α ∈ (0, π/2) and β ∈ (0, π) are uniform random numbers.The corresponding incident Stokes vector with respect to the xy axes (i.e., the y = 0 plane) is I i = (1, cos 2α, sin 2α cos 2β , sin 2α sin 2β ) T .The incident Stokes vector with respect to the φ = 0 plane is given by I i = R(φ )I i .The ensemble average of I i (I i ) T over the random variates α and β yields Noting that the inverse of the matrix I i (I i ) T is given by we obtain from Eq. ( 12).The backscattering Mueller matrix is then computed from the reduced Mueller matrix by Eq. ( 10).The computed backscattering Mueller matrix is displayed in Fig. (3).Each matrix element is given as a two-dimensional image of the surface, 20l s × 20l s in size, with the laser being incident in the center.The displayed Mueller matrix has been normalized by the maximum light intensity of the (1, 1) element of M bs .Only 7 elements of the Mueller matrix are independent.The symmetrical relation between different elements of the Mueller matrix has been considered by Kattawa et.al. [19,31].The symmetry of the backscattering Mueller matrix in Fig. (3) agrees with Kattawa et.al. [13,19,31] and Berezhnyy et.al. [32].The elements M bs 14 and M bs 41 vanish as predicted by the theory [13,19].It is more instructive, though, to display the reduced Mueller matrix M bs 0 (ρ).Each element of the reduced Mueller matrix is given as a one-dimensional curve versus the distance ρ from the origin in Fig. (4).The reduced backscattering Mueller matrix is found to be 2 × 2 block diagonal.The first two elements of the Stokes vector I 1,2 are then decoupled from the rest two elements of the Stokes vector I 3,4 .This means, for example, the backscattered light due to a normally incident light linearly polarized in the φ = 0 plane (I i = (1, 1, 0, 0) T ) has no circular polarization component (the fourth element of Stokes vector I o = M bs 0 (ρ)I i is always zero) with respect to the φ = 0 plane.The circular polarization component of this backscattered light is no longer zero with respect to a different reference frame rather than the φ = 0 plane.

Conclusion
In conclusion, we have presented a Monte Carlo method based on tracing the electric field for simulating polarized light propagation in turbid media.The Electric Field Monte Carlo method has been successfully used to study the backscattering speckle pattern and the backscattering Mueller matrix of an aqueous suspension of polystyrene spheres in a slab geometry.The tracing of the electric field in simulation makes the Electric Field Monte Carlo method possible to simulate also coherent properties of light.The EMC source code will be put in the public

Fig. 1 .
Fig.1.A photon moving along s is scattered to s with a scattering angle θ and an azimuthal angle φ inside a local coordinate system spanned by orthonormal bases (m, n, s). e 1,2 and e 1,2 are the unit vectors parallel and perpendicular to the current scattering plane spanned by s and s prior to and after scattering.The local coordinate system (m, n, s) is rotated to (m , n , s ) after scattering.

e 2 = e 2 after
scattering, respectively.The local coordinate system (m, n, s) is rotated to (m , n , s ) whose m = e 1 and n = e 2 after scattering.The incident electric field E = E 1 m + E 2 n is scattered to E = E 1 m + E 2 n whose parallel and perpendicular components are given by E 1 = S 2 E • e 1 and E 2 = S 1 E • e 2 , respectively.

Fig. 2 .
Fig. 2. (a) Speckle pattern formed by the angular-resolved backscattering light.(b) Normalized speckle I x / I x follows a negative exponential distribution.

Fig. 3 .
Fig. 3. Backscattering Mueller matrix for the slab.All 4 × 4 matrix element are displayed as a two-dimensional image of the surface, 20l s × 20l s in size, with the laser being incident in the center.The displayed Mueller matrix has been normalized by the maximum light intensity of the (1, 1) element.The parameters of the slab is given in the text.

Fig. 4 .
Fig. 4. Reduced backscattering Mueller matrix for the slab.All 4 × 4 elements of the reduced Mueller matrix is displayed as a one-dimensional curve versus the distance ρ/l s from the origin.The reduced backscattering Mueller matrix is 2 × 2 block diagonal.