Integral-Equation Formulations of Plasmonic Problems in the Visible Spectrum and Beyond

Abstract

Computational modeling of nano-plasmonic structures is essential to understand their electrodynamic responses before experimental efforts in measurement setups. Similar to the other ranges of the electromagnetic spectrum, there are alternative methods for the numerical analysis of nano-plasmonic problems, while the optics literature is dominated by differential equations that require discretizations of the host media with artificial truncations. These approaches often need serious assumptions, such as periodicity, infinity, or self-similarity, in order to reduce the computational load. On the other hand, surface integral equations based on integro-differential operators can bring important advantages for accurate and efficient modeling of nano-plasmonic problems with arbitrary geometries. Electrical properties of materials, which may be obtained either experimentally or via physical modeling, can easily be inserted into integral-equation formulations, leading to accurate predictions of electromagnetic responses of complex structures. This chapter presents the implementation of such accurate, efficient, and reliable solvers based on appropriate combinations of surface integral equations, discretizations, numerical integrations, fast algorithms, and iterative techniques. As a case study, nanowire transmission lines are investigated in wide-frequency ranges, demonstrating the capabilities of the developed implementations.

1. Introduction

As in all areas of electrodynamics, numerical study of plasmonic problems is essential to understand interactions between electromagnetic waves and matter at the higher range of the spectrum. Applications include nanowires for negative refraction, imaging, and super-resolution [1, 2], and nanoantennas for energy harvesting, single-molecule sensing, and optical links [3–9], to name a few. At optical frequencies, some metals are known to possess strong plasmonic properties [10] that are crucial for a majority of such applications, while their accurate analysis requires more than perfectly conducting models that are common in radio and microwave regimes. In the infrared region, it may not be obvious when perfect conductivity or impedance approximation methods can safely be used. Hence, it is desirable to extend the plasmonic-modeling capabilities across wide ranges of frequencies until they converge to the other forms. While, in the literature, experimental studies are often supported by differential solvers, their applicability to complex problems is usually limited to small-scale and/or simplified models due to well-known drawbacks, such as need for space (host-medium) discretizations that are accompanied with artificial truncations. Major tools of computational electromagnetics, that is, surface integral equations [11, 12] employing integro-differential operators, are recently applied to plasmonic problems with promising results for realistic simulations of complex structures [13–23]. In fact, surface integral equations need only the discretization of boundaries between different media, which usually correspond to the surface of the plasmonic object. In addition to homogeneous bodies, they are also applicable to piecewise homogeneous cases, making it possible to analyze structures with coexisting multiple materials [24].

Using surface integral equations, it is possible to solve plasmonic problems involving finite models with arbitrary geometries, without periodicity, self-similarity, and infinity assumptions. When the object is large in terms of wavelength, fast and efficient methods, such as the multilevel fast multipole algorithm (MLFMA) [25], are available to accelerate solutions [26–28]. For plasmonic modeling, effective permittivity values with negative real parts are required, while they are already available via theoretical and experimental studies [10]. In the phasor domain with time-harmonic sources, which is considered in this chapter, permittivity is a simulation parameter with a fixed value at a given frequency. Then, frequency sweeps can be performed by using the discrete values of the permittivity with respect to frequency. As theoretical models, Drude (D) or Lorentz-Drude (LD) models are commonly used. While these models (especially the Lorentz-Drude model) provide reliable permittivity values in wide-frequency ranges, they deviate from experimental data at higher frequencies of the optical spectrum. From the perspective of surface integral equations, it does not matter where the permittivity values are obtained from. Besides, there is a great flexibility in geometric modeling, allowing sharp edges and corners, tips, and subwavelength details [29]. On top of these, the background of surface integral equations provides self-consistency and accuracy-check mechanisms, such as based on the equivalence theorem, enabling accuracy analysis without resorting to alternative solvers [30].

From numerical point of view, surface integral equations bring their own challenges when they are applied to plasmonic problems. In free space, plasmonic objects are naturally high-contrast problems [15], leading to difficulties in maintaining the accuracy and/or efficiency. Considering the equivalence theorem, ideal mesh size for surface formulations can be selected based on wavenumber of the host medium, where the impressed sources are located [26]. Therefore, the source of the inaccuracy is not directly the discretization size, but a combination of geometric deviation (for smooth objects), numerical integration, and imbalanced contributions from inner/outer media. Efficiency of iterative solutions may also deteriorate due to imbalanced matrix blocks that lead to ill-conditioned matrix equations [31]. On the other hand, numerical challenges are not only due to the high contrasts of plasmonic objects. The effective permittivity of a plasmonic medium is typically negative, which becomes increasingly large at lower frequencies. In numerical solutions, integro-differential operators become localized with exponentially decaying Green’s function. This localization is responsible for the evolution of plasmonic formulations into perfectly conducting types, while this process may not be achieved smoothly in discrete forms. Some traditional formulations break down due to dominant inner contributions, which are difficult to compute accurately [32], if not impossible. Classical singularity extractions may fail to provide smooth integrands, leading to increasingly inaccurate near-zone interactions. While all formulations may be improved by manipulating integrations into more suitable forms, our focus is to develop new formulations that reduce into perfectly conducting formulations in the limit. All results presented in this chapter are obtained by such a stabilized integral-equation formulation, namely a modified combined tangential formulation (MCTF), which provides accurate results using the conventional Rao-Wilton-Glisson (RWG) discretizations [33].

The chapter is organized as follows. In Section 2, we present surface integral equations, with the emphasis on MCTF. Discretization is presented in Section 3, including implementation details that may be followed by the readers to develop their own solvers. MLFMA is further discussed in Section 4, demonstrating how to accelerate numerical solutions. Finally, we present an extensive case study, involving nanowire transmission lines in a wide range of frequency to illustrate the significant differences between the analytical models and measurement data for the permittivity values. In the following, time-harmonic electrodynamic problems are considered with exp(− iωt) time dependency, where i² = −1 and ω = 2πf is the angular frequency.

2. Surface integral equations

For deriving surface formulations, we consider a plasmonic object with permittivity/permeability (${\epsilon }_p/{\mu }_p$) located in unbounded free space with permittivity/permeability (${\epsilon }_o/{\mu }_o$). Alternative surface integral equations can be obtained by considering the boundary conditions on the surface of the object. In a general form, we have

$\left[\begin{array}{cc}{\mathcal{Z}}_{11}& {\mathcal{Z}}_{12}\\ {\mathcal{Z}}_{21}& {\mathcal{Z}}_{22}\end{array}\right]\cdot \left[\begin{array}{c}J\\ M\end{array}\right](r)=\left[\begin{array}{c}a{\displaystyle \widehat{n}}\times {\displaystyle \widehat{n}}\times E^{\text{inc}}-e{\displaystyle \widehat{n}}\times H^{\text{inc}}\\ c{\displaystyle \widehat{n}}\times {\displaystyle \widehat{n}}\times H^{\text{inc}}+g{\displaystyle \widehat{n}}\times E^{\text{inc}}\end{array}\right](r),$E1

where $J={\displaystyle \widehat{n}}\times H$ and $M=-{\displaystyle \widehat{n}}\times E$ are the equivalent currents written in terms of the tangential electric field intensity E and the magnetic field intensity H on the closed surface ($r\in S$). In the above, ${\displaystyle \widehat{n}}$ is the unit vector outward the object, and $E^{\text{inc}}$ and $H^{\text{inc}}$ are the incident electric and magnetic fields, respectively, created by impressed sources located in the host medium. At an observation point on a locally planar surface (solid angle = 2π), the combined operators can be written as

