## Abstract

A method for directly simulating coherent backscattering of polarized light by a turbid medium has been developed based on the Electric field Monte Carlo (EMC) method. Electric fields of light traveling in a pair of time-reversed paths are added coherently to simulate their interference. An efficient approach for computing the electric field of light traveling along a time-reversed path is derived and implemented based on the time-reversal symmetry of electromagnetic waves. Coherent backscattering of linearly and circularly polarized light by a turbid medium containing Mie scatterers is then investigated using this method.

© 2008 Optical Society of America

## 1. Introduction

The study of light propagation in turbid media is a fundamental problem in applied science that has vast applications in many different fields [1–4]. Inside turbid media, light is multiply scattered and the behavior of light propagation is usually described by radiative transfer, where the interference effects of light are ignored. The interference of light, however, does survive inside the coherent backscattering cone around the exact backscattering direction for incident coherent light [5, 6]. This is the celebrated phenomenon of coherent backscattering (CBS), also known as weak localization of light.

Formally, CBS is the consequence of constructive interference of photons which travel in a pair of time-reversed paths. The phase difference between two partial waves traveling in a pair of time-reversed paths annihilates when the backscattered light remits in the direction close to that of exact backscattering. Inside the backscattering cone, the intensity of light needs to be written as |**E**
_{out}+**E**
^{Reν}
_{out}|^{2}, where **E**
_{out},**E**
^{Reν}
_{out} are the electric fields of the two partial waves, roughly double the intensity |**E**
_{out}|^{2}+|**E**
^{Reν}
_{out}|^{2} when backscattered light is assumed to be incoherent.

Most conventional Monte Carlo simulations trace the Stokes vector [7–11]:

where

and *E _{x}*,

_{y}are the complex electric field components perpendicular to each other and the propagation direction of light. The superscript

*T*refers to the transpose of the matrix and 〈 〉 denotes ensemble average. Recently, the Electric field Monte Carlo (EMC) method was developed by us to simulate polarized light propagating through turbid media [12,13]. The electric field, rather than the Stokes vector, of light is traced in EMC. The complete phase of light is accrued from both light scattering by particles and propagation through the host medium. EMC is hence most amenable to study coherence phenomena of multiple scattering light, in particular, CBS.

In this paper, we first derive a relation between the electric fields of light traveling in a pair of time-reversed paths based on the time reversal symmetry of electromagnetic waves. This is done in section 2 and the electric field of light traveling in the time-reversed path is shown to be computed from that of light traveling in the forward path, greatly reducing the computation time. The coherent backscattering of circularly polarized light and linearly polarized light is then studied in section 3. Their differences in the dependence of the coherent backscattered light on the azimuthal angle and the dependence of the coherent enhancement on the polar angle are investigated. Concluding remarks are presented in section 4.

## 2. Theory

The underlying equation for the updating rule of the electric field of scattered light [14] in EMC is:

where **E** is the complex electric field with components parallel and perpendicular to the previous scattering plane, **E**′ is the complex electric field with components parallel and perpendicular to the present scattering plane, *S* is the amplitude scattering matrix dependent on the scattering angle *θ* between the incoming propagation direction and outgoing propagation direction, and *R* is the rotation matrix dependent on the angle *ϕ* between the incoming perpendicular electric field component and the outgoing perpendicular electric field component. The rotation matrix *R* rotates the reference frame *ϕ* degrees azimuthally to align the incoming perpendicular electric field component to the normal of the present scattering plane. Equation (1) can be explicitly written as the following for Mie scatterers:

where *E*
_{1} and *E*
_{2} are the complex parallel and perpendicular components of the incoming electric field and *E*′_{1} and *E*′_{2} are the complex electric field components after the scattering event. The properties of the scattering particle are contained in the *S* matrix: *S*
_{1} and *S*
_{2} dependent on the refractive indices inside and outside the particle as well as the size of the particle and the wavelength of the incident light.

## 2.1 Electric field of the partial waves in the forward and time-reversed paths

The local coordinate system described by three vectors **m**, **n**, and **s** is rotated after every scattering event. Here **m** is the direction of the parallel electric field component and **n** is the direction of the perpendicular electric field component. The vector **n** is perpendicular to the scattering plane spanned by the previous propagation direction and the current propagation direction, **s**.

