## Abstract

We report what we believe to be the first simulation of enhanced backscattering (EBS) of light by numerically solving Maxwell’s equations without heuristic approximations. Our simulation employs the pseudospectral time-domain (PSTD) technique, which we have previously shown enables essentially exact numerical solutions of Maxwell’s equations for light scattering by millimeter-volume random media consisting of micrometer-scale inhomogeneities. We show calculations of EBS peaks of random media in the presence of speckle; in addition, we demonstrate speckle reduction using a frequency-averaging technique. More generally, this new technique is sufficiently robust to permit the study of EBS phenomena for random media of arbitrary geometry not amenable to simulation by other approaches, especially with regard to extension to full-vector electrodynamics in three dimensions.

©2005 Optical Society of America

## 1. Introduction

Recently, the phenomenon of enhanced backscattering (EBS) of light, also known as coherent backscattering (CBS), has attracted significant attention in the study of light propagation in random media [1–5]. This phenomenon is due to the interference effect of coherent light scattered by random media, where light interferes between time-reversed light paths, resulting in an enhanced intensity cone in the backscattering direction.

Conventionally, the EBS peak can be indirectly simulated using Monte Carlo models [6, 7] by statistically treating the problem of light interaction within random media as series of independent scattering events. By approximating light propagation as a diffusion problem, the photon radial probability distribution is obtained, which can then be used to yield the EBS peak indirectly by means of a Fourier transform. We note that Monte Carlo simulations are often performed beyond the diffusion approximation, and thus could have applicability to the vectorial case.

An alternative technique presented in [8] investigated the multiple scattering of scalar waves in diffusive media by means of the radiative transfer equation. This approach does not rely on the diffusion approximation, and becomes asymptotically exact when the scattering mean free path is much larger than the wavelength. An important element discussed in detail in [8] was the solution of the Milne equation, which provides the coherent backscattering cone without relying on the diffusion approximation.

A worthwhile goal is the development of predictive techniques for full-vector optical wave interactions with potentially dense, three-dimensional (3-D) random media. Here, the constituent scatterers can be inhomogeneous, arbitrarily shaped, strongly coupled in the near field, and span size scales ranging from subwavelength to many wavelengths. Existing analytical and numerical techniques may not be sufficiently robust to deal with this level of complexity [9–11]. Thus, a rigorous simulation technique based on first-principles electromagnetic theory (i.e., the full set of Maxwell’s equations) is preferred.

In pursuit of this goal, we report in this paper what we believe to be the first direct, numerical simulation of the EBS phenomenon by solving Maxwell’s equations for millimeter-volume random media consisting of micrometer-scale inhomogeneities. We employ the pseudospectral time-domain (PSTD) technique [12–14]. Furthermore, we present a frequency-averaging method that can significantly suppress speckles in EBS simulations.

With the numerical methods reported in this paper, EBS phenomena can be simulated without heuristic approximations for two-dimensional random media having constituent scatterers of arbitrary shape, composition, size, and spacing. While the specific results reported in this paper are for two-dimensional problems, it should be understood that our method is directly extendible and applicable to full-vector solutions of Maxwell’s equations in three dimensions. In fact, we have developed computer programs for this purpose, and will report fully three-dimensional results in subsequent papers.

## 2. Methods

With the PSTD technique, the spatial derivatives of the electric and magnetic fields that are required to implement Maxwell’s curl equations are obtained using the differentiation theorem for Fourier transforms:

where **F** and **F**
^{-1} denote, respectively, the forward and inverse discrete Fourier transforms, and *k*̃_{x} is the Fourier transform variable representing the *x*-component of the numerical wavevector. The spatial derivatives {(*∂V/∂x*)_{i}} can be calculated in one step. In multiple dimensions, this process is repeated for each cut parallel to the major axes of the space lattice. PSTD techniques have been shown to possess spectral accuracy; that is, errors due to spatial sampling decrease exponentially as the meshing density increases beyond the Nyquist rate. With trigonometric basis functions, this permits the PSTD meshing density to approach two samples per wavelength in each spatial dimension.

Briefly summarizing key details of our algorithm, we eliminate the wraparound caused by the periodicity in the discrete Fourier transforms used to implement PSTD by using the anisotropic perfectly matched layer (APML) absorbing boundary condition [15]. Incident wave excitation is provided by the scattered-field technique [16]. The surfaces of the scattering shapes are approximated by staircasing at the grid resolution. In this manner, we have demonstrated PSTD results for plane-wave scattering by canonical shapes that are accurate to better than 1 dB over dynamic ranges exceeding 50 dB. This holds for the scattering intensity observed at all possible angles, including backscattering.