${\mathcal{Z}}_{11}=-{\displaystyle \widehat{n}}\times {\displaystyle \widehat{n}}\times (a{\eta }_o{\mathcal{T}}_o+b{\eta }_p{\mathcal{T}}_p)+{\displaystyle \widehat{n}}\times (e{\mathcal{K}}_{\text{PV}\mathrm{,}o}-f{\mathcal{K}}_{\text{PV}\mathrm{,}p})-(e+f)\mathcal{I}/2$E2

${\mathcal{Z}}_{12}={\displaystyle \widehat{n}}\times {\displaystyle \widehat{n}}\times (a{\mathcal{K}}_{\text{PV}\mathrm{,}o}+b{\mathcal{K}}_{\text{PV}\mathrm{,}p})-(a-b){\displaystyle \widehat{n}}\times \mathcal{I}/2+{\displaystyle \widehat{n}}\times (e{\eta }_o^{-1}{\mathcal{T}}_o-f{\eta }_p^{-1}{\mathcal{T}}_p)$E3

${\mathcal{Z}}_{21}=-{\displaystyle \widehat{n}}\times {\displaystyle \widehat{n}}\times (c{\mathcal{K}}_{\text{PV}\mathrm{,}o}+d{\mathcal{K}}_{\text{PV}\mathrm{,}p})+(c-d){\displaystyle \widehat{n}}\times \mathcal{I}/2-{\displaystyle \widehat{n}}\times (g{\eta }_o{\mathcal{T}}_o-h{\eta }_p{\mathcal{T}}_p)$E4

${\mathcal{Z}}_{22}=-{\displaystyle \widehat{n}}\times {\displaystyle \widehat{n}}\times (c{\eta }_o^{-1}{\mathcal{T}}_o+d{\eta }_p^{-1}{\mathcal{T}}_p)+{\displaystyle \widehat{n}}\times (g{\mathcal{K}}_{\text{PV}\mathrm{,}o}-h{\mathcal{K}}_{\text{PV}\mathrm{,}p})-(g+h)\mathcal{I}/\mathrm{2,}$E5

where $\left\{a\mathrm{,}b\mathrm{,}c\mathrm{,}d\mathrm{,}e\mathrm{,}f\mathrm{,}g\mathrm{,}h\right\}$ are generalized coefficients. In the above, ${\eta }_o=\sqrt{{\mu }_o}/\sqrt{{\epsilon }_o}$ is the intrinsic impedance of the host medium, whereas ${\eta }_p=\sqrt{{\mu }_p}/\sqrt{{\epsilon }_p}$ is the complex intrinsic impedance of the plasmonic object. The integro-differential and identity operators are derived as

${\mathcal{T}}_u\left\{X\right\}(r)=i{k}_u{\displaystyle {\int }_S}d{r}^{\prime }[X(r^{\prime })+\frac{1}{k_u^2}{\nabla }^{\prime }\cdot X(r^{\prime })\nabla ]g_u(r\mathrm{,}{r}^{\prime })$E6

${\mathcal{K}}_{PV\mathrm{,}u}\left\{X\right\}(r)={\displaystyle {\int }_{\text{PV}\mathrm{,}S}}d{r}^{\prime }X(r^{\prime })\times {\nabla }^{\prime }{g}_u(r\mathrm{,}{r}^{\prime })$E7

$\mathcal{I}\left\{X\right\}(r)=X(r)$E8

for $r\in S$, where $\text{PV}$ indicates the principal value of the integral, $\nabla ={\displaystyle \widehat{x}}\partial /\partial x+{\displaystyle \widehat{y}}\partial /\partial y+{\displaystyle \widehat{z}}\partial /\partial z$ is the differential operator, $g_u(r\mathrm{,}{r}^{\prime })=\exp (i{k}_u|r-r^{\prime }|)/(4\pi |r-r^{\prime }|)$ is the homogeneous-space Green’s function, and $k_u=2\pi /{\lambda }_u=\omega \sqrt{{\mu }_u{\epsilon }_u}$ is the wavenumber for $u=\left\{o\mathrm{,}p\right\}$.

The conventional formulations can be obtained by setting the generalized coefficients to suitable values such that the outer and inner problems are coupled while the internal resonances are removed. By using nonzero values for {$e\mathrm{,}f\mathrm{,}g\mathrm{,}h$} while setting $\left\{a\mathrm{,}b\mathrm{,}c\mathrm{,}d\right\}$ to zero leads to N-formulations, such as the Müller formulation and the combined normal formulation [12]. These formulations contain the identity operator $\mathcal{I}$, which usually dominates the matrix equations when Galerkin discretization is used. Therefore, matrix equations derived from N-formulations are generally easier to solve iteratively. On the other hand, T-formulations are obtained by selecting $\left\{a\mathrm{,}b\mathrm{,}c\mathrm{,}d\right\}$ nonzero, while inserting zero values for $\left\{e\mathrm{,}f\mathrm{,}g\mathrm{,}h\right\}$. The Poggio-Miller-Chang-Harrington-Wu-Tsai formulation [34] and the combined tangential formulation [12] are among the well-known T-formulations. As opposed to N-formulations, T-formulations contain either the rotational identity operator ${\displaystyle \widehat{n}}\times \mathcal{I}$ or no identity operator at all (when $a=b$ and $c=d$). Hence, using a Galerkin discretization, T-formulations do not contain a dominant identity operator and they produce matrix equations that are potentially ill-conditioned. Finally, when a mixture of coefficients are used from the sets $\left\{a\mathrm{,}b\mathrm{,}c\mathrm{,}d\right\}$ and {e, f, g, h}, mixed formulations are obtained. For example, the JM combined-field integral equation [35] is a mixed formulation when all coefficients are nonzero. Obviously, mixed formulations always contain a dominant identity operator (due to either $\mathcal{I}$ or ${\displaystyle \widehat{n}}\times \mathcal{I}$).

Discretization is an important stage of numerical solutions. All formulations described above can be discretized in different ways such that the derived matrix equations can be well conditioned, and, at the same time, they may produce accurate results. On the other hand, using a Galerkin scheme employing the same set of basis and testing functions, N-formulations and mixed formulations usually produce better-conditioned matrix equations than T-formulations, as mentioned above. In addition, when low-order discretizations are used, the existence of a dominant identity operator is critical in terms of accuracy. It is well known that a discretized identity operator acts like a discretized integro-differential operator with a Dirac-delta kernel [36]. Therefore, a low-order discretization of the identity operator may produce large errors, leading to inaccurate results if the operator is directly tested such that it dominates the matrix equation. RWG discretizations of N-formulations and mixed formulations have this serious drawback, making them less preferred (despite their faster iterative solutions) in comparison to T-formulations in many applications. The tradeoff between the efficiency and accuracy has been resolved in many studies [37] by improving the accuracy of N-formulations and mixed formulations via alternative discretizations and/or by improving the efficiency of T-formulations via preconditioning.

In the context of plasmonic problems, further challenges appear in surface formulations. First, considering that their permittivity values can be written as ${\epsilon }_p={\epsilon }_o(-{\epsilon }_R+i{\epsilon }_I)$, where both ${\epsilon }_R$ and ${\epsilon }_I$ are positive, plasmonic objects are naturally high-contrast structures in free space (except for very high frequencies for which $-{\epsilon }_R\to 1$). Then, the matrix equations derived from surface formulations can be unbalanced, leading to efficiency and/or accuracy problems. For planar discretizations of curved surfaces, fine discretizations are needed to capture the geometry of the object. At lower frequencies of the optical range, ${\epsilon }_R$ can be very large (as large as 1000 and beyond) such that the localization of the operators as ${\mathcal{T}}_p\to -\mathcal{I}/2$ and ${\mathcal{K}}_{\text{PV}\mathrm{,}p}-\mathcal{I}/2\to -\mathcal{I}/2$ when ${\epsilon }_R\to \infty$ leads to numerical problems if the blocks are not weighted properly (that occurs in many conventional formulations). While the well-known perfectly conducting models may be used at lower frequencies, it may not be obvious where the plasmonic model can be omitted for a given structure. Hence, it is desirable to extend the applicability of the surface integral equations in wide-frequency ranges until other kinds of approaches can safely be used. In a recent study, we show that a new tangential formulation, namely MCTF, provides reliable and convergent solutions in wide ranges of frequencies of the optical spectrum [32]. Considering the general form, MCTF is obtained by using $a=b=1$ and $c=d={\eta }_o{\eta }_p$, while setting $e=f=g=h=0$. Therefore, we obtain