As is displayed in Figs. 1 and 2, in the forward path, light is scattered to the **s**
_{1} direction after the first scattering event; *ϕ*
_{0} is the angle formed between the incoming perpendicular electric field component and **n**
_{1}=**s**
_{0}×**s**
_{1}/|**s**
_{0}×**s**
_{1}. After a total of *n* scattering events, light is scattered into the direction **s**
_{n} and escapes the medium. *ϕ*
_{n-1} is the angle formed between **n**
_{n-1}=**s**
_{n-2}×**s**
_{n-1}/|**s**
_{n-2}×**s**
_{n-1}| and **n**
_{n}=**s**
_{n-1}×**s**
_{n}/|**s**
_{n-1}×**s**
_{n}|. When the path is reversed, the first scattering event occurs at the *n ^{th}* scatterer, where light is scattered from

**s**

_{0}to -

**s**

_{n-1}.

*ϕ*

_{n}is the angle formed between the incoming perpendicular electric field component and

**n**′=

**s**

_{0}×(-

**s**

_{n-1})/|

**s**

_{0}×(-

**s**

_{n-1})| and

*ϕ*′

_{n}is the angle formed between

**n**′ and -

**n**

_{n-1}=(-

**s**

_{n-1})×(-

**s**

_{n-2})/|(-

**s**

_{n-1})×(-

**s**

_{n-2})|. After the

*ϕ*′

_{n}rotation, the second scattering event in the time-reversed path occurs at the

*n*-1 scatterer, followed by a rotation of

*ϕ*

_{n-2}azimuthally along -

**s**

_{n-2}from -

**n**

_{n-1}to -

**n**

_{n-2}. This continues until the last site that scatters light in the reverse path, which is the first site in the forward path. At this last scattering event, light is scattered into

**s**

_{n}from -

**s**

_{1}.

*ϕ*

_{1}is the angle between -

**n**

_{2}=(-

**s**

_{2})×(-

**s**

_{1})/|(-

**s**

_{2})×(-

**s**

_{1})| and -

**n**

_{1}=(-

**s**

_{1})×(-

**s**

_{0})/|(-

**s**

_{1})×(-

**s**

_{0})| and

*ϕ*′

_{1}is the angle between -

**n**

*1*and

**n**″=(-

**s**

_{1})×

**s**

_{n}/|(-

**s**

_{1})×

**s**

_{n}|. The angle

*ϕ*

_{1}′=0 if

**s**

_{n}=-

**s**

_{0}. After a total of

*n*scattering events the electric field in the forward path is hence given by:

and the electric field in the reverse path is given by:

Here the superscript of *S* denotes the site at which the scattering event takes place.

## 2.2 Special azimuthal rotations in the time-reversed path

Two special rotations are needed at the scatterer rer “*n*” (the first scatterer in the reverse path) in order to align the reference frame to the scattering plane spanned by -**s**
_{n-1} and -**s**
_{n-2}. Let **m**
_{0} and **n**
_{0} be the directions of the incoming parallel electric field and the incoming perpendicular electric field, respectively. In the reverse path, the first two rotations can be described by the following expression:
$\left(\begin{array}{c}{\mathbf{m}}_{0}\\ {\mathbf{n}}_{0}\\ {\mathbf{s}}_{0}\end{array}\right)\stackrel{R\left({\varphi}_{n}\right)}{\to}\stackrel{{S}^{\left(n\right)}}{\to}\left(\begin{array}{c}\mathbf{m}\prime \\ \mathbf{n}\prime \\ -{\mathbf{s}}_{n-1}\end{array}\right)\stackrel{R\left({\varphi \prime}_{n}\right)}{\to}\left(\begin{array}{c}{\mathbf{m}}_{n-1}\\ -{\mathbf{n}}_{n-1}\\ -{\mathbf{s}}_{n-1}\end{array}\right).$

The first rotation, *R*(*ϕ*
_{n}), aligns the perpendicular electric field of the incident beam to the normal of the scattering plane. The second rotation, *R*(*ϕ*′_{n}), aligns **m**′ to **m**
_{n-1}, and **n**′ to -**n**
_{n-1}. The direction, -**n**
_{n-1}, is the normal of the scattering plane for the upcoming scattering event at the second scatterer “*n*-1” in the reverse path.

At the last scattering event in the reverse path, a special rotation is also required. That is: $\left(\begin{array}{c}{\mathbf{m}}_{1}\\ {-\mathbf{n}}_{1}\\ -{\mathbf{s}}_{1}\end{array}\right)\stackrel{R\left({\varphi \prime}_{1}\right)}{\to}\left(\begin{array}{c}\mathbf{n}\u2033\times \left(-{\mathbf{s}}_{1}\right)\\ \mathbf{n}\u2033\\ -{\mathbf{s}}_{1}\end{array}\right)\stackrel{{S}^{\left(1\right)}}{\to}\left(\begin{array}{c}\mathbf{m}\u2033\\ \mathbf{n}\u2033\\ {\mathbf{s}}_{n}\end{array}\right).$

