## Abstract

In this paper, we introduce an efficient numerical method based on surface integral equations to characterize the scattering of an arbitrarily incident Gaussian beam by arbitrarily shaped particles with multiple internal inclusions. The incident Gaussian beam is described by the Davis–Barton fifth-order approximation in combination with rotation Euler angles. For numerical purposes, the surfaces of the host particle and the inclusions are modeled using small triangular patches and the established surface integral equations are discretized with the method of moments. The resultant matrix equation is solved by using a parallel implementation of conjugate gradient method on distributed-memory architectures. Some numerical results are included to illustrate the validity and capability of the developed method. These results are also expected to provide useful insights into the scattering of Gaussian beam by composite particles.

©2012 Optical Society of America

## 1. Introduction

The scattering of light by various particles is an active and important subject of research with myriad practical applications in fields ranging from atmospheric optics to optical particle sizing to aerosol detection and Raman scattering diagnostics. With regard to the internal structure of the particles, the problem of light scattering from composite particles consisting of one or multiple inclusions within a host particle, for incidence of a plane wave or a focused Gaussian beam, has attracted much interest in recent years [1–24]. This is partially due to the fact that such particles are very common in our daily life as well as in the research for industry or environment concerns. Specific examples include raindrops, dust grains, ice crystals, daily cosmetics, medicinal sprays, biological cells and aerosols in the atmosphere.

During the years the scattering of plane wave by various composite particles has been investigated extensively by many researchers, see [1–15] to quote a few of them. For the case of an incident focused Gaussian beam, an early study was carried out by Khaled *et al*. [16]. In that paper, they applied the T-matrix method to examine the scattering of an off-axis focused Gaussian beam by a spherical particle with a concentric spherical inclusion. Later, within the framework of the generalized Lorenz-Mie theory (GLMT) [17], Gouesbet and Gréhan [18], Han *et al.* [19], Yan *et al.* [20], and Wang *et al.* [21, 22] investigate the light scattering by a spherical particle with an eccentrically located spherical inclusion illuminated by a Gaussian beam with arbitrary incidence. Subsequently, Zhang and Liao [23] employed the GLMT to study the scattering of Gaussian beam by a spherical particle with a spheroidal inclusion. Yan *et al.* [24] also studied the case of a spheroidal particle with a spherical inclusion. Despite some studies, as reviewed herein, have been carried out on the Gaussian beam scattering by several kinds of composite particles, these studies mainly focused on the cases of a sphere or a spheroid with only one internal inclusion. To the best of our knowledge, the scattering of a focused Gaussian beam by arbitrarily shaped particles with multiple internal inclusions of arbitrary shape has not been reported.

In this paper, we present an efficient numerical method for the treatment of Gaussian beam scattering by arbitrarily shaped host particles containing multiple arbitrarily shaped inclusions. Specifically, the arbitrarily incident Gaussian beams are described by using the Davis-Barton fifth-order approximate expressions [25] in combination with rotation Euler angles [26, 27] and the scattering problems involving arbitrarily shaped particles with multiple inclusions of arbitrary shape are formulated by utilizing the surface integral equations [28–32]. The solution procedure uses the well-known method of moments (MOM) [33] with Galerkin’s type testing. In implementation, the surfaces of the host particle and each inclusion are modeled by using small triangular patches, and the unknown equivalent electromagnetic currents are expanded in Rao-Wilton-Glisson (RWG) basis functions [28]. The resulting matrix equations are solved iteratively by employing a parallel implementation of conjugate gradient method (CGM) [32] on distributed-memory architectures.

This paper is organized as follows. In Section 2, a detailed description of the arbitrarily incident Gaussian beam is given. Section 3 is devoted to describe an extension of the surface integral equation method for the simulation of light scattering by arbitrarily shaped particles with multiple internal inclusions. Section 4 presents some numerical results to demonstrate the validity and capability of the present method. Finally, Section 5 concludes the paper.

