1. 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–4]. Due to the lack of analytical solutions to radiative transfer[5] of polarized light within a bounded medium, numerical solutions[6–8] of the radiative transfer equation and Monte Carlo simulations[9–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 I=〈|E_l |²+|E_r |²〉, Q=〈|E_l |²-|E_r |²〉, U=〈$E_l^{*}$ E_r +E_l$E_r^{*}$ 〉, V=-i〈$E_l^{*}$ E_r -E_l$E_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–18]. Some recent work includes [13, 14, 16–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.

2. 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:

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

prior to scattering and

e′ ₂=e ₂

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

 

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.

Download Full Size | PDF

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

with

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 |E ₁|²+|E ₂|²=1 is assumed unity in Eq. (8). The scattered light intensity |E′ ₁|²+|E′ ₂ |² is then conserved (equal to unity) 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 ⁽ⁿ⁾, s ⁽ⁿ⁾) after consecutive rotations, the final complex electric field components $E_{1,2}^{\left(n\right)}$ 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 E _d =[$E_1^{\left(n\right)}$ m ⁽ⁿ⁾+$E_2^{\left(n\right)}$ n ⁽ⁿ⁾]exp(ikl) 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 ² 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 uniformly 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 F′(θ)=|S ₂|²+|S ₁|²/2 such 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.

3. Results and Discussion

3.1. 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–28]) and follows a Rayleigh distribution[29] $p\left(I_{\alpha }\right)=\frac{1}{〈I_{\alpha }〉}\exp (-I_{\alpha }⁄〈I_{\alpha }〉)$ 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⁸ 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.

 

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

Download Full Size | PDF

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

3.2. 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⁸ 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⁸ 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_0^{\text{bs}}$(ρ)≡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,cos2α, sin2α cos2β, sin2α sin2β) ^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_{14}^{\text{bs}}$ and $M_{41}^{\text{bs}}$ vanish as predicted by the theory[13, 19].

 

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.

Download Full Size | PDF

It is more instructive, though, to display the reduced Mueller matrix $M_0^{\text{bs}}$ (ρ). 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_0^{\text{bs}}$ (ρ)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.

4. 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 domain[33].

 

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.

Download Full Size | PDF