By employing PSTD, the EBS phenomenon in random media can be solved without heuristic approximations, such as treating light as a diffusion problem or treating the random media as a cluster of point scatterers. Further, in PSTD simulations, the scattered light intensity for a broadband spectrum of wavelengths (λ_{0} ranging from 1 µm up to 600 µm) can be obtained in a single simulation.

## 3. PSTD Simulation of EBS

We report the initial application of PSTD to model two-dimensional (2-D) transverse-magnetic EBS scattering of light by an 800×400 µm rectangular cluster of scatterers. As depicted in Fig. 1, the cluster consists of *N* randomly positioned, closely packed, infinitely long, dielectric cylinders with refractive index n=1.25. The cluster geometry is created by randomly positioning scatterers in free space, subject to the constraint of a minimum distance of 2.4 µm allowed between cylinders.

For the incident wavelength λ_{0}=1 µm, the transport mean free path, *l _{s}*’, is 37.74 µm and 5.59 µm, for

*N*=10,000 and

*N*=20,000, respectively. We use a PSTD grid having a uniform spatial resolution of 0.33 µm, equivalent to 0.42

*λ*

_{d}(

*λ*

_{d}: optical wavelength in dielectric material) at frequency 300 THz for a cylinder refractive index

*n*=1.25. In each simulation, the cluster is illuminated by a coherent plane wave at a 15° incident angle to avoid specular reflection. Both the incident light and backscattered light are polarized perpendicular to the plane of incidence, equivalent to collinear detection in EBS experiments.

Due to the limitations of finite computer memory, we necessarily can model only a finite random medium region. However, we have determined that we can efficiently account for the finite size of our modeling region via convolution of comparative benchmark analytical results with an appropriate windowing function (representing the finite aperture of the modeled region).

Figure 2(a) illustrates sample data for the scattered light intensity calculated in a single PSTD simulation. Here, 0° corresponds to direct backscattering at 15° from the normal. Since the incident light source is coherent, the scattered light intensity contains considerable speckle. In order to suppress the speckle, the scattered light intensity is ensemble-averaged over several simulations, each corresponding to a different random arrangement of cylinders within the rectangular cluster. As shown in Fig. 2(b), after averaging over 40 PSTD simulations, the speckle contribution is significantly suppressed; and the formerly obscured EBS peak becomes visible.

Figure 2(c) illustrates how frequency-averaging can further suppress speckle. Due to coherent interference effects, speckle occurs at different angular locations for different wavelengths. Here, speckle is significantly reduced by averaging over 50 PSTD-computed scattered intensities corresponding to slightly different incident frequencies within a small frequency range centered at *f*
_{0}. This is similar to experimental observations of EBS using non-monochromatic illumination with a temporal coherence length of 10 µm. (The 50 PSTD-computed scattered intensities correspond to 50 frequencies evenly spaced between 1.05**f*
_{0} and 0.95**f*
_{0}. An estimated 16 modes are averaged in the process to suppress speckle.) As shown in Fig. 2(c), following frequency-averaging, the EBS peak is clearly visible, showing a significant correlation with experimental results reported by Tomita and Ikari [17]. We believe this to be the first report of direct observation of the EBS peak from numerical experiments based upon Maxwell’s equations, without employing indirect transformations or heuristic approximations.

Figure 3 compares the angular distribution of the PSTD-computed EBS peak with that obtained using standard EBS theory based on the diffusion approximation [18]:

where *l*=transport mean free path length; *q*=2*πθ/λ* ; z_{0}=2*l* 3; * denotes a convolution; and *SF(θ)* is the far-field scattering function for a homogeneous slab of the same size illuminated by a plane wave. We observe a good agreement between the PSTD simulations and the benchmark EBS theory.

## 6. Summary and discussion

We have reported a new PSTD simulation technique for the EBS phenomenon which permits numerical solution of Maxwell’s equations for millimeter-volume random media consisting of micrometer-scale inhomogeneities of arbitrary shape, size, spacing, and composition. In addition, we have shown that speckle arising from these simulations can be reduced significantly by a frequency-averaging method. We believe this to be the first direct numerical simulation of the EBS phenomenon from Maxwell’s equations without employing indirect transformations or heuristic approximations.

While the specific results reported in this paper are for two-dimensional problems that potentially can be solved using existing approaches, it should be understood that our method is directly extendible and applicable to full-vector solutions of Maxwell’s equations in three dimensions. In fact, we have developed computer programs for this purpose, and will report fully three-dimensional results in subsequent papers utilizing the resources currently available on the NSF TeraGrid. A broad implication of this work is the emerging feasibility of computational biophotonics based upon the full-vector Maxwell’s equations in three dimensions.