## 2. Description of the incident Gaussian beam

For the purpose of the present numerical study, a detailed mathematical description of the electromagnetic field components of the incident Gaussian beam which accurately satisfies Maxwell’s equations is required. When the beam waist radius is much greater than the wavelength, the Davis first-order Gaussian beam description [34] has been found to give good results. However, for tightly focused beams, the Davis first-order Gaussian beam description satisfies Maxwell’s equations less accurately and the numerical calculations using the method described below likewise become less accurate. To accurately describe the incident Gaussian beam, particularly for tightly focused conditions, one should consider the use of higher-order approximate expressions. It is also possible to use the localized beam model developed by Gouesbet *et al* [17]. In this study, we adopt the Davis-Barton fifth-order approximation [25] in combination with rotation Euler angles [26, 27] to describe the arbitrarily incident focused Gaussian beam. Mathematically, details of the description are the following. Referring to Fig. 1
, a focused Gaussian beam with its electric field polarized in the$u$ direction at the waist is assumed to propagate along the$w$axis in its own Cartesian coordinate system${O}_{G}uvw$. The beam waist center${O}_{G}$is located at $\left({x}_{0},{y}_{0},{z}_{0}\right)$ in the Cartesian coordinate system of the particle$Oxyz$.The frame system$Oxyz$can be obtained from the beam system${O}_{G}uvw$by rotations through Euler angles$\left(\alpha ,\beta ,\gamma \right)$ followed by a translation of $\left({x}_{0},{y}_{0},{z}_{0}\right)$. To be specific, we write the rotation relationship between these two coordinate systems as

On the bases of such a transformation, the beam description of ${E}^{inc}\left({E}_{u},{E}_{v},{E}_{w}\right)$ and ${H}^{inc}\left({H}_{u},{H}_{v},{H}_{w}\right)$in the beam coordinate system${O}_{G}uvw$can be transformed to their counterparts${E}^{inc}\left({E}_{x},{E}_{y},{E}_{z}\right)$and${H}^{inc}\left({H}_{x},{H}_{y},{H}_{z}\right)$in the particle coordinate system $Oxyz$via the following formula:

For a fifth-order Davis-Barton beam [25], the electromagnetic field components in the ${O}_{G}uvw$system read as

## 3. Surface integral equation method

Now, let us consider the problem of Gaussian beam scattering by an arbitrarily shaped particle with multiple internal inclusions, as illustrated in Fig. 2 . The host particle is characterized by permittivity and permeability$\left({\epsilon}_{p},{\mu}_{p}\right)$, all the internal inclusions are assumed to be formed by homogeneous dielectric medium and the$ith\left(i=1,2,\cdots ,m\right)$ inclusion is characterized by permittivity and permeability$\left({\epsilon}_{i},{\mu}_{i}\right)$, with$m$being the total number of the inclusions. For simplicity, the surrounding medium is considered to be free space with parameters${\epsilon}_{0}$and${\mu}_{0}$. Let${S}_{p}$denote the surface of the host particle and ${S}_{i}$denote the surface of the$ith$inclusion. Also, the free space region is denoted as ${\Omega}_{0}$, the region occupied by the $ith$ internal inclusion is denoted as${\Omega}_{i}$, and the region occupied by the host particle except those occupied by all the internal inclusions is denoted as${\Omega}_{p}$. To solve this problem by means of the surface integral equation method, we introduce equivalent electromagnetic currents$\left({J}_{p},{{\rm M}}_{p}\right)$on${S}_{p}$ and$\left({J}_{i},{{\rm M}}_{i}\right)$on${S}_{i}\left(i=1,2,\cdots ,m\right)$. Through the Stratton-Chu representation formulas, the scattered fields in each region can be expressed in terms of the equivalent electromagnetic currents. Specifically, the scattered fields in region${\Omega}_{0}$, due to the equivalent electric and magnetic currents ${J}_{p}$ and ${{\rm M}}_{p}$on${S}_{p}$, can be expressed as

