## Abstract

The performance of most widespread surface integral equation (SIE) formulations with the method of moments (MoM) are studied in the context of plasmonic materials. Although not yet widespread in optics, SIE-MoM approaches bring important advantages for the rigorous analysis of penetrable plasmonic bodies. Criteria such as accuracy in near and far field calculations, iterative convergence and reliability are addressed to assess the suitability of these formulations in the field of plasmonics.

©2012 Optical Society of America

## 1. Introduction

The proper characterization of the plasmonic optical response of metals constitutes an outstanding task in nanoscience and nanotechnology. The field confinement properties of collective electron excitations at metal interfaces, surpassing the diffraction limit, give rise to interesting applications in waveguiding or sensing, among others [1–3]. The most important metals for plasmonics in visible and near-infrared bands are silver and gold. The main attraction for silver is its low losses, while for gold is its chemical stability. Nevertheless, other metals such as copper or aluminum are also considered, the last with a good plasmon resonance in the ultraviolet. When dealing with metallic nanoparticles (sub-wavelength scale) the excitation of surface plasmons by light is denoted as localized surface plasmon resonance (LSPR). Light intensity enhancement is a very important characteristic of LSPRs and therefore, a key aspect in nanoantennas analysis and design at optical frequencies [4–6].

The interest of getting control over light propagation has also given a boost to the design and engineering of photonic crystals. Increasing the refractive index contrast by means of metallic features in photonic crystals looking for enhancing or modifying the performance of conventional optical components and nanostructures is an interesting and active research field [7–9].

The wide range of applications and exploitation possibilities of optical plasmonic properties of metals demand numerical techniques which provide accurate modeling and analysis of problems involving metal-dielectric interfaces in visible and near-infrared bands. Since optical response of metals can be well-described in the classical framework based on Maxwell’s equations [1], differential formulations as the finite difference time domain (FDTD) [10] or the finite integration technique [11] have been commonly used to solve this kind of problems due to their easy implementation from differential equations. Mainly due to the high computational requirements of the differential techniques, alternative tools demanding a lower number of unknowns for a given problem such as surface integral equation (SIE) formulations solved by the well-known method of moments (MoM) [12] have increased their presence in the context of plasmonics [13–15].

A study of the performance of most widespread SIE formulations was carried out in [16] for metallic, dielectric, piecewise dielectric and composite objects and later on in [17] for dielectric objects in the framework of MoM acceleration techniques. A similar study involving SIE formulations applied to left-handed metamaterials (LHM’s) has been accomplished recently in [18] by the authors, and later on in [19] using acceleration algorithms. In this work, the comparison of these formulations properties is established in the interesting context of plasmonic materials.

The paper is structured as follows: the SIE formulations under consideration are presented in Section 2. A complete set of near and far field results obtained with these formulations are contrasted with analytical references in Section3, in order to assess their accuracy and their iterative performance. Finally, the concluding remarks are gathered in Section 4.

## 2. Surface integral equation formulations

Let us consider the combination of normal and tangential equations derived from the boundary conditions imposed separately to the electric and magnetic fields in order to derive a set of stable and well-tested SIEs for (piecewise) homogeneous real conductors and/or dielectric objects. The following proposal of combination of the tangential electric field integral equation (T-EFIE), the tangential magnetic field integral equation (T-MFIE), the normal electric field integral equation (N-EFIE) and the normal magnetic field integral equation (N-MFIE) has proven to be stable [16]:

*η*is the intrinsic impedance in region

_{i}*R*(

_{i}*R*

_{1}and

*R*

_{2}are the exterior and the interior regions of the material, respectively) and

*a*,

_{i}*b*,

_{i}*c*and

_{i}*d*are the coefficients to handle the different formulations (see Table 1 ). Tangential-EFIE/MFIE and normal-EFIE/MFIE are derived from the equivalence theorem and the boundary conditions for the electric/magnetic fields tangential to the boundary surfaces. The tangential equations are formulated according to $-\widehat{n}\times \widehat{n}\times F\equiv {F}_{tan}$, where $\widehat{n}$ is the vector normal to the objects boundary surfaces and interfaces,

_{i}**F**is the electric/magnetic field, and subindex

*tan*denotes tangential component. Alternatively, in the normal equations the boundary conditions are imposed for the rotated tangential components of the fields according to $\widehat{n}\times F$. The details of the general SIE-MoM formulation for the analysis of multiple plasmonic media can be looked up in [15].

The comparative study presented here involves different known formulations that can be obtained by selecting the complex combination parameters *a _{i}*,

*b*,

_{i}*c*and