${\mathcal{Z}}_{11}^{\text{MCTF}}=-{\displaystyle \widehat{n}}\times {\displaystyle \widehat{n}}\times ({\eta }_o{\mathcal{T}}_o+{\eta }_p{\mathcal{T}}_p)$E9

${\mathcal{Z}}_{12}^{\text{MCTF}}={\displaystyle \widehat{n}}\times {\displaystyle \widehat{n}}\times ({\mathcal{K}}_{\text{PV}\mathrm{,}o}+{\mathcal{K}}_{\text{PV}\mathrm{,}p})$E10

${\mathcal{Z}}_{21}^{\text{MCTF}}=-{\displaystyle \widehat{n}}\times {\displaystyle \widehat{n}}\times {\eta }_o{\eta }_p({\mathcal{K}}_{\text{PV}\mathrm{,}o}+{\mathcal{K}}_{\text{PV}\mathrm{,}p})$E11

${\mathcal{Z}}_{22}^{\text{MCTF}}=-{\displaystyle \widehat{n}}\times {\displaystyle \widehat{n}}\times ({\eta }_p{\mathcal{T}}_o+{\eta }_o{\mathcal{T}}_p).$E12

It can be observed that MCTF is completely free of the identity operator, and it can be shown that it smoothly turns into the electric-field integral equation for perfectly conducting objects as the frequency drops and ${\epsilon }_R$ goes to infinity. In the following, we consider numerical solutions of plasmonic problems formulated with MCTF.

3. Discretization

Similar to the diversity of surface integral equations, discretization can be performed in alternative ways. Using a Galerkin scheme, the basis and testing functions are selected as the same set of $N$ functions locally defined on the surface. As a popular choice for triangular discretizations, which is also considered in this chapter, the RWG functions are defined as [33]