By enforcing boundary condition on the surface of the host particle, namely that the total tangential fields across the interface are continuous, we can obtain two boundary integral equations as follows:

At this point, the integral equations established above can be solved numerically by using the MOM [33]. In implementation, the surfaces of the host particle and the inclusions are modeled by using small triangular patches, and the unknown equivalent electromagnetic currents are expanded in vector basis functions associated with the edges of the triangular patches. For triangular patches, the most suitable basis functions are the so-called RWG basis functions developed by Rao, Wilton and Glisson [28]. Mathematical representation of the RWG basis function associated with the$nth$edge is given as

Therefore, Eq. (31) can take the form

The solution to Eq. (33) could be obtained by a direct method or an iterative matrix solver. However, traditional methods incur very high computational cost and memory requirements with the increasing of the unknowns. In addition, the conventional approaches to computing the matrix elements consume a considerable portion of the total solution time. One approach to overcome these difficulties is to adopt the parallel computation. In this paper, we use a parallel implementation of the CGM to solve the resultant matrix equation. The parallel program is developed on distributed-memory architectures using message passing interface. The coefficient matrix is decomposed by the row and stored in distributed memory of the processor. The detailed description of the parallel CGM is given in the authors’ previous paper [32] and is not repeated here. Once the unknown coefficients for the equivalent electric and magnetic currents are solved from Eq. (33), the far-zone scattered fields ${E}_{far}^{sca}$ and ${H}_{far}^{sca}$ radiated by the equivalent currents ${J}_{p}$ and${{\rm M}}_{p}$can be obtained by

## 4. Numerical results and discussion

In this section, we present some of our preliminary numerical results for the scattering of Gaussian beam by arbitrarily shaped particles with multiple internal inclusions. The results to be presented are in terms of the differential scattering cross section (DSCS) in the particle system$Oxyz$, the E-plane corresponds to the$xOz$-plane and the H-plane corresponds to the $yOz$-plane. The DSCS is defined as

Without loss of generality, the amplitude of the electric field for the incident beam is assumed to be unity. In what follows, all the computations are performed on a HP Workstation. It has 16 nodes and each node is Intel Xeon MP 3.0 GHz, 4 GB memory. Iterations are carried out until the residual error is reduced to below 0.001.

In order to demonstrate the validity of the present method for the case of Gaussian beam illumination, we first consider a spherical particle with a concentric spherical inclusion. The radii of the host sphere and the inclusion are$r=1.0\lambda $and ${r}_{1}=0.5\lambda $with refractive indices of $m=1.33$and${m}_{1}=1.55$, respectively, with$\lambda $being the wavelength of the incident beam. The particle is illuminated by an obliquely incident Gaussian beam with ${\omega}_{0}=1.0\lambda $. The beam waist center is located at the origin of the particle system with ${x}_{0}={y}_{0}={z}_{0}=0$ and the rotation Euler angles are specified as $\alpha ={0}^{o},\beta ={45}^{o}$and$\gamma ={0}^{o}$. For numerical solution, the surfaces of the host particle and the inclusion are respectively discretized into 3788 and 1024 triangular patches using the mesh density of 12 parts per wavelength. As a result, a total of 14436 unknowns are generated. Performing the parallel CGM on the HP Workstation with 40 processors, 112 s is needed to solve the resultant matrix equation. Figure 4 presents the computed DSCSs as a function of the scattering angle in both the E-plane and the H-plane. The results obtained from the GLMT are also presented. As is evident from the figure, good agreement is obtained between the results obtained from the two methods.