_{i}*d*according to Table 1. The Poggio-Miller-Chang-Harrington-Wu-Tsai (PMCHWT) formulation [20], the combined tangential formulation (CTF) [16], the normal Müller formulation [21, 22], the combined normal formulation (CNF) [16], and the electric and magnetic current combined field integral equation (JMCFIE) [23] are considered.

_{i}## 3. Assessment of numerical results

The accuracy of the aforementioned MoM-based formulations when dealing with plasmonic media in near and far field computations is studied in this section. All the analyses have been focused on spheres in order to properly compare the obtained results with analytical solutions provided by the Mie’s series [24]. Gold, silver and aluminum spheres have been considered and they have been modeled employing the well-known Rao-Wilton-Glisson (RWG) basis functions [25]. The optical properties (relative dielectric permittivity, *ε _{r}*) of these metals at the selected operating wavelength,

*λ*

_{0}= 548.6 nm, have been extracted from [26, 27] (

*ε*= −5.8−

_{r}*j*2.1 for gold,

*ε*= −12.8−

_{r}*j*0.4 for silver, and

*ε*= −35.2−

_{r}*j*9.82 in the case of aluminum). The Galerkin’s testing method has been assumed in this work, meaning that the same RWG basis and testing functions are used. Singular integrals have been accurately evaluated by means of the analytical extraction procedures of [28–31] and Gaussian quadrature rules of 3 points/triangle have been applied for the numerical integration of smooth varying integrands. The MoM solution has been obtained by LU matrix factorization in all cases except for the last subsection concerning the iterative behavior of the considered formulations.

#### 3.1 Near field accuracy validations

The total electric near field has been computed using the five MoM-based integral formulations registered in Table 1 for gold, silver and aluminum spheres with radius *λ*_{0}/2, *λ*_{0}/4, *λ*_{0}/8 and *λ*_{0}/16. The spheres are centered at the origin and illuminated by an $\widehat{x}$polarized plane wave impinging in the $\widehat{z}$ direction (*θ* ^{inc} = 180°). The near fields have been calculated on centered square *XY* planes (z = 0) of side 4*r* with a 100$\times $100 points resolution, *r* being the radius of the spheres. The normalized root mean square (RMS) errors with respect to the Mie’s series results are shown in Fig. 1
and Fig. 2
versus the number of unknowns (mesh size). The error is calculated as follows:

*E**and*

_{Mie}

*E**are the total near fields corresponding to Mie’s series and each formulation, respectively, and*

_{f}*N*is the number of near field samples. Looking at these figures, the first important result is that, in all cases, CTF provides poor results that do not improve with the refinement of the mesh (for small spheres it even get worse). This behavior of CTF in the context of plasmonic materials contrasts with the performance observed in [18, 19] when dealing with LHMs, where CTF showed good accuracy in general lossy cases. It also contrasts with its behavior in conventional dielectrics [16, 17], where it provides greater accuracy than the other tested formulations. Despite combining only tangential equations, just like CTF, PMCHWT works well for these problems. Given that the main difference between PMCHWT and CTF is the set of parameters used to combine the equations, the source of the CTF failure in plasmonics probably lies in the T-EFIE and T-MFIE weighting by coefficients 1/

*η*

_{2}and

*η*

_{2}respectively, which have large imaginary parts in plasmonic media, instead of using 1 and 1 like in PMCHWT. The remaining formulations present in general a good response, being the accuracy of the JMCFIE formulation slightly worse especially for larger spheres with coarser mesh sizes, as can be observed in Fig. 1.

It is worth mentioning that the larger errors observed for low number of unknowns in the smallest spheres are not only due to the formulation, but also to the loss of sphericity given the small number of flat triangular facets used to describe the geometry.

The total near field predicted by the PMCHWT formulation inside and outside the spheres of radius *λ*_{0}/2 and *λ*_{0}/8 is represented in Fig. 3
and Fig. 4
, respectively, to illustrate the great agreement with the Mie’s series solution. It must be pointed out that the agreement is obtained even for a sphere with a strong field inside such as the *λ*_{0}/8 radius sphere of Fig. 4, whose dimensions are close to the resonant size as can be clearly derived from the representation of the absorption, scattering and extinction efficiencies of Fig. 7
(Section 3.2). The results provided by Müller, CNF and JMCFIE formulations (not shown) are very similar, especially those corresponding to Müller and CNF, since results from JMCFIE are somewhat noisy in the case of coarser discretizations. Otherwise, the near field predictions of CTF (not shown) are significantly inaccurate in all cases, with asymmetries and fictitious hot spots near the metallo-dielectric interface.