## Acknowledgments

The authors thank Hariharan Subramanian for insightful comments. In addition, the authors thank the NIH National Cancer Institute Contract Grants 5R01-CA085991 and 5R01-HD044015, NSF Grant BES-0238903, and NSF TeraGrid Grant No. MCB040062N for their support of this research. This work was performed in part under the auspices of the U.S. DOE by UC, LLNL under Contract W-7405-Eng-48. Young L. Kim was supported in part by the cancer research and prevention foundation. Snow H. Tseng’s email address is snow@ece.northwestern.edu.

## References and links

**1. **E. Akkermans and G. Montambaux, “Mesoscopic physics of photons” J. Opt. Soc. Am. B **21**, 101–112 (2004). [CrossRef]

**2. **Y.L. Kim, Y. Liu, V.M. Turzhitsky, H.K. Roy, R.K. Wali, and V. Backman, “Coherent backscattering spectroscopy” Opt. Lett. **29**, 1906–1908 (2004). [CrossRef] [PubMed]

**3. **R. Sapienza, S. Mujumdar, C. Cheung, A.G. Yodh, and D. Wiersma, “Anisotropic weak localization of light” Phys. Rev. Lett.92 (2004). [CrossRef] [PubMed]

**4. **G. Labeyrie, D. Delande, C.A. Muller, C. Miniatura, and R. Kaiser, “Coherent backscattering of light by cold atoms: Theory meets experiment” Europhys. Lett. **61**, 327–333 (2003). [CrossRef]

**5. **I.V. Meglinski, V.L. Kuzmin, D.Y. Churmakov, and D.A. Greenhalgh, “Monte Carlo simulation of coherent effects in multiple scattering” Proc. of the Royal Society of London Series A-Mathematical Physical and Engineering Sciences , **461**: 43–53 (2005). [CrossRef]

**6. **R. Lenke, R. Tweer, and G. Maret, “Coherent backscattering of turbid samples containing large Mie spheres” J. Opt. A **4**, 293–298 (2002). [CrossRef]

**7. **V.L. Kuzmin and I.V. Meglinski, “Coherent multiple scattering effects and Monte Carlo method” JETP Lett. **79**, 109–112 (2004). [CrossRef]

**8. **E. Amic, J.M. Luck, and T.M. Nieuwenhuizen, “Anisotropic multiple scattering in diffusive media” J. Phys. A: Math. Gen. **29**, 4915–4955 (1996). [CrossRef]

**9. **J. Pearce and D.M. Mittleman, “Propagation of single-cycle terahertz pulses in random media” Opt. Lett. **26**, 2002–2004 (2001). [CrossRef]

**10. **M. Haney and R. Snieder, “Breakdown of wave diffusion in 2D due to loops” Phys. Rev. Lett.91, (2003). [CrossRef] [PubMed]

**11. **L. Marti-Lopez, J. Bouza-Dominguez, J.C. Hebden, S.R. Arridge, and R.A. Martinez-Celorio, “Validity conditions for the radiative transfer equation” J. Opt. Soc. Am. A **20**, 2046–2056 (2003). [CrossRef]

**12. **Q.H. Liu, “Large-scale simulations of electromagnetic and acoustic measurements using the pseudospectral time-domain (PSTD) algorithm” IEEE Transactions on Geoscience and Remote Sensing , **37**, 917–926 (1999). [CrossRef]

**13. **Q.H. Liu, “The PSTD algorithm: A time-domain method requiring only two cells per wavelength” Microwave Opt. Technol. Lett. **15**, 158–165 (1997). [CrossRef]

**14. **S.H. Tseng, J.H. Greene, A. Taflove, D. Maitland, V. Backman, and J.T. Walsh, “Exact solution of Maxwell’s equations for optical interactions with a macroscopic random medium” Opt. Lett. **29**, 1393–1395 (2004); **30**: 56–57 (2005). [CrossRef] [PubMed]

**15. **S.D. Gedney, “An anisotropic perfectly matched layer absorbing media for the truncation of FDTD lattices” IEEE transactions on Antennas and Propagation , **44**: 1630–1639 (1996). [CrossRef]

**16. **A. Taflove and S.C. Hagness, *Computational Electrodynamics: the finite-difference time-domain method*. 2000: Artech House. 852.

**17. **M. Tomita and H. Ikari, “Influence of Finite Coherence Length of Incoming Light on Enhanced Backscattering” Physical Review B , **43**: 3716–3719 (1991). [CrossRef]

**18. **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**, 1471–1474 (1986). [CrossRef] [PubMed]