Now, we consider the problem of Gaussian beam scattering by a spheroidal particle with two identical spherical inclusions, as illustrated in Fig. 5 . The semimajor axis and the semiminor axis of the host spheroid are$a=2.0\lambda $and$b=1.0\lambda $, respectively, and the refractive index is$m=1.414$. Both inclusions are located on the$z$-axis in the particle system with${z}_{1}=-1.0\lambda $and${z}_{2}=1.0\lambda $, respectively. Each spherical inclusion has a radius of${r}_{i}=0.5\lambda $and the refractive index of${m}_{i}=2.0$. The beam waist is centered at $\left({x}_{0},{y}_{0},{z}_{0}\right)=\left(0.0,0.0,0.0\right)$, and the angle set of the beam is$\alpha =\beta =\gamma ={0}^{o}$. The computed DSCSs for different values of the beam waist radius are given in Fig. 6 . For comparison purpose, the results for the case of plane wave illumination obtained with the commercial soft FEKO are given in the same figure. It is obvious that the DSCS for Gaussian beams is smaller than that for a plane wave because of the influence of the beam shape coefficients. As expected, the results in the case of Gaussian beam incidence with a relatively large waist radius of${\omega}_{0}=20.0\lambda $are in excellent agreement with the results in the case of plane wave incidence. This further demonstrates the validity the present method.

Next, we present three examples to illustrate the capabilities of the present method. The first one is a spherical particle with four inclusions of different shape, as depicted in Fig. 7 . The radius of the host sphere is$r=2.0\lambda $and the refractive index is$m=1.5$. These four internal inclusions are assumed to be formed by identical homogeneous medium characterized by refractive index${m}_{i}=1.8$. The first inclusion is a sphere with radius${r}_{1}=0.5\lambda $, the second one is a cylinder with radius${r}_{2}=0.25\lambda $and height${h}_{2}=1.0\lambda $, the third one is an oblate spheroid with semi-axes${a}_{3}=0.5\lambda $and${b}_{3}=0.4\lambda $, and the fourth inclusion is a cone with radius${r}_{4}=0.5\lambda $and height${h}_{4}=1.0\lambda $. The positions of the internal inclusions with respect to the global coordinate system of the host particle is specified by the coordinates $\left({x}_{1},{y}_{1},{z}_{1}\right)=\left(-1.0,0.0,0.0\right)\lambda $,$\left({x}_{2},{y}_{2},{z}_{2}\right)=\left(0.0,0.0,-1.0\right)\lambda $, $\left({x}_{3},{y}_{3},{z}_{3}\right)=\left(1.0,0.0,0.0\right)\lambda $and $\left({x}_{4},{y}_{4},{z}_{4}\right)=\left(0.0,0.0,0.0\right)\lambda $. The particle is illuminated by a focused Gaussian beam at incidence angles$\alpha =\beta =\gamma ={0}^{o}$with beam waist${\omega}_{0}=1.5\lambda $. The location of the beam waist center is a varied parameter. Figure 8 shows the effects of the beam waist center positioning on the DSCS distributions for both the E and H-planes. It is found that the position offset of the beam waist center makes a small difference to the DSCS by reducing it for most but not all angles. This occurs because the beam waist center deviate the center of the particle leading to the change of the distribution of the scattering intensity.

The second example considered is the scattering of a focused Gaussian beam by a cubic particle containing 27 randomly distributed spherical inclusions, as shown in Fig. 9 . All the inclusions are assumed to be uniform and the positions are generated from the Monte Carlo method [39] with fractional volume$f=6.0\%$. The side length of the host cube is $l=3.0\lambda $ and the refractive index is$m=1.2-i0.2$. Each spherical inclusion has a radius ${r}_{i}=0.25\lambda $ and the refractive index$m=1.5-i0.1$. Both the centers of the host cube and the beam waist are located at the origin of the particle system with${x}_{0}={y}_{0}={z}_{0}=0$. The incident focused Gaussian beam has a beam waist of ${\omega}_{0}=2.0\lambda $and the incident angle is taken in the E-plane, that is$\alpha =\gamma ={0}^{o}$. Thus, the E-plane scattering property is what we are most interested in. Figure 10 presents angular distributions of the DSCS in the E-plane for the Gaussian beam with different values of incident angle$\beta $. It can be observed that the scattering angle corresponding to the extreme value of the DSCS is just coincident with the incident direction of the beam, which indicates that most of the scattering energy concentrates forward. This is identical with the general idea of light scattering theory.