#### 3.2 Far field validations

The extinction efficiency (normalized cross-section) [24, 32] has been computed using the five integral formulations under consideration for the same spheres used in the previous section. The RMS error of this parameter calculation taking Mie’s series result as the reference has been calculated according (3) and it is represented in Fig. 5
and Fig. 6
. As occurred in the near field computations, it can be seen in both figures that the extinction efficiency calculated by CTF shows high error rates. In contrast, PMCHWT error levels seem to be the most stable values regardless of the material, the size of the sphere and the discretization, demonstrating that it is a very reliable formulation in all cases tested. Regarding CNF, an improvement of its results can be appreciated as the sphere size decreases (Fig. 6). For larger spheres (Fig. 5), however, CNF is less accurate than PMCHWT (except for the aluminum *λ*_{0}/4 sphere). Otherwise, normal Müller and JMCFIE formulations generally have certain lack of regularity, giving small errors in some cases but significantly large in others. Nonetheless, the errors of these formulations are well below CTF errors. As in the near field case, the larger errors for low number of unknowns are also attributable to the loss of sphericity due to the small number of flat triangular facets used to describe the geometry.

The previous results have been obtained for a fixed frequency. Next, the absorption, scattering and extinction efficiencies [24, 32] corresponding to each formulation are shown for different gold spheres of increasing size versus the product of the free-space wavenumber *k*_{0} and the radius *r*. The wavelength (so the permittivity) has been fixed at *λ*_{0} = 548.6 nm and the sphericity of the geometries has been guaranteed by selecting a mesh size of *λ*_{0}/80 for the two smallest spheres, *λ*_{0}/40 for the next four samples, and *λ*_{0}/20 for the last four cases. The efficiencies are displayed together with the Mie’s series results in Fig. 7. CTF shows once again the worst accuracy, while the agreement with the reference of the rest of formulations is quite good, along the lines of the preceding results.

Finally, in order to investigate the ability of the different formulations to accurately handle the plasmon resonance, Fig. 8
depicts the absorption, scattering and extinction efficiencies, as well as the near-field normalized RMS error, as a function of the excitation wavelength around its resonance for a gold sphere with radius 74.712 nm. This sphere corresponds to the case *k*_{0}*r* = 0.8557 in Fig. 7, and it is resonant around 548.6 nm. A very fine wavelength sweep from 500 nm to 600 nm in steps of 0.5 nm has been considered. The dielectric permittivity of gold at these wavelengths has been modeled from [26, 27] using linear interpolation for the intermediate (unavailable) samples. Looking at Fig. 8, it can be observed that all the formulations are free of interior resonance corruption. In accordance with previous results, CTF is in general less accurate, but this cannot be attributed to the internal resonance problem. The observed resonance-free behavior of the formulations could be expected since all of them are combinations of electric and magnetic field integral equations. Otherwise, the analytical extraction procedures of [28–31] used for the accurate evaluation of singular integrals undoubtedly contribute to a large extent to the solid behavior shown by the implementation in all the cases tested.

#### 3.3 Iterative performance

The information obtained from the previous analyses, which have been obtained through matrix factorization and direct solution, is completed next with the observation of the iterative behavior of the five studied formulations. The *λ*_{0}/2 radius sphere was analyzed using MoM solved iteratively by the generalized minimal residual method (GMRES) [33] with a restart parameter of 20 and a maximum number of external iterations of 500, stopped if the normalized residue of 10^{−6} is achieved.

The iterative performance of the formulations is shown in Fig. 9 , where the convergence is illustrated depicting the residue progress versus the number of external iterations required. Solution of CTF shows a bad iterative behavior. Again, this contrasts with its behavior in conventional dielectrics and LHMs [16–19], but it is coherent with the lack of precision shown by this formulation in the previous sections for this kind of problems. Regarding PMCHWT, the known numerical imbalance of the diagonal matrix blocks in this formulation leads to slow convergence and inaccurate results when solving iteratively. This behavior was expected since it was previously reported; see for example [16, 34–36]. A preconditioning technique would be required to rise to the direct solving accuracy. Investigation on proper preconditioning strategies for this and other formulations in the context of plasmonics is left for future work. The rest of formulations show quite good iterative convergence, the fastest one corresponding to the JMCFIE formulation, as happened when dealing with LHMs [18].

## 4. Summary