${\mathit{f}}_n(r)=\{\begin{array}{c}\begin{array}{rr}\hfill \frac{l_n}{2A_{n1}}(r-r_{n1}),& \hfill r\in S_{n1}\\ \hfill \frac{l_n}{2A_{n2}}(r_{n2}-r),& \hfill r\in S_{n2}\\ \hfill \mathrm{0,}& \hfill r\notin S_n\mathrm{.}\end{array}\end{array}$E13

Each RWG function is located on a pair of triangles sharing an edge. In the above, $l_n$ represents the length of the main edge, $A_{n1}$ and $A_{n2}$ are, respectively, the areas of the first ($S_{n1}$) and the second ($S_{n2}$) triangles, and $r_{n1}$ and $r_{n2}$ represent the coordinates of the nodes opposite of the edge. The RWG functions are divergence conforming and their divergence is finite everywhere, that is,

$\nabla \cdot f_n(r)=\{\begin{array}{c}\begin{array}{rr}\hfill \frac{l_n}{A_{n1}}\mathrm{,}& \hfill r\in S_{n1}\\ \hfill -\frac{l_n}{A_{n2}}\mathrm{,}& \hfill r\in S_{n2}\\ \hfill \mathrm{0,}& \hfill r\notin S_n\mathrm{,}\end{array}\end{array}$E14

while the charge neutrality is satisfied locally as $A_{n1}{l}_n/A_{n1}-A_{n2}{l}_n/A_{n2}=0$.

By selecting the basis and testing functions ($b_n$ and $t_m$ for $\left\{n\mathrm{,}m\right\}=\{1,2,\dots \mathrm{,}N\}$) as the same set of the RWG functions, MCTF can be discretized as

$\left[\begin{array}{cc}{{\displaystyle \overline{Z}}}_{11}^{\text{MCTF}}& {{\displaystyle \overline{Z}}}_{12}^{\text{MCTF}}\\ {{\displaystyle \overline{Z}}}_{21}^{\text{MCTF}}& {{\displaystyle \overline{Z}}}_{22}^{\text{MCTF}}\end{array}\right]\cdot \left[\begin{array}{c}{a}_J\\ a_M\end{array}\right]=\left[\begin{array}{c}{w}_1^{\text{MCTF}}\\ w_2^{\text{MCTF}}\end{array}\right]\mathrm{,}$E15

where $a_J$ and $a_M$ are vectors containing complex coefficients to expand the current densities. The matrix elements and the elements of the right-hand-side vector are derived as

${{\displaystyle \overline{Z}}}_{11}^{\text{MCTF}}={\eta }_o{{\displaystyle \overline{T}}}_o^T+{\eta }_p{{\displaystyle \overline{T}}}_p^T$E16

${{\displaystyle \overline{Z}}}_{12}^{\text{MCTF}}=-{{\displaystyle \overline{K}}}_{\text{PV}\mathrm{,}o}^T-{{\displaystyle \overline{K}}}_{\text{PV}\mathrm{,}p}^T$E17

${{\displaystyle \overline{Z}}}_{21}^{\text{MCTF}}={\eta }_o{\eta }_p\left({{\displaystyle \overline{K}}}_{\text{PV}\mathrm{,}o}^T+{{\displaystyle \overline{K}}}_{\text{PV}\mathrm{,}p}^T\right)$E18

${{\displaystyle \overline{Z}}}_{22}^{\text{MCTF}}={\eta }_p{{\displaystyle \overline{T}}}_o^T+{\eta }_o{{\displaystyle \overline{T}}}_p^T$E19

and

$w_1^{\text{MCTF}}=-{\displaystyle {\int }_{S_m}}dr{t}_m(r)\cdot E^{\text{inc}}(r)$E20

$w_2^{\text{MCTF}}=-{\eta }_o{\eta }_p{\displaystyle {\int }_{S_m}}dr{t}_m(r)\cdot H^{\text{inc}}(r),$E21

respectively. Furthermore, the discretized operators can be written as

${{\displaystyle \overline{T}}}_u^T[m\mathrm{,}n]=i{k}_u{\displaystyle {\int }_{S_m}}dr{t}_m(r)\cdot {\displaystyle {\int }_{S_n}}dr\prime g_u(r\mathrm{,}r\prime )b_n(r\prime )+\frac{i}{k_u}{\displaystyle {\int }_{S_m}}dr{t}_m(r)\cdot {\displaystyle {\int }_{S_n}}dr\prime \nabla g_u(r\mathrm{,}r\prime ){\nabla }^{\prime }\cdot b_n(r\prime )$E22

${{\displaystyle \overline{K}}}_{\text{PV}\mathrm{,}u}^T[m\mathrm{,}n]={\displaystyle {\int }_{S_m}}dr{\mathit{t}}_m(r)\cdot {\displaystyle {\int }_{\text{PV}\mathrm{,}{S}_n}}dr\prime {\mathit{b}}_n(r\prime )\times {\nabla }^{\prime }{g}_u(r\mathrm{,}r\prime ),$E23

where the integrals are evaluated on the supports of the testing and basis functions ($S_m$ and $S_n$).

At this stage, we can consider the interaction of two half RWG functions associated with the ath triangle of the mth edge and bth triangle of the nth edge, respectively ($\left\{a\mathrm{,}b\right\}=\{1,2\}$). One can obtain

$\begin{array}{c}{{\displaystyle \overline{T}}}_u^T[m\mathrm{,}n\mathrm{,}a\mathrm{,}b]=\frac{{\gamma }_{ma}{\gamma }_{nb}{l}_m{l}_n}{4}i{k}_u\frac{1}{A_{ma}}{\displaystyle {\int }_{S_{ma}}}dr(r-r_{ma})\cdot \frac{1}{A_{nb}}{\displaystyle {\int }_{S_{nb}}}dr\prime (r\prime -r_{nb})g_u(r\mathrm{,}r\prime )\\ -{\gamma }_{ma}{\gamma }_{nb}{l}_m{l}_n\frac{i}{k_u}\frac{1}{A_{ma}}{\displaystyle {\int }_{S_{ma}}}dr\frac{1}{A_{nb}}{\displaystyle {\int }_{S_{nb}}}dr\prime g_u(r\mathrm{,}r\prime )\end{array}$E24

${{\displaystyle \overline{\mathit{K}}}}_{\text{PV}\mathrm{,}u}^T[m\mathrm{,}n\mathrm{,}a,b]=\frac{{\gamma }_{ma}{\gamma }_{nb}{l}_m{l}_n}{4}\frac{1}{A_{ma}}{\displaystyle {\int }_{S_{ma}}}dr(r-r_{ma})\cdot (r-r_{nb})\times \frac{1}{A_{nb}}{\displaystyle {\int }_{\text{PV}\mathrm{,}{S}_{nb}}}dr\prime {\nabla }^{\prime }{g}_u(r\mathrm{,}r\prime ),$E25

where ${\gamma }_{nb}\mathrm{,}{\gamma }_{ma}=\pm 1$, depending on the direction of the basis and testing functions on triangles. For the integrations on the testing and basis triangles, alternative methods can be used. Applying Gaussian quadrature is common in the literature, if the singularity of Green’s function is extracted from the inner integrals. In any case, the integration methods used on the testing and basis triangles do not have to be the same, that is, different sampling schemes can be used. For the sake of brevity, we consider a single-point testing scheme by using the center point of each triangle $r_{ma}^{\text{cr}}$, leading to

$\begin{array}{ll}{{\displaystyle \overline{\mathit{T}}}}_u^T[m\mathrm{,}n\mathrm{,}a\mathrm{,}b]\hfill & =\frac{{\gamma }_{ma}{\gamma }_{nb}{l}_m{l}_n}{4}i{k}_u(r_{ma}^{\text{cr}}-r_{ma})\cdot \frac{1}{A_{nb}}{\displaystyle {\int }_{S_{nb}}}dr\prime (r\prime -r_{nb})g_u(r_{ma}^{\text{cr}}\mathrm{,}r\prime )\hfill \\ \hfill & -{\gamma }_{ma}{\gamma }_{nb}{l}_m{l}_n\frac{i}{k_u}\frac{1}{A_{nb}}{\displaystyle {\int }_{S_{nb}}}dr\prime g_u(r_{ma}^{\text{cr}}\mathrm{,}r\prime )\hfill \end{array}$E26

$\begin{array}{ll}{{\displaystyle \overline{T}}}_u^T[m\mathrm{,}n\mathrm{,}a\mathrm{,}b]\hfill & =\frac{{\gamma }_{ma}{\gamma }_{nb}{l}_m{l}_n}{4}i{k}_u({\mathit{\rho }}_{ma}^{\text{cr}}-{\mathit{\rho }}_{ma})\cdot \frac{1}{A_{nb}}{\displaystyle {\int }_{S_{nb}}}dr\prime (\mathit{\rho }\prime -{\mathit{\rho }}_{ma}^{\text{cr}})g_u(r_{ma}^{\text{cr}}\mathrm{,}r\prime )\hfill \\ \hfill & +\frac{{\gamma }_{ma}{\gamma }_{nb}{l}_m{l}_n}{4}i{k}_u({\mathit{\rho }}_{ma}^{\text{cr}}-{\mathit{\rho }}_{ma})\cdot ({\mathit{\rho }}_{ma}^{\text{cr}}-{\mathit{\rho }}_{nb})\frac{1}{A_{nb}}{\displaystyle {\int }_{S_{nb}}}dr\prime g_u(r_{ma}^{\text{cr}}\mathrm{,}r\prime )\hfill \\ \hfill & -{\gamma }_{ma}{\gamma }_{nb}{l}_m{l}_n\frac{i}{k_u}\frac{1}{A_{nb}}{\displaystyle {\int }_{S_{nb}}}dr\prime g_u(r_{ma}^{\text{cr}}\mathrm{,}r\prime )\hfill \end{array}$E27

${{\displaystyle \overline{K}}}_{\text{PV}\mathrm{,}u}^T[m\mathrm{,}n\mathrm{,}a,b]=\frac{{\gamma }_{ma}{\gamma }_{nb}{l}_m{l}_n}{4}(r_{ma}^{\text{cr}}-r_{ma})\cdot (r_{ma}^{\text{cr}}-r_{nb})\times \frac{1}{A_{nb}}{\displaystyle {\int }_{\text{PV}\mathrm{,}{S}_{nb}}}dr\prime {\nabla }^{\prime }{g}_u(r_{ma}^{\text{cr}}\mathrm{,}r\prime ),$E28

where $\{\mathit{\rho }\prime \mathrm{,}{\mathit{\rho }}_{ma}\mathrm{,}{\mathit{\rho }}_{nb}\mathrm{,}{\mathit{\rho }}_{ma}^{\text{cr}}\}$ represent the projections of $\left\{r^{\prime }\mathrm{,}{r}_{ma}\mathrm{,}{r}_{nb}\mathrm{,}{r}_{ma}^{\text{cr}}\right\}$ onto the basis plane.

It is generally more efficient to compute the interactions via triangle by triangle (rather than RWG by RWG) since common integrals related to a basis triangle can be evaluated once and used in multiple interactions related to the triangle. For MCTF, interactions are calculated (for $a=\mathrm{1,2}$ and $b=\mathrm{1,2}$, and $u=o\mathrm{,}p$) as

${{\displaystyle \overline{Z}}}_{11}^{\text{MCTF}}[m\mathrm{,}n]\leftarrow \frac{{\gamma }_{ma}{\gamma }_{nb}{l}_m{l}_n}{4}i{k}_u{\eta }_u\left\{({\mathit{\rho }}_{ma}^{\text{cr}}-{\mathit{\rho }}_{ma})\cdot \left[I_{ma\mathrm{,}nb\mathrm{,}u}^B+({\mathit{\rho }}_{ma}^{\text{cr}}-{\mathit{\rho }}_{nb})I_{ma\mathrm{,}nb\mathrm{,}u}^A\right]-\frac{4}{k_u^2}{I}_{ma\mathrm{,}nb\mathrm{,}u}^A\right\}$E29

${{\displaystyle \overline{Z}}}_{22}^{\text{MCTF}}[m\mathrm{,}n]\leftarrow \frac{{\gamma }_{ma}{\gamma }_{nb}{l}_m{l}_n}{4}i{k}_u\frac{{\eta }_o{\eta }_p}{{\eta }_u}\left\{({\mathit{\rho }}_{ma}^{\text{cr}}-{\mathit{\rho }}_{ma})\cdot \left[I_{ma\mathrm{,}nb\mathrm{,}u}^B+({\mathit{\rho }}_{ma}^{\text{cr}}-{\mathit{\rho }}_{nb})I_{ma\mathrm{,}nb\mathrm{,}u}^A\right]-\frac{4}{k_u^2}{I}_{ma\mathrm{,}nb\mathrm{,}u}^A\right\}$E30

${{\displaystyle \overline{Z}}}_{12}^{\text{MCTF}}[m\mathrm{,}n]\leftarrow -\frac{{\gamma }_{ma}{\gamma }_{nb}{l}_m{l}_n}{4}(r_{ma}^{\text{cr}}-r_{ma})\cdot \left[(r_{ma}^{\text{cr}}-r_{nb})\times I_{ma\mathrm{,}nb\mathrm{,}u}^{C\mathrm{,}\text{PV}}\right]$E31

${{\displaystyle \overline{Z}}}_{21}^{\text{MCTF}}[m\mathrm{,}n]\leftarrow \frac{{\gamma }_{ma}{\gamma }_{nb}{l}_m{l}_n}{4}{\eta }_o{\eta }_p(r_{ma}^{\text{cr}}-r_{ma})\cdot \left[(r_{ma}^{\text{cr}}-r_{nb})\times I_{ma\mathrm{,}nb\mathrm{,}u}^{C\mathrm{,}\text{PV}}\right]\mathrm{,}$E32

where $\leftarrow$ indicates the update operation. Each matrix element is obtained by combining the contributions of four triangle-triangle interactions. By using triangle-triangle interactions, a basis integral ($I_{ma\mathrm{,}nb\mathrm{,}u}^A$, $I_{ma\mathrm{,}nb\mathrm{,}u}^B$, or $I_{ma\mathrm{,}nb\mathrm{,}u}^{C\mathrm{,}\text{PV}}$) are used in nine different RWG-RWG interactions. These common integrals (with singularity extractions) can be listed as

$\begin{array}{l}{I}_{ma\mathrm{,}nb\mathrm{,}u}^A=\frac{1}{A_{nb}}{\displaystyle {\int }_{S_{nb}}}dr\prime g_u(r_{ma}^{\text{cr}}\mathrm{,}r\prime )=\frac{1}{A_{nb}}{\displaystyle {\int }_{S_{nb}}}dr\prime \left[g_u(r_{ma}^{\text{cr}}\mathrm{,}r\prime )-\frac{1}{4\pi |r_{ma}^{\text{cr}}-r\prime |}\right]\\ +\frac{1}{A_{nb}}{\displaystyle {\int }_{S_{nb}}}dr\prime \frac{1}{4\pi |r_{ma}^{\text{cr}}-r\prime |}\end{array}$E33

$\begin{array}{l}{I}_{ma\mathrm{,}nb\mathrm{,}u}^B=\frac{1}{A_{nb}}{\displaystyle {\int }_{S_{nb}}}dr\prime (\mathit{\rho }\prime -{\mathit{\rho }}_{ma}^{\text{cr}})g_u(r_{ma}^{\text{cr}}\mathrm{,}r\prime )=\frac{1}{A_{nb}}{\displaystyle {\int }_{S_{nb}}}dr\prime (\mathit{\rho }\prime -{\mathit{\rho }}_{ma}^{\text{cr}})\left[g_u(r_{ma}^{\text{cr}}\mathrm{,}r\prime )-\frac{1}{4\pi |r_{ma}^{\text{cr}}-r\prime |}\right]\\ +\frac{1}{A_{nb}}{\displaystyle {\int }_{S_{nb}}}dr\prime (\mathit{\rho }\prime -{\mathit{\rho }}_{ma}^{\text{cr}})\frac{1}{4\pi |r_{ma}^{\text{cr}}-r\prime |}\end{array}$E34

$\begin{array}{l}{I}_{ma\mathrm{,}nb\mathrm{,}u}^{C\mathrm{,}\text{PV}}=\frac{1}{A_{nb}}{\displaystyle {\int }_{\text{PV}\mathrm{,}{S}_{nb}}}dr\prime {\nabla }^{\prime }{g}_u(r_{ma}^{\text{cr}}\mathrm{,}r\prime )=\frac{1}{A_{nb}}{\displaystyle {\int }_{\text{PV}\mathrm{,}{S}_{nb}}}dr\prime {\nabla }^{\prime }\left[g_u(r_{ma}^{\text{cr}}\mathrm{,}r\prime )-\frac{1}{4\pi |r_{ma}^{\text{cr}}-r\prime |}\right]\\ +\frac{1}{A_{nb}}{\displaystyle {\int }_{\text{PV}\mathrm{,}{S}_{nb}}}dr\prime \nabla \prime \left(\frac{1}{4\pi |r_{ma}^{\text{cr}}-r\prime |}\right)\mathrm{.}\end{array}$E35

Using the same convention and single-point testing, the elements of the right-hand-side vectors are evaluated as

$w_1^{\text{MCTF}}[m]\leftarrow -\frac{{\gamma }_{ma}{l}_m}{2}(r_{ma}^{\text{cr}}-r_{ma})\cdot E^{\text{inc}}(r_{ma}^{\text{cr}})$E36

$w_2^{\text{MCTF}}[m]\leftarrow -\frac{{\gamma }_{ma}{l}_m}{2}{\eta }_o{\eta }_p(r_{ma}^{\text{cr}}-r_{ma})\cdot H^{\text{inc}}(r_{ma}^{\text{cr}})$E37

for $a=\{1,2\}$.

Matrix equations obtained as summarized in this section can be solved in different ways, particularly via iterative techniques accelerated via fast algorithms. Once the current coefficients $a_J$ and $a_M$ are found, electric and magnetic fields can be obtained at any location inside or outside the object. Using the RWG functions, secondary fields can be written as

$\begin{array}{c}{E}^{\text{sec}}(r)={\displaystyle {\displaystyle \sum _{n=1}^N}}{\displaystyle {\displaystyle \sum _{b=1}^2}}{a}_J[n]\frac{{\gamma }_{nb}{l}_n}{2}i{k}_u{\eta }_u\{(\mathit{\rho }-{\mathit{\rho }}_{nb})I_{nb\mathrm{,}u}^A(r)+I_{nb\mathrm{,}u}^B(r)-\frac{2}{k_u^2}{I}_{nb\mathrm{,}u}^C(r)\}\\ -{\displaystyle {\displaystyle \sum _{n=1}^N}}{\displaystyle {\displaystyle \sum _{b=1}^2}}{a}_M[n]\frac{{\gamma }_{nb}{l}_n}{2}(r-r_{nb})\times I_{nb\mathrm{,}u}^C(r)\end{array}$E38

$\begin{array}{c}{H}^{\text{sec}}(r)={\displaystyle {\displaystyle \sum _{n=1}^N}}{\displaystyle {\displaystyle \sum _{b=1}^2}}{a}_M[n]\frac{{\gamma }_{nb}{l}_n}{2}\frac{i{k}_u}{{\eta }_u}\{(\mathit{\rho }-{\mathit{\rho }}_{nb})I_{nb\mathrm{,}u}^A(r)+I_{nb\mathrm{,}u}^B(r)-\frac{2}{k_u^2}{I}_{nb\mathrm{,}u}^C(r)\}\\ +{\displaystyle {\displaystyle \sum _{n=1}^N}}{\displaystyle {\displaystyle \sum _{b=1}^2}}{a}_J[n]\frac{{\gamma }_{nb}{l}_n}{2}(r-r_{nb})\times I_{nb\mathrm{,}u}^C(r),\end{array}$E39

where

$I_{nb\mathrm{,}u}^A(r)=\frac{1}{A_{nb}}{\displaystyle {\int }_{S_{nb}}}dr\prime g_u(r\mathrm{,}r\prime )$E40

$I_{nb\mathrm{,}u}^B(r)=\frac{1}{A_{nb}}{\displaystyle {\int }_{S_{nb}}}dr\prime (\mathit{\rho }\prime -\mathit{\rho })g_u(r\mathrm{,}r\prime )$E41

$I_{nb\mathrm{,}u}^C(r)=\frac{1}{A_{nb}}{\displaystyle {\int }_{S_{nb}}}dr\prime {\nabla }^{\prime }{g}_u(r\mathrm{,}r\prime ).$E42

Similar to the matrix elements, a triangle loop (rather than an RWG loop) can be used to efficiently perform the near-field computations. If the observation point r is close to the surface of the object, singularity extractions must be used for accurate integrations. If the medium parameters are set to ${\epsilon }_p$ and ${\mu }_p$, the computations above lead to inner electromagnetic fields, while the fields outside the surface vanish due to the equivalence theorem. In fact, this can be used to assess the accuracy of numerical solutions, since any nonzero field outside corresponds to a numerical error. Similarly, using ${\epsilon }_o$ and ${\mu }_o$, inner fields must be zero, while secondary fields are obtained outside. Then, the total fields outside the object can be obtained as

$E(r)=E^{\text{inc}}(r)+E^{\text{sec}}(r)$E43

$H(r)=H^{\text{inc}}(r)+H^{\text{sec}}(r).$E44

4. Matrix-vector multiplications with MLFMA

Plasmonic problems often involve large structures in terms of wavelength. In addition, typical ${\lambda }_o/10$ triangulations may not be sufficient to obtain accurate results, and dense discretizations are usually needed, leading to a large number of unknowns. Since direct solutions (e.g., Gaussian elimination) of the resulting matrix equations may not be feasible, fast iterative solvers are required for efficient analysis of plasmonic structures in reasonable processing times and using available memory. MLFMA is an efficient algorithm that can be used to perform fast matrix-vector multiplications with $\mathcal{O}(N\log N)$ complexity for an $N\times N$ dense matrix equation derived from an electrodynamic problem [25, 26]. Hence, MLFMA can be used within a Krylov subspace algorithm, such as the generalized minimal residual (GMRES) method, for efficient iterative solutions.

MLFMA is well known in the literature as a method with controllable accuracy. In practice, however, its accuracy heavily depends on the expansion method. In the most standard form, plane waves are used to diagonalize the addition theorem for Green’s function. Then, the interaction distances, hence, the recursive clustering of the electrodynamic interactions, are limited by a low-frequency breakdown. For example, two to three digits of accuracy (1% and 0.1% maximum relative error) using a one-box-buffer scheme need a minimum box size of around ${\lambda }_u$. It is possible to use smaller boxes and/or to achieve higher accuracy, if alternative expansion tools [38, 39], such as a direct application of multipoles [40] or evanescent waves [41], are employed. In this chapter, where numerical solutions are performed with maximum 1% error, we restrict ourselves to the plane-wave expansion.

Using plane waves, Green’s function is decomposed as

$g_u(r\mathrm{,}r\prime )=\frac{\exp (i{k}_u|r-r\prime |)}{4\pi |r-r\prime |}=\frac{\exp (i{k}_u|w+v|)}{4\pi |w+v|}\approx \frac{i{k}_u}{{(4\pi )}^2}{\displaystyle \int }{d}^2{\displaystyle \widehat{k}}\beta (k_u\mathrm{,}v){\alpha }_{\tau }(k_u\mathrm{,}w)$E45

for $w=\left|w\right|>v=\left|v\right|$, where $k_u={\displaystyle \widehat{k}}{k}_u$, $d^2{\displaystyle \widehat{k}}=d\theta d\varphi \sin \theta$, and

$\beta (k_u\mathrm{,}v)=\exp (i{k}_u\cdot v)$E46

${\alpha }_{\tau }(k_u\mathrm{,}w)={\displaystyle {\displaystyle \sum _{t=0}^{\tau }}}{(i)}^t(2t+1)h_t^{(1)}(k_{u}w)P_t({\displaystyle \widehat{k}}\cdot {\displaystyle \widehat{w}})$E47

are diagonal shift and translation operators, respectively. It is remarkable that, as a result of the factorization, the shift vector $v$ and the translation vector $w$, which satisfy $w+v=r-r\prime$, are separated. In addition, with the help of the diagonalization, sampling of the shift and translation operators leads to diagonal matrices, as the shift or translation of a plane wave in a given direction does not contribute to plane waves in other directions. In the above, the translation operator is written in terms of the Legendre polynomials P_t and the spherical Hankel function of the first kind $h_t^{(1)}$, while $\tau$ is the truncation number that can be found via excess-bandwidth formulas [42]. We note that the derivatives of Green’s function can also be obtained as

$\nabla g_u(r\mathrm{,}r\prime )\approx \frac{i{k}_u}{{(4\pi )}^2}{\displaystyle \int }{d}^2{\displaystyle \widehat{k}}(i{k}_u)\beta (k_u\mathrm{,}v){\alpha }_{\tau }(k_u\mathrm{,}w)$E48

$\nabla {\nabla }^{\prime }{g}_u(r\mathrm{,}r\prime )\approx \frac{i{k}_u}{{(4\pi )}^2}{\displaystyle \int }{d}^2{\displaystyle \widehat{k}}(k_u{k}_u)\beta (k_u\mathrm{,}v){\alpha }_{\tau }(k_u\mathrm{,}w).$E49

These expressions can directly be used to factorize the discretized operators by replacing Green’s function with the diagonalized forms. In the context of MCTF, we have

${{\displaystyle \overline{T}}}_u^T[m\mathrm{,}n\mathrm{,}a\mathrm{,}b]={\left(\frac{i{k}_u}{4\pi }\right)}^2{\displaystyle \int }{d}^2{\displaystyle \widehat{k}}{R}_{ma}^{\mathcal{T}}(k_u\mathrm{,}{r}_C)\cdot {\alpha }_{\tau }(k_u\mathrm{,}{r}_C-r_{C^{\prime }})S_{nb}(k_u\mathrm{,}{r}_{C^{\prime }})$E50

${{\displaystyle \overline{K}}}_{\text{PV}\mathrm{,}u}^T[m\mathrm{,}n\mathrm{,}a\mathrm{,}b]={\left(\frac{i{k}_u}{4\pi }\right)}^2{\displaystyle \int }{d}^2{\displaystyle \widehat{k}}{R}_{ma}^{\mathcal{K}}(k_u\mathrm{,}{r}_C)\cdot {\alpha }_{\tau }(k_u\mathrm{,}{r}_C-r_{C^{\prime }})S_{nb}(k_u\mathrm{,}{r}_{C^{\prime }}),$E51

where $r_C$ and $r_{C^{\prime }}$ are testing and basis centers, respectively. Using the RWG functions, the radiation and receiving patterns of the half basis and testing functions are derived as

$S_{nb}(k_u\mathrm{,}{r}_{C^{\prime }})=\frac{{\gamma }_{nb}{l}_n}{2}({{\displaystyle \overline{I}}}_{3\times 3}-{\displaystyle \widehat{k}}{\displaystyle \widehat{k}})\cdot \frac{1}{A_{nb}}{\displaystyle {\int }_{S_{nb}}}dr\prime \beta (k_u\mathrm{,}{r}_{C^{\prime }}-r\prime )(r\prime -r_{nb})$E52

$R_{ma}^{\mathcal{T}}(k_u\mathrm{,}{r}_C)=\frac{{\gamma }_{ma}{l}_m}{2}({{\displaystyle \overline{I}}}_{3\times 3}-{\displaystyle \widehat{k}}{\displaystyle \widehat{k}})\cdot \frac{1}{A_{ma}}{\displaystyle {\int }_{S_{ma}}}dr\beta (k_u\mathrm{,}r-r_C)(r-r_{ma})$E53

${\mathit{R}}_{ma}^{\mathcal{K}}(k_u\mathrm{,}{r}_C)=-\frac{{\gamma }_{ma}{l}_m}{2}{\displaystyle \widehat{k}}\times \frac{1}{A_{ma}}{\displaystyle {\int }_{S_{ma}}}dr\beta (k_u\mathrm{,}r-r_C)(r-r_{ma})=-{\displaystyle \widehat{k}}\times R_{ma}^{\mathcal{T}}(k_u\mathrm{,}{r}_C),$E54

where ${{\displaystyle \overline{I}}}_{3\times 3}={\displaystyle \widehat{k}}{\displaystyle \widehat{k}}+{\displaystyle \widehat{\mathit{\theta }}}{\displaystyle \widehat{\mathit{\theta }}}+{\displaystyle \widehat{\mathit{\varphi }}}{\displaystyle \widehat{\mathit{\varphi }}}$. Using a single-point integration, the patterns can be calculated as

$S_{nb}(k_u\mathrm{,}{r}_{C^{\prime }})=\frac{{\gamma }_{nb}{l}_n}{2}\beta (k_u\mathrm{,}{r}_{C^{\prime }}-r_{nb}^{\text{cr}})({{\displaystyle \overline{I}}}^{3\times 3}-{\displaystyle \widehat{k}}{\displaystyle \widehat{k}})\cdot (r_{nb}^{\text{cr}}-r_{nb})$E55

$R_{ma}^{\mathcal{T}}(k_u\mathrm{,}{r}_C)=\frac{{\gamma }_{ma}{l}_m}{2}\beta (k_u\mathrm{,}{r}_{ma}^{\text{cr}}-r_C)({{\displaystyle \overline{I}}}^{3\times 3}-{\displaystyle \widehat{k}}{\displaystyle \widehat{k}})\cdot (r_{ma}^{\text{cr}}-r_{ma})$E56

$R_{ma}^{\mathcal{K}}(k_u\mathrm{,}{r}_C)=-{\displaystyle \widehat{k}}\times R_{ma}^{\mathcal{T}}(k_u\mathrm{,}{r}_C),$E57

where $r_{ma}^{\text{cr}}$ and $r_{nb}^{\text{cr}}$ represent the centers of the associated testing and basis triangles, respectively. Then, the radiation/receiving patterns of the full RWG functions can be obtained as $S_n=S_{n1}+S_{n2}$ and $R_m^{\mathcal{K},\mathcal{T}}=R_{m1}^{\mathcal{K}\mathrm{,}\mathcal{T}}+R_{m2}^{\mathcal{K}\mathrm{,}\mathcal{T}}$ by combining the contributions of the half functions. These patterns, as well as the truncation operator, are sampled on the unit sphere, where the sampling scheme is a matter choice depending on the implementation.

In a standard implementation of MLFMA, the object is placed inside a computational cubic box, which is divided into sub-boxes until the smallest possible box size determined by the desired accuracy. Empty boxes that do not contain a part of the object (discretized surface) are omitted directly and are not divided further. This way, it is possible to construct a tree structure (consisting of $L=\mathcal{O}(\log N)$ levels) involving nonempty boxes at different levels with $\mathcal{O}(N)$ complexity. Using the child/parent relationship between the boxes, the stages of a matrix-vector multiplication, namely aggregation, translation, and disaggregation, are as follows.

In an aggregation stage, radiated fields of boxes are computed from bottom to top. At the lowest level, we have

$a[n]S_n(k_u\mathrm{,}{r}_{C^{\prime }})\to S_{C^{\prime }}(k_u\mathrm{,}{r}_{C^{\prime }}),(b_n\in C^{\prime }),$E58

where the coefficients provided by the iterative solver are used to weight the contributions of the basis functions to the overall radiation patterns of the boxes $C^{\prime }$ at the lowest level. At higher levels ($l=\mathrm{2,3,}\dots \mathrm{,}L$), aggregation is performed recursively as

$\beta (k_u\mathrm{,}{r}_{P\left\{C^{\prime }\right\}}-r_{C^{\prime }})S_{C^{\prime }}(k_u\mathrm{,}{r}_{C^{\prime }})\to S_{P\left\{C^{\prime }\right\}}(k_u\mathrm{,}{r}_{P\left\{C^{\prime }\right\}}),$E59

where $P\left\{C^{\prime }\right\}$ represents the parent of $C^{\prime }$. Due to the exponential shifts from different locations within a box, the radiated fields become more oscillatory as the box size gets larger. Hence, the sampling rate must be increased, generally with $\mathcal{O}(k_u^2{D}^2)$ where $D$ is the box size.

After completing an aggregation stage, the radiated fields are translated between the boxes at the same level. For $l=\mathrm{1,2,}\dots \mathrm{,}L$, this can be written as

${\alpha }_{\tau }(k_u\mathrm{,}{r}_C-r_{C^{\prime }})S_{C^{\prime }}(k_u\mathrm{,}{r}_{C^{\prime }})\to G_C(k_u\mathrm{,}{r}_C),(C^{\prime }\in F\{C\}),$E60

where $F\left\{C\right\}$ represents the far-zone boxes for a given box $C$. It is remarkable that F{C} contains $\mathcal{O}$(1) elements since interactions between too far boxes, for example, C and $C^{\prime }$ at level l, are made at a higher level $(l^{\prime }>l)$. Using a one-box-buffer scheme, the condition for translation is that the boxes should not intersect at a surface, line, or corner, while their parents must intersect at a surface, line, or corner.

In a translation stage, incoming fields are collected at the box centers, but they are only partial data, since the total incoming fields at the center of a box contain contribution from its parent (if exists) due to the translations at higher levels. Therefore, a disaggregation stage is performed recursively for $l=L-\mathrm{1,}L-\mathrm{2,}\dots \mathrm{,1}$ as

$G_C(k_u\mathrm{,}{r}_C)+\beta \left(k_u\mathrm{,}{r}_C-r_{P\left\{C\right\}}\right)G_{P\left\{C\right\}}^{+}\left(k_u\mathrm{,}{r}_{P\left\{C\right\}}\right)\to G_C^{+}(k_u\mathrm{,}{r}_C).$E61

At the lowest level, the testing functions receive the incoming fields as

${\left(\frac{i{k}_u}{4\pi }\right)}^2{\displaystyle \int }{d}^2{\displaystyle \widehat{k}}{R}_m^{\mathcal{K}\mathrm{,}\mathcal{T}}(k_u\mathrm{,}{r}_C)\cdot G_C^{+}(k_u\mathrm{,}{r}_C)\to y^{\text{FF}}[m],(t_m\in C),$E62

where

$y^{\text{FF}}[m]={\displaystyle {\displaystyle \sum _{n=1}^N}}{{\displaystyle \overline{Z}}}^{\text{FF}}[m\mathrm{,}n]a[n]$E63

for $m=\mathrm{1,2,}\dots \mathrm{,}N$. The overall matrix-vector multiplication is completed by also considering the near-field interactions (that cannot be calculated via aggregation-translation-disaggregation stages) as

$y[m]=y^{\text{FF}}[m]+{{\displaystyle \overline{Z}}}^{\text{NF}}[m\mathrm{,}n]a[n].$E64

Using MLFMA, each matrix-vector multiplication can be performed in $\mathcal{O}(N\log N)$ time and using $\mathcal{O}(N\log N)$ memory.

For plasmonic objects with high negative permittivity values, electromagnetic interactions decay quickly with respect to the distance between the observation and source points. For a given accuracy, interactions at long distances can be omitted since the inner and outer interactions are combined in the surface formulations and outer interactions (related to the free space) dominate the related matrix elements [30]. The threshold distance for this purpose can also be found by considering the exponential behavior of the decay for large imaginary values of the wavenumber. This way, the processing time for the matrix-vector multiplication can significantly be reduced. As the negative permittivity increases, far-zone interactions related to the inner medium may completely vanish, leaving only near-zone interactions. In the limit, near-zone interactions further reduce into self interactions of basis/testing functions, leading to the Gram matrix to represent the inner medium.

5. Case example: numerical simulations of nanowires

Using surface integral equations and MLFMA, electrical properties of a plasmonic object are simply parameters, which can be used as variables in the implementations. For the electrical properties, that is, permittivity and permeability, alternative choices, including measurement data and those based on certain models for the materials, can be used. As an example, Figure 1 presents the relative permittivity of silver (Ag) with respect to frequency from 200 to 1600 THz. In addition to measurement data [10], Drude (D) and Lorentz-Drude (LD) models are used to predict the real and imaginary parts of the relative permittivity. It can be observed that the real part of the permittivity has large negative values at the lower (infrared) frequencies and it increases smoothly toward unity as the frequency increases to the visible range and beyond. For imaginary values, which represent ohmic losses, we observe varying values between 0.01 and 10, while large discrepancies exist between measurement and D/LD models (especially considering the logarithmic scale of the y-axis). These discrepancies are responsible for different results in the simulations of plasmonic problems presented below.

Figure 1.

Real and imaginary parts of the relative permittivity of Ag with respect to frequency. In addition to measurement data [10], values based on Drude (D) and Lorentz-Drude (LD) models are depicted.

As a case study, we consider transmission though a pair of Ag nanowires described in Figure 2. The length of the nanowires is 5μm and each nanowire has 0.1×0.1μm (square) cross section. The distance between the nanowires is also 0.1μm. The transmission line is excited by a pair of Hertzian dipoles oriented in the opposite directions and located at 0.2μm distance from the nanowires. Figure 3 presents the electromagnetic response of the transmission line in the infrared frequencies from 250 to 430 THz. The power density in dB scale (dBW/m²) in the vicinity of the nanowires is depicted (normalized to 0 dB and using 40-dB dynamic range), when LD model and measurement data are used for the permittivity values. It can be observed that the electromagnetic power is effectively transmitted from the source region (right) to the transmission region (left). Coupling to the free space at the end of the line leads to two beams with decaying amplitudes due to propagation. Comparing the results, we observe very good agreement between the power density values when LD and measurement permittivity values are used. Considering Figure 1, the negative real permittivity dominates the response of the nanowires at these frequencies.

Figure 2.

A transmission line involving two Ag nanowires of length 5μm. The nanowires are excited by a pair of dipoles located at 0.2μm distance.

Figure 3.

Power density in the vicinity of the nanowire system (Figure 2) from 250 to 430 THz. Numerical results obtained by using permittivity values derived from the LD model (left column) and those based on the measurement data (right column) are compared.

Figure 4 presents similar results when the frequency is in the visible range. In this case, there are significant discrepancies between the power density values when the LD model and the measurement data for the permittivity are used. This is mainly due to the higher values for the imaginary permittivity predicted by the LD model. As the frequency increases, using the LD model, the transmission ability of the nanowire system deteriorates significantly. Specifically, the power density on the surfaces of the nanowires decreases and the transmitted power toward the left-hand side of the nanowires diminishes, leading to progressively weaker beams. It is remarkable that, using the measurement data that may be more accurate description of Ag, the transmission ability of the nanowire system is still at high levels, indicating that the transmission line operates as desired. These results may explain some of the contradictory results (especially simulations vs. measurements) for the nanowire and similar plasmonic systems investigated in the visible spectrum.

Figure 4.

Power density in the vicinity of the nanowire system (Figure 2) from 450 to 630 THz. Numerical results obtained by using permittivity values derived from the LD model (left column) and those based on the measurement data (right column) are compared.

As depicted in Figure 5, nanowires cannot maintain a good transmission ability as the frequency increases. Using the measurement data, the transmission of the nanowire system deteriorates significantly at the higher frequencies of the visible spectrum (e.g., at 770 THz). At 750 THz, the power density drops to less than −40 dB after a few μm along the nanowires. We note that the effective length of the nanowires increases with the frequency. For example, at 250 THz, the length of the nanowires is approximately $4.17{\lambda }_o$, while it is around $12.5{\lambda }_o$ at 750 THz. In addition, the effective distance between the sources and the nanowires increases. However, investigating the power values on the nanowire surfaces close to the source, it is obvious that the poor power transmission cannot be explained only with the increasing effective lengths at the higher frequencies. Since the power cannot be coupled to the free space, the power density along the nanowires possesses an oscillatory behavior. At the end of the visible spectrum, the discrepancy between the results obtained by using the LD and measurement values decreases, both predicting reduced interaction between the sources and nanowires.

Figure 5.

Power density in the vicinity of the nanowire system (Figure 2) from 650 to 830 THz. Numerical results obtained by using permittivity values derived from the LD model (left column) and those based on the measurement data (right column) are compared.

Figure 6 presents the results even at higher (lower-ultraviolet) frequencies. In this range, the nanowires are not expected to demonstrate transmission abilities, as predicted by both LD model and measurement data for the permittivity values. At lower frequencies of the range, the nanowires are more visible close to the source region, while, as the frequency increases, their effects diminish and the power distribution becomes close to that of two dipoles in free space. Figures 7 and 8 present the summary of input/output of the transmission line, for the LD model and measurement data, respectively, from 450 to 750 THz. For the input, the power density is sampled at 30 nm distance from the nanowires on a horizontal line from −1 to 1μm. The double-peak pattern due to two dipoles in opposite directions is clearly visible, with some variations due to reflections from the nanowires. For the output, samples are selected again on a horizontal line from −1 to 1μm in the transmission side at 40 μm distance from the nanowires. Using the LD model, the output pattern deteriorates significantly as the frequency increases. Using measured permittivity values, however, the double-peak pattern is effectively maintained for most frequencies until 750 THz, at which the transmission fails. Figure 9 presents the average input/output graphics, confirming consistency between the LD model and measurement data at lower and higher frequencies. On the other hand, at some frequencies in the visible range, there is more than 30 dB difference between the predicted output levels.

Figure 6.

Power density in the vicinity of the nanowire system (Figure 2) from 850 to 1000 THz. Numerical results obtained by using permittivity values derived from the LD model (left column) and those based on the measurement data (right column) are compared.

Figure 7.

Power density on the input and output sides of the nanowire system (Figure 2) from 450 to 750 THz. For the input and output, the samples are selected on 2μm lines at 30 and 40 nm distances from the nanowires. LD model is used for the permittivity values.

Figure 8.

Power density on the input and output sides of the nanowire system (Figure 2) from 450 to 750 THz. For the input and output, the samples are selected on 2μm lines at 30 and 40 nm distances from the nanowires. Measurement data are used for the permittivity values.

Figure 9.

Average input and output power density values for the nanowire system (Figure 2) from 150 to 1000 THz. Despite the consistency of the inputs, significant discrepancies in the output values obtained when LD model and measurement data for the permittivity values are used from 450 to 750 THz.

6. Concluding remarks

Surface integral equations combined with iterative algorithms employing MLFMA provide accurate solutions of nano-plasmonic problems without resorting to fundamental assumptions, such as periodicity and infinity. Three-dimensional and finite structures, which are typically of tens of wavelengths, but at the same time containing small details, can be investigated both precisely and efficiently. In addition to the visible ranges, the developed solvers are very beneficial at higher frequencies, where the discrepancy between the experimental results and theoretical predictions, such as based on the Drude and Lorentz-Drude models, increases. Surface formulations enable trivial integration of electrical parameters, allowing for fast tuning of the numerical results with the increasingly precise measurements. On the other hand, such a reliable simulation environment can be constructed only with appropriate combinations of surface integral equations, discretizations, numerical integrations, fast algorithms, and iterative techniques, as shown in this chapter. We present how to construct such an implementation with all details from formulations to iterative solutions using MLFMA, along with a set of results involving a nanowire transmission line in a wide range of frequencies to demonstrate the capabilities of the developed solvers.

Acknowledgments

This work was supported by the Scientific and Technical Research Council of Turkey (TUBITAK) under the Research Grant 113E129 and by the Turkish Academy of Sciences (TUBA) in the framework of the Young Scientist Award Program.

Authors contributions

Ö.E. developed the theory and formulations, B.K. and Ö.E. implemented the solvers. A.Ç. conducted the numerical experiments.