In the last example, we apply the present method to simulate the scattering of an off-axis obliquely incident Gaussian beam by a hexagonal prism containing a fractal aggregate that consists of 100 identical spherical particles, as illustrated in Fig. 11 . Internal aggregate positions were generated by using the cluster-cluster aggregation (CCA) algorithm described by Mackowski [40], with prefactor constant${k}_{f}=1.19$and fractal dimension${D}_{f}=1.82$. The radius and the height of the host hexagonal prism are$r=2.0\lambda $and$h=4.0\lambda $, respectively, and the refractive index is$m=1.313$. Each primary particle of the internal fractal aggregate has radius${r}_{i}=0.2\lambda $and refractive index$m=1.50-i0.65$. The beam waist is centered at$\left({x}_{0},{y}_{0},{z}_{0}\right)=\left(-1.0,-1.0,-1.0\right)\lambda $ with a beam waist radius of ${\omega}_{0}=1.5\lambda $, and the rotation Euler angles are specified as $\alpha ={0}^{o},\beta ={45}^{o}$and$\gamma ={0}^{o}$. Figure 12 presents the simulated DSCSs as a function of the scattering angle in both the E-plane and the H-plane.

## 5. Conclusion

In this paper, we present an efficient numerical method for the simulation of light scattering by arbitrarily shaped particles with multiple internal inclusions illuminated by a focused Gaussian beam. Specifically, the Davis-Barton fifth-order approximate expressions in combination with rotation Euler angles is employed to describe the arbitrarily incident Gaussian beam. The surface integral equations are applied to formulate the scattering problems involving arbitrarily shaped particles with multiple internal inclusions and are numerically discretized by the MoM with the RWG basis functions defined on triangular patches. To reduce computational burden, the resultant matrix equation is solved by a parallel implementation of the CGM on distributed-memory architectures. The present method is validated and its capability illustrated in several characteristic examples. This work is helpful for further research on the scattering of an arbitrarily incident focused Gaussian beam by composite particles.

## References and links

**1. **F. Borghese, P. Denti, R. Saija, and O. I. Sindoni, “Optical properties of spheres containing a spherical eccentric inclusion,” J. Opt. Soc. Am. A **9**(8), 1327–1335 (1992). [CrossRef]

**2. **N. C. Skaropoulos, M. P. Ioannidou, and D. P. Chrissoulidis, “Indirect mode-matching solution to scattering from a dielectric sphere with an eccentric inclusion,” J. Opt. Soc. Am. A **11**(6), 1859–1866 (1994). [CrossRef]

**3. **F. Borghese, P. Denti, and R. Saija, “Optical properties of spheres containing several spherical inclusions,” Appl. Opt. **33**(3), 484–493 (1994). [CrossRef] [PubMed]

**4. **G. Videen, D. Ngo, P. Chylek, and R. G. Pinnick, “Light scattering from a sphere with an irregular inclusion,” J. Opt. Soc. Am. A **12**(5), 922–928 (1995). [CrossRef]

**5. **D. Ngo, G. Videen, and P. Chýlek, “A FORTRAN code for the scattering of EM waves by a sphere with a nonconcentric spherical inclusion,” Comput. Phys. Commun. **99**(1), 94–112 (1996). [CrossRef]

**6. **A. Macke, M. I. Mishchenko, and B. Cairns, “The influence of inclusions on light scattering by large ice particles,” J. Geophys. Res. **101**(D18), 23311–23316 (1996). [CrossRef]

**7. **M. I. Mishchenko, J. W. Hovenier, and L. D. Travis, *Light Scattering by Nonspherical Particles:Ttheory, Measurements, and Applications* (Academic, San Diego, 2000).