Five extended MoM-SIE approaches have been applied in the analysis of spheres made of gold, silver and aluminum, which are materials with very different electrical properties that are widely used in diverse plasmonic applications. In view of the results obtained and contrasted with an analytical reference, the first important conclusion is that CTF fails for the analysis of plasmonic problems. This contrasts with the great accuracy provided by this formulation in the analysis of conventional dielectrics and left-handed metamaterials, especially when losses are involved.

The next important observation is that PMCHWT has shown to be a very reliable formulation in plasmonics when using a direct solution. Its error levels seem to be the most stable values both in near and far field calculations, regardless of the material, the size of the sphere, and the geometry mesh refinement. Given that both CTF and PMCHWT formulations combine only tangential integral equations, the accuracy shown by PMCHWT highlights the great influence of the combination parameters on the level of accuracy of the formulations (in addition to its known strong influence on convergence). The main (and well-known) inconvenient of PMCHWT is its lack of convergence in the framework of iterative solutions. Investigation on proper preconditioners may overcome this lack in the future. Regarding the rest of formulations considered, they have shown good accuracy levels in general. The fast iterative convergence of Müller and JMCFIE formulations observed yet in other contexts is also present in the results of this work.

Finally, considering all the above comments we conclude that PMCHWT is preferable for plasmonic problems where factorization and direct solution is feasible. For larger problems, where iterative and/or fast solutions are mandatory, JMCFIE would be a more eligible formulation to use without preconditioning.

## Acknowledgments

This work was supported by the Spanish Government and FEDER (projects TEC2011-28784-C02-01, TEC2011-28784-C02-02, CONSOLIDER-INGENIO2010 CSD2008-00068).

## References and links

**1. **S. A. Maier, *Plasmonics: Fundamentals and Applications* (Springer, 2007).

**2. **J. A. Schuller, E. S. Barnard, W. Cai, Y. C. Jun, J. S. White, and M. L. Brongersma, “Plasmonics for extreme light concentration and manipulation,” Nat. Mater. **9**(3), 193–204 (2010). [CrossRef] [PubMed]

**3. **S. Kawata, ed., *Near-Field Optis and Surface Plasmon Polaritons* (Springer, 2010).

**4. **P. Mühlschlegel, H.-J. Eisler, O. J. F. Martin, B. Hecht, and D. W. Pohl, “Resonant optical antennas,” Science **308**(5728), 1607–1609 (2005). [CrossRef] [PubMed]

**5. **L. Novotny and N. F. van Hulst, “Antennas for Light,” Nat. Photonics **5**(2), 83–90 (2011). [CrossRef]

**6. **P. Biagioni, J. S. Huang, and B. Hecht, “Nanoantennas for visible and infrared radiation,” Physics Optics arXiv:1103.1568v1, (2011).

**7. **J. Zhou, Y. Zhou, S. L. Ng, H. X. Zhang, W. X. Que, Y. L. Lam, Y. C. Chan, and C. H. Kam, “Three-dimensional photonic band gap structure of a polymer-metal composite,” Appl. Phys. Lett. **76**, 3337–3339 (2000).

**8. **B. T. Schwartz and R. Piestun, “Total external reflection from metamaterials with ultralow refractive index,” J. Opt. Soc. Am. B **20**(12), 2448–2453 (2003). [CrossRef]

**9. **M. Salaün, B. Corbett, S. B. Newcomb, and M. E. Pemble, “Fabrication and characterization of three-dimensional silver/air inverted opal photonic crystals,” J. Mater. Chem. **20**(36), 7870–7874 (2010). [CrossRef]

**10. **A. Taflove and M. E. Brodwin, “Numerical solution of steadystate electromagnetic scattering problems using the timedependent Maxwell’s equations,” IEEE Trans. Microw. Theory Tech. **23**(8), 623–630 (1975). [CrossRef]

**11. **T. Weiland, “A discretization method for the solution of Maxwell’s equations for six-component fields,” Arch. Elektron. Übertragungstech. **31**, 116–120 (1977).

**12. **R. F. Harrington, *Field Computation by Moment Method* (IEEE Press, 1993).

**13. **A. M. Kern and O. J. F. Martin, “Surface integral formulation for 3D simulations of plasmonic and high permittivity nanostructures,” J. Opt. Soc. Am. A **26**(4), 732–740 (2009). [CrossRef] [PubMed]

**14. **B. Gallinet, A. M. Kern, and O. J. F. Martin, “Accurate and versatile modeling of electromagnetic scattering on periodic nanostructures with a surface integral approach,” J. Opt. Soc. Am. A **27**(10), 2261–2271 (2010). [CrossRef] [PubMed]