Prior to the last scattering event at the scatterer “1,” the local coordinate system is: (**m**
_{1},-**n**
_{1},-**s**
_{1})^{T}. The rotation *R*(*ϕ*′_{1}) aligns the perpendicular electric field component to the normal of the upcoming scattering plane spanned by -**s**
_{1} and **s**
_{n}. **m**″ and **n**″ are the directions of the outgoing parallel and perpendicular electric field components, respectively.

## 2.3 Relation between electric fields in the forward and time-reversed paths

Let us take the middle portions of Eqs. (3) and (4) and call them *T* and *T ^{Rev}*, respectively:

and

*T* and *T*
^{Reν} are related as will be shown. Indeed, the transpose of *T* is given by

Noting that

and the scattering amplitude matrices in the forward path and its time reversal is related by [15]:

where $Q=\left[\begin{array}{cc}1& 0\\ 0& -1\end{array}\right],$

we find:

and:

Finally, the electric field in the forward path is:

and in the reverse path:

## 3. Results and discussion

We performed EMC simulations of coherent backscattering of normally incident circularly polarized light and linearly polarized light by a turbid medium. Each simulation used 125,000 photon packets using the partial photon technique[12] and took less than twenty-four hours on a single 2.33 GHz Intel Xeon core to simulate the electric field in both the forward and reverse paths covering the backscattering directions at 73 azimuthal detection angles and 76 zenith detection angles. The turbid medium is composed of a water suspension of polystyrene spheres whose diameter is 0.1 micrometers and refractive index is 1.58984. The wavelength of the incident light is 0.5145 micrometers. The thickness of the turbid medium is 20 scattering mean free paths.

## 3.1 Intensity and enhancement of circularly polarized light around the exact backscattering direction

Figure 3 illustrates coherent backscattering of normally incident circularly polarized light. The top row displays the backscattering intensities of incident coherent light of circular polarization (left: I_{+}+I_{-}, middle: I_{+}, right: I_{-}). I_{±} is the intensity of the backscattered light of the same (or opposite) helicity as that of the incident beam. The second row displays the corresponding intensities for incident incoherent light of circular polarization. As one can see from Fig. 3, concentric circles of equal intensity are visible with or without coherence. The intensity also drops off as the zenith angle θ_{b} of detection increases in all cases. This azimuthal angle ϕ_{b} independence is expected as circularly polarized light has by definition centric symmetry.

As displayed in Fig. 4, for coherent backscattered intensities, I_{+} is greater than I_{-} close to exact backscattering direction and is smaller than I_{-} for all other angles; for incoherent backscattered intensities, I_{+} is less than I_{-} for all θ_{b} simulated. Enhancement is higher for I_{+} near the exact escape direction but falls off more sharply than I_{_} as θ_{b} increases. Incoherent backscattered light is stronger in the negative helicity channel than in the positive helicity one (I_{+}>I_{_}) because the helicity asymmetry is negative for backscattered light by small Mie scatterers in a turbid medium [16]. Multiply scattered light loses coherence and helicity simultaneously upon scattering. The enhancement factor is larger for the helicity preserved channel than the helicity flipped channel near the exact backscattering direction. This boosts the intensity of the coherent backscattered I_{+} beyond I_{_} within a narrow angular range around the exact backscattering direction.

## 3.2 Intensity and enhancement of linearly polarized light around the exact backscattering direction

Figure 5 illustrates the intensity of backscattered light for normally incident linearly polarized light around the exact backscattered direction. The first row in Fig. 5 is, from left to right, backscattered I_{x}+I_{y}, I_{x}, and I_{y} for incident coherent light polarized along the *x* direction. The second row is the same except the incident light is incoherent. The third row is the same as the first row but the zenith detection angle goes from 0 degrees to 2.25 degrees instead. Figure 6 displays the angular profiles for coherent backscattering light, incoherent backscattering light, and the enhancement factor along three directions corresponding to the azimuthal angle of 0, 45, and 90 degrees, respectively.