**8. **D. R. Prabhu, M. Davies, and G. Videen, “Light scattering calculations from oleic-acid droplets with water inclusions,” Opt. Express **8**(6), 308–313 (2001). [CrossRef] [PubMed]

**9. **M. P. Ioannidou and D. P. Chrissoulidis, “Electromagnetic-wave scattering by a sphere with multiple spherical inclusions,” J. Opt. Soc. Am. A **19**(3), 505–512 (2002). [CrossRef] [PubMed]

**10. **T. Weigel, J. Schulte, and G. Schweiger, “Inelastic scattering on particles with inclusions,” J. Opt. Soc. Am. A **22**(6), 1048–1052 (2005). [CrossRef] [PubMed]

**11. **A. Doicu, T. Wriedt, and Y. A. Eremin, *Light Scattering by Systems of Particles* (Springer, Berlin, 2006).

**12. **A. P. Moneda and D. P. Chrissoulidis, “Dyadic Green’s function of a sphere with an eccentric spherical inclusion,” J. Opt. Soc. Am. A **24**(6), 1695–1703 (2007). [CrossRef] [PubMed]

**13. **D. K. Wu and Y. P. Zhou, “Forward scattering light of droplets containing different size inclusions,” Appl. Opt. **48**(15), 2957–2965 (2009). [CrossRef] [PubMed]

**14. **M. Mikrenska and P. Koulev, “Simulation of light scattering by large particles with randomly distributed spherical or cubic inclusions,” J. Quant. Spectrosc. Radiat. Transf. **110**(14–16), 1411–1417 (2009). [CrossRef]

**15. **S. Xian-Ming, W. Hai-Hua, L. Wan-Qiang, and S. Jin, “Light scattering by a spherical particle with multiple densely packed inclusions,” Chin. Phys. B **18**(3), 1040–1044 (2009). [CrossRef]

**16. **E. E. Khaled, S. C. Hill, and P. W. Barber, “Light scattering by a coated sphere illuminated with a Gaussian beam,” Appl. Opt. **33**(15), 3308–3314 (1994). [CrossRef] [PubMed]

**17. **G. Gouesbet and G. Gréhan, *Generalized Lorenz*-*Mie Theories* (Springer, Berlin, 2011).

**18. **G. Gouesbet and G. Gréhan, “Generalized Lorenz-Mie theory for a sphere with an eccentrically located spherical inclusion,” J. Mod. Opt. **47**(5), 821–837 (2000).

**19. **G. X. Han, Y. P. Han, J. Y. Liu, and Y. Zhang, “Scattering of an eccentric sphere arbitrarily located in a shaped beam,” J. Opt. Soc. Am. B **25**(12), 2064–2072 (2008). [CrossRef]

**20. **B. Yan, X. E. Han, and K. F. Ren, “Scattering of a shaped beam by a spherical particle with an eccentric spherical inclusion,” J. Opt. A **11**(1), 015705 (2009). [CrossRef]

**21. **J. J. Wang, G. Gouesbet, Y. P. Han, and G. Gréhan, “Study of scattering from a sphere with an eccentrically located spherical inclusion by generalized Lorenz-Mie theory: internal and external field distribution,” J. Opt. Soc. Am. A **28**(1), 24–39 (2011). [CrossRef] [PubMed]

**22. **J. J. Wang, G. Gouesbet, G. Gréhan, Y. P. Han, and S. Saengkaew, “Morphology-dependent resonances in an eccentrically layered sphere illuminated by a tightly focused off-axis Gaussian beam: parallel and perpendicular beam incidence,” J. Opt. Soc. Am. A **28**(9), 1849–1859 (2011). [CrossRef]

**23. **H. Y. Zhang and T. Q. Liao, “Scattering of Gaussian beam by a spherical particle with a spheroidal inclusion,” J. Quant. Spectrosc. Radiat. Transf. **112**(9), 1486–1491 (2011). [CrossRef]

**24. **B. Yan, H. Y. Zhang, and C. H. Liu, “Scattering of Gaussian beam by a spheroidal particle with a spherical inclusion at the center,” Opt. Commun. **284**(16–17), 3811–3815 (2011). [CrossRef]

**25. **J. P. Barton and D. R. Alexander, “Fifth-order corrected electromagnetic field components for a fundamental Gaussian beam,” J. Appl. Phys. **66**(7), 2800–2802 (1989). [CrossRef]

**26. **A. R. Edmonds, *Angular Momentum in Quantum Mechanics* (Princeton University Press, Princeton, 1957).

**27. **G. Gouesbet, J. J. Wang, and Y. P. Han, “Transformations of spherical beam shape coefficients in generalized Lorenz-Mie theories through rotations of coordinate systemsIII. Special Euler angles,” Opt. Commun. **283**(17), 3235–3243 (2010). [CrossRef]

**28. **S. M. Rao, D. R. Wilton, and A. W. Glisson, “Electromagnetic scattering by surfaces of arbitrary shape,” IEEE Trans. Antenn. Propag. **30**(3), 409–418 (1982). [CrossRef]

**29. **S. M. Rao, C. C. Cha, R. L. Cravey, and D. Wilkes, “Electromagnetic scattering from arbitrary shaped conducting bodies coated with lossy materials of arbitrary thickness,” IEEE Trans. Antenn. Propag. **39**(5), 627–631 (1991). [CrossRef]

**30. **P. Ylä-Oijala, M. Taskinen, and J. Sarvas, “Surface integral equation method for general composite metallic and dielectric structures with junctions,” Prog. Electromagn. Res. **52**, 81–108 (2005). [CrossRef]

**31. **J. Rivero, J. M. Taboada, L. Landesa, F. Obelleiro, and I. García-Tuñón, “Surface integral equation formulation for the analysis of left-handed metamaterials,” Opt. Express **18**(15), 15876–15886 (2010). [CrossRef] [PubMed]

**32. **Z. W. Cui, Y. P. Han, and M. L. Li, “Solution of CFIE-JMCFIE using parallel MOM for scattering by dielectrically coated conducting bodies,” J. Electromagn. Waves Appl. **25**(2), 211–222 (2011). [CrossRef]

**33. **R. F. Harrington, *Field Computation by Moment Methods* (Macmillan, New York, 1968).

**34. **L. W. Davis, “Theory of electromagnetic beams,” Phys. Rev. A **19**(3), 1177–1179 (1979). [CrossRef]

**35. **Z. W. Cui, Y. P. Han, and H. Y. Zhang, “Scattering of an arbitrarily incident focused Gaussian beam by arbitrarily shaped dielectric particles,” J. Opt. Soc. Am. B **28**(11), 2625–2632 (2011). [CrossRef]

**36. **Z. W. Cui, Y. P. Han, and Q. Xu, “Numerical simulation of multiple scattering by random discrete particles illuminated by Gaussian beams,” J. Opt. Soc. Am. A **28**(11), 2200–2208 (2011). [CrossRef] [PubMed]

**37. **D. A. Dunavant, “High degree efficient symmetrical Gaussian quadrature rules for the triangle,” Int. J. Numer. Methods Eng. **21**(6), 1129–1148 (1985). [CrossRef]

**38. **R. D. Graglia, “On the numerical integration of the linear shape functions times the 3-D green’s function or its gradient on a plane triangle,” IEEE Trans. Antenn. Propag. **41**(10), 1448–1455 (1993). [CrossRef]

**39. **L. Tsang, J. A. Kong, K. H. Ding, and C. O. Ao, *Scattering of Electromagnetic Waves: Numerical Simulations* (Wiley, 2001).

**40. **D. W. Mackowski, “A simplified model to predict the effects of aggregation on the absorption properties of soot particles,” J. Quant. Spectrosc. Radiat. Transf. **100**(1–3), 237–249 (2006). [CrossRef]