## 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*^{2} = −1 and *ω = 2πf* is the angular frequency.

## 2. Surface integral equations

For deriving surface formulations, we consider a plasmonic object with permittivity/permeability (

where ** E** and the magnetic field intensity

**on the closed surface (**

*H**π*), the combined operators can be written as

where

for

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, 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

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

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

## 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

Each RWG function is located on a pair of triangles sharing an edge. In the above,

while the charge neutrality is satisfied locally as

By selecting the basis and testing functions (

where

and

respectively. Furthermore, the discretized operators can be written as

where the integrals are evaluated on the supports of the testing and basis functions (

At this stage, we can consider the interaction of two half RWG functions associated with the *a*th triangle of the *m*th edge and *b*th triangle of the *n*th edge, respectively (

(24) |

where

(26) |

(27) |

where

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

where

(33) |

(34) |

(35) |

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

for

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

(38) |

(39) |

where

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

## 4. Matrix-vector multiplications with MLFMA

Plasmonic problems often involve large structures in terms of wavelength. In addition, typical

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

Using plane waves, Green’s function is decomposed as

for

are diagonal shift and translation operators, respectively. It is remarkable that, as a result of the factorization, the shift vector *P*_{t} and the spherical Hankel function of the first kind

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

where

where

where

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

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

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

where

After completing an aggregation stage, the radiated fields are translated between the boxes at the same level. For

where *F*{*C*} contains *C* and *l*, are made at a higher level

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

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

where

for

Using MLFMA, each matrix-vector multiplication can be performed in

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.

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^{2}) 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 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.

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

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

## 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.