Like the case of circularly polarized light with and without coherence, the linearly polarized light shows symmetry about the *x*,*y* axes, but it no longer has circular symmetry. The intensity of the backscattered light of *x* polarization from the incident coherent or incoherent *x*-polarized light is, in general, greater than that of *y* polarization because the preference in linear polarization only gets lost after multiple scattering. As with circularly polarized coherent light, backscattering of linearly polarized coherent light is enhanced relative to that of incoherent light within the backscattering cone. I_{x} is enhanced greater than I_{y} because photons contributing to I_{x} preserve the original linear polarization and suffer less scattering events than photons contributing to I_{y} and hence maintain coherence better. Unlike the case involving circularly polarized light, the magnitude of the intensity of linearly polarized light, I_{x}+I_{y}, is no longer circular symmetric. This shows that there is a dependence on the azimuthal angle ϕ_{b} for backscattering of linearly polarized light. As to understand why the patterns for I_{x}’s and I_{x}+I_{y}’s are elongated along the *y*-axis and the pattern for I_{y}’s is elongated along the *x*-axis, we point out that light tends to be scattered preferably into directions out of the plane of its polarization when scattered by a Mie scatterer. Thus the *x*-component is elongated along the *y*-axis and squeezed in along the *x*-axis and the *y*-component is elongated along the *x*-axis and squeezed in along the *y*-axis. Both coherent and incoherent I_{x}’s increase slightly with the zenith angle, mostly owing to the influence of single scattered light.

Closer to the exact backscattering direction (bottom row of Fig. 5), backscattering of coherent linearly polarized light appears to display different symmetries. I_{x}+I_{y} and I_{x} appear to now be elongated along the *x*-axis and squeezed slightly at the *y*-axis due to the much stronger enhancement factor for backscattered light remitting from along the *x*-axis [see Fig. 6(f)]. Also, I_{y} displays interesting 4-fold “X” symmetry closer to the exact escape direction, which is characteristic of light being multiply scattered by Rayleigh-like particles.

## 4. Conclusion

In conclusion, the Electric field Monte Carlo method for coherent backscattering of light has been developed. This approach simulates coherent backscattering of light directly by adding coherently the electric fields of backscattered light propagating in a pair of time-reversed paths. The conventional simulation of coherent backscattering of light based on the Fourier transform of the spatial distribution of incoherent backscattering light is valid for scalar waves but introduces appreciable errors for polarized light as additional amplitude and phase differences between the two partial waves propagating in a pair of time-reversed paths appear. These errors originate from the depolarization of light. The EMC approach properly takes into full account of these additional complexities and a detailed study of this issue will be published elsewhere.

The computation of the electric field of light traveling the reverse path was simplified by using the time-reversal symmetry of electromagnetic waves that allow for the reuse of the majority of the information obtained in the forward path and as a result save a great amount of computational time. The EMC method was then used to investigate coherent backscattering of linearly or circularly polarized light from a turbid medium containing Mie scatterers. The coherent backscattering patterns were shown to depend strongly on the polarization of the incident light and detection conditions.

## Acknowledgments

This work is supported by Cottrel College Science Award.

## 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**, 38–40 (1995). [CrossRef]

**3. **S. K. Gayen and R. R. Alfano, “Emerging optical biomedical imaging techniques,” Opt. Photonics News 7, 17–22 (1996).

**4. **S. R. Arridge, “Optical tomography in medical imaging,” Inverse Probl. **15**, R41–R93 (1999). [CrossRef]

**5. **P. Wolf and G. Maret, “Weak localization and coherent backscattering of photons in disordered media,” Phys. Rev. Lett. **55**, 2696–2699 (1985). [CrossRef] [PubMed]

**6. **M. P. Van Albada and Ad Lagendijk, “Observation of weak localization of light in a random medium,” Phys. Rev. Lett. **55**, 2692–2695 (1985). [CrossRef] [PubMed]

**7. **I. Lux and L. Koblinger, *Monte Carlo Particle Transport Methods: Neutron and Photon Calculations* (CRC Press, Boca Raton, Fla., 1991).

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

**9. **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**, 400–412 (2001). [CrossRef]

**10. **B. Kaplan, G. Ledanois, and B. Villon, “Mueller matrix of dense polystyrene latex sphere suspensions: Measurements and Monte Carlo simulation,” Appl. Opt. **40**, 2769–2777 (2001). [CrossRef]

**11. **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). [CrossRef] [PubMed]

**12. **M. Xu, “Electric field Monte Carlo simulation of polarized light propagation through turbid media,” Opt. Express **12**, 6530–6539 (2004). [CrossRef] [PubMed]

**13. **K. G. Philips, M. Xu, S. K. Gayen, and R. R. Alfano, “Time-resolved ring structures of circularly polarized beams backscattered from forward scattering media,” Opt. Express **13**, 7954–7969 (2005). [CrossRef]

**14. **H. C. van de Hulst, *Light Scattering by Small Particles* (Dover, New York, 1981).

**15. **D. S. Saxon, “Tensor scattering matrix for the electromagnetic field,” Phys. Rev. **100**, 1771–1775 (1955). [CrossRef]

**16. **M. Xu and R. R. Alfano, “Circular polarization memory of light,” Phys. Rev. E 72 , 065061(R) (2005).