**15. **J. M. Taboada, J. Rivero, F. Obelleiro, M. G. Araújo, and L. Landesa, “Method-of-moments formulation for the analysis of plasmonic nano-optical antennas,” J. Opt. Soc. Am. A **28**(7), 1341–1348 (2011). [CrossRef] [PubMed]

**16. **P. Ylä-Oijala, M. Taskinen, and S. Järvenpää, “Surface integral equation formulations for solving electromagnetic scattering problems with iterative methods,” Radio Sci. **40**(6), RS6002 (2005). [CrossRef]

**17. **Ö. Ergül and L. Gürel, “Comparison of integral-equation formulations for the fast and accurate solution of scattering problems involving dielectric objects with the multilevel fast multipole algorithm,” IEEE Trans. Antenn. Propag. **57**(1), 176–187 (2009). [CrossRef]

**18. **M. G. Araújo, J. M. Taboada, J. Rivero, and F. Obelleiro, “Comparison of Surface Integral Equations for Left-Handed Materials,” Prog. Electromagn. Res. **118**, 425–440 (2011). [CrossRef]

**19. **Ö. Ergül, “Fast and Accurate Analysis of Homogenized Metamaterials With the Surface Integral Equations and the Multilevel Fast Multipole Algorithm,” IEEE Antennas Wirel. Propag. Lett. **10**, 1286–1289 (2011). [CrossRef]

**20. **J. R. Mautz and R. F. Harrington, “Electromagnetic scattering from a homogeneous material body of revolution,” Arch. Elektr. Uebertrag. **33**, 71–80 (1979).

**21. **C. Müller, *Foundations of the Mathematical Theory of Electromagnetic Waves* (Springer, 1969).

**22. **P. Ylä-Oijala and M. Taskinen, “Well-conditioned Müller formulation for electromagnetic scattering by dielectric objects,” IEEE Trans. Antenn. Propag. **53**(10), 3316–3323 (2005). [CrossRef]

**23. **P. Ylä-Oijala and M. Taskinen, “Application of combined field integral equation for electromagnetic scattering by dielectric and composite objects,” IEEE Trans. Antenn. Propag. **53**(3), 1168–1173 (2005). [CrossRef]

**24. **C. F. Bohren and D. R. Huffman, *Absorption and Scattering of Light by Small Particles* (John Wiley, 1983).

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

**26. **P. B. Johnson and R. W. Christy, “Optical constants of the noble metals,” Phys. Rev. B **6**(12), 4370–4379 (1972). [CrossRef]

**27. **G. Hass and J. E. Waylonis, “Optical Constants and Reflectance and Transmittance of Evaporated Aluminum in the Visible and Ultraviolet,” J. Opt. Soc. Am. **51**(7), 719–722 (1961). [CrossRef]

**28. **D. R. Wilton, S. M. Rao, A. W. Glisson, D. H. Schaubert, O. M. Al-Bundak, and C. M. Butler, “Potential integrals for uniform and linear source distributions on polygonal and polyhedral domains,” IEEE Trans. Antenn. Propag. **32**(3), 276–281 (1984). [CrossRef]

**29. **R. E. Hodges and Y. Rahmat-Samii, “The evaluation of MFIE integrals with the use of vector triangle basis functions,” Microw. Opt. Technol. Lett. **14**(1), 9–14 (1997). [CrossRef]

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

**31. **P. Ylä-Oijala and M. Taskinen, “Calculation of CFIE impedance matrix elements with RWG and$\widehat{n}\times $ RWG functions,” IEEE Trans. Antenn. Propag. **51**(8), 1837–1846 (2003). [CrossRef]

**32. **A. Ishimaru, *Electromagnetic Wave Propagation, Radiation and Scattering* (Prentice-Hall, 1991).

**33. **Y. Saad and M. Schultz, “Gmres: A generalized minimal residual algorithm for solving nonsymmetric linear systems,” SIAMJ. Sci. Statist. Comput. **7**(3), 856–869 (1986). [CrossRef]

**34. **T. W. Lloyd, J. M. Song, and M. Yang, “Numerical study of surface integral formulations for low-contrast objects,” IEEE Antennas Wirel. Propag. Lett. **4**(1), 482–485 (2005). [CrossRef]

**35. **P. Ylä-Oijala and M. Taskinen, “Improving conditioning of electromagnetic surface integral equations using normalized field quantities,” IEEE Trans. Antenn. Propag. **55**(1), 178–185 (2007). [CrossRef]

**36. **Y. A. Liu and W. C. Chew, “Stability of surface integral equation for left-handed materials,” IEEE Trans. Microw. Antennas Propag. **1**(1), 84–89 (2007). [CrossRef]