Computational Electromagnetics in Plasmonics

1. Introduction

In electromagnetics, numerical techniques have been essential in the development of new technology in the last two decades. The rapidly growing computer capacity and calculation speeds make accurate solutions of very complex problems feasible. This has been especially true in the design of microwave and millimeter wave components and antennas. Whereas 30 years ago, the design of an antenna was based on simple analytical models, or trial and error strategies, nowadays, simulations seem to be as crucial to the design as real measurements.

The situation is quite different in plasmonics. Plasmonics is a quite novel research field and the application of computational electromagnetics in plasmonics can be categorized as “very recent”. There are many challenges that still need to be faced and “missing links” that have to be solved.

The plasmonic structures targeted are structures in the order of magnitude of a wavelength at plasmonic frequencies, i.e. at near IR and optical frequencies, and beyond. Although this frequency range is totally different from the traditional range where computational tools have been developed, i.e. the microwave range, in most cases no special modeling techniques have to be used. By far most plasmonic topologies reported in literature have been analyzed / designed with the well-known numerical techniques implemented within in-house developed or commercial software packages. This means that in this chapter the major numerical techniques can be overviewed in a general sense, referring to standard literature. These techniques will not be derived or explained here in full detail. Instead, this chapter focuses on those aspects that come into the picture when the structure is plasmonic.

After the section on techniques and tools, this chapter will focus on the performance of these techniques and tools for plasmonic structures. This is done through an overview of benchmarks available in literature and by considering a few thoroughly analyzed structures. Missing links will be pointed out and suggestions will be given for the future.

2. Fundamental physical modeling differences

This section discusses the differences between classical topologies handled with computational electromagnetics and plasmonic topologies. Classical means structures at much lower frequencies, for example at microwave frequencies, where computational electromagnetics is a well-developed mature field.

2.1. Frequency

First of all, it is essential to point out that the interaction between light and plasmonic structures in the frequency bands considered can still be analyzed with a high degree of accuracy using classical electromagnetic theory.Although the frequency is orders of magnitude higher in plasmonics compared to microwaves, Mawell’s laws are the same. They are linear, and thus scalable. As long as quantum effects do not have to be taken into account, which is still the case for the plasmonic applications considered, since the structures are not that small [1], there is no fundamental problem. The underlying formulation of Maxwell’s equations can remain unaffected.

The fact that at this small scale, no quantum effects have to be taken into account is really a crucial observation. It means that the concept of a “scatterer”, a device able to be excited by electromagnetic waves rather than particles, still works. Basically, the coupling between an electromagnetic (light-) wave and a plasmonic scatterer is thus the same as it is at microwave frequencies, and can be studied in the same way.

2.2. Volumetric currents

A consequence of the high frequencies is the small scale: the elementary building blocks of the topologies considered are at nanoscale. There are two important issues related to this. First, nanoscale fabrication technology of today is only able to generate 3D type structures (i.e. volumes). Thin 2D sheets, where the thickness of the metal is orders of magnitude smaller than the transversal dimensions of the pattern, as commonly used at microwave frequencies (for example a conducting strip or patch), are not possible. Second, for nanostructures operating at plasmonic frequencies the skin depth may be comparable with the structural dimensions, so that currents do flow over the complete volume and the device has to be described with volumetric currents. A surface current description is not sufficient.

2.3. Material characteristics

At microwave frequencies, in most cases the material properties are constant over the frequency bands considered. Also, apart from the losses, most metals behave more or less in the same way, i.e. as good conductors. At IR and optical frequencies however, most metals have very dispersive properties. Permittivities and conductivities may vary orders of magnitude over these frequency bands. The real part of the permittivity may even be negative. In Fig. 1, the permittivity of Cu and Ag are depicted, both real and imaginary part, which corresponds to conductivity. These data were obtained through experimental ellipsometry. Note that the difference between the two metals is enormous. It is evident that this strong variation has a serious impact on the behavior of a device over the frequency band. The variability of material characteristics thus has to be taken into account fully into the modeling.

Figure 1.

Top: permittivity of Cu; bottom: permittivity of Ag.

3. Modeling techniques

In this section the full wave modeling techniques are introduced and categorized on the basis of their solution method: Finite Elements (FE), Finite Differences in the Time Domain (FDTD), Finite Integration Technique (FIT), and Integral Equations (IE) solved by the Method of Moments (MoM). Based on their theoretical specificities, the application of each method in the case of plasmonics is discussed.

The cradle of computational electromagnetics can be found in the microwave research community. Since this community traditionally is dealing with structures in the order of wavelengths, right from its beginning days, it had no choice than to try to rigorously solve Maxwell’s equations. Most standard reference works on full-wave computational electromagnetics by consequence can be found within this community. A history and a comprehensive overview of the different numerical techniques and of their application in computational electromagnetics (CEM) may be found in [2]. Recent work on the last developments in CEM [3] concentrates on the two main approaches of differential and integral methods. A Good review and perspectives concerning the relationship between differential and integral equations (IE) modeling is recommended in a review paper by Miller [4]. A thorough discussion on the different techniques with clarifying examples is also given in [5].

3.1. Differential equation techniques

In time domain, Maxwell’s equations are

$$\nabla \times E=-\frac{dB}{dt}$$E1

(1)

$$\nabla \times H=\frac{dD}{dt}+J$$E2

(2)

with E and B the electric field and the magnetic induction, respectively, H and D the magnetic field and the electric induction, respectively, and J the electric current flowing. In free space, and in a homogeneous, isotropic, time-invariant, linear medium (most materials behave like this in the microwave frequency range), the following relations hold

$$D=\epsilon E$$E3

(3)

$$B=\mu H$$E4

(4)

They are called the constitutive relations, with$$\epsilon$$ and $$\mu$$ the permittivity and permeability of the medium surrounding the observation point considered. It is crucial to point out that in plasmonics in general (3) and (4) cannot be used. For example, the very dispersive properties of metals at optical frequencies (and beyond) prohibit this for these materials.

In time domain, the time variation has to be determined. In frequency domain, it is assumed that all field quantities are varying in a sinusoidal way. This means that the variation with time is known, and only the variation of the fields in space has to be determined. Using complex notation for the frequency domain with a $$\exp (j\omega t)$$ time dependency, Maxwell’s equations become

$$\nabla \times E=-j\omega B$$E5

(5)

$$\nabla \times ({\mu }_{}^{-1}B)=j\omega \epsilon E+J$$E6

(6)

with $$\omega =2\pi \cdot$$frequency the pulsation. Note that in (5) and (6) $$\epsilon$$and $$\mu$$ are used. In general, here they are depending on frequency, and in this way dispersion is fully taken into account (see also section 2.3). In the field of plasmonics, this is an important difference between computational tools in time and in frequency domain. Computational tools in frequency domain may use directly the well-known concepts of permittivity and permeability, albeit that they become frequency dependent. Computational tools in time domain cannot use these concepts directly. As will be shown later, this complicates things.

From a mathematical perspective, Maxwell’s equations are differential equations relating vector fields to each other.Differential equation methods in Computational Electromagnetics are methods that directly consider Maxwell’s equations (or theHelmholtz wave equations derived from them), with little analytical preprocessing. Basically, these differential equations, which are valid in any point of space and time,are solved by approximating them by difference equations, which are valid in a discrete set of points in space and time. This is done by chopping up space (and time if time domain is considered) in little pieces in which the field variation has a pre-described profile. This reduces the problem of fields varying over space (and time) to a discrete (matrix) problem that can be handled on a computer. The differences between the several differential equation methods are related to the different ways in which space (and time) can be chopped up.

Since the number of unknowns is proportional to the volume and the resolution considered, differential equation methods are particularly suitable for modeling small full three-dimensional volumes that have complex geometrical details, for example smaller closed-region problems involving inhomogeneous media [6]. Intrinsically, differential equations are less suited for open problems. The reason is that in principle they require a discretization of the entire space under consideration. This space is limited in case of closed problems, but corresponds to infinite space in case of open problems. In practice, this problem is solved by the introduction of techniques like Absorbing Boundary Conditions, and Perfectly Matched Layers (PML) [7]. They mimic the wave arriving at the boundary as travelling further to infinity. The quality of these truncating techniques nowadays is very high so that, in practice, the intrinsic problem with open structures has been overcome, albeit it in an approximate numerical way.

The most popular differential equation-based methods are the Finite Element Method (FEM), for example utilized in Ansoft’s HFSS software package, and the Finite-Difference Time Domain method (FDTD), which is employed for example by CST’s Time Domain transient solver (in the particular case of Cartesian grids), and by Lumerical.

3.1.1. The Finite Element Method (FEM) [7], [8]

FEM is a method based on solving partial differential equations. It is most commonly formulated based on a variational expression. It subdivides space in elements, for example tetrahedra. Fields inside these elements are expressed in terms of a number of basic functions, for example polynomials. These expressions are inserted into the functional of the equations, and the variation of the functional is made zero. This yields a matrix eigenvalue equation whose solution yields the fields at the nodes. FEM gives rise to a very sparse matrix equation, which can be solved using dedicated matrix algebra technology, leading to very fast solution times, considering the huge number of unknowns.

Its first formulations were developed as matrix methods for structural mechanics. This lead to the idea to approximate solids and Courant (1942) introduced an assembly of triangular elements and the minimum of potential energy to torsion problems [9]. The first paper on the application of FEM to electrical problems appeared in 1968 [10]. An extensive review on the history of FEM in electromagnetics was published in an issue of the Antennas and Propagation Magazine [11]. FEM normally is formulated in the frequency domain, i.e. for time-harmonic problems. This means that, as for IE-MoM, the solution has to be calculated for every frequency of interest.

Numerous references can be given developing, explaining, and using FEM. A good book to start with is [8]. A software tool using FEM and very widely spread is Ansoft HFSS.

3.1.2. The Finite-Difference Time-Domain technique (FDTD) [12], [13], [14]

The nature of Maxwell's differential equations is that the time derivative of the H-field is dependent on the curl of the E-field, and the time derivative of the H-field is dependent on the curl of the E-field. These basic properties result in the core FDTD time-stepping relation that, at any point in space, an updated value of an E/H-field in time is dependent on the stored value of the E/H-field and the numerical curl of the local distribution of the H/E-field in space. The numerical translation into a time-stepping algorithm was introduced by Yee in 1966. Indeed, swapping between E-field and H-field updates allows to define a marching-on-in-time process wherein sampled fields of the continuous electromagnetic waves under consideration are used. These waves can be seen to propagate in the Yee lattice, a numerical three-dimensional space lattice comprised of a multiplicity of Yee cells, see Fig. 2. More specifically, Yee proposed a leapfrog scheme for marching-on in time wherein the E-field and H-field updates are staggered so that E-field updates are observed midway during each time-step between successive H-field updates, and vice versa. A huge advantage of FDTD is that this explicit time-stepping scheme avoids the need to solve simultaneous equations, so no matrix inversions are necessary. It also yields dissipation-free numerical wave propagation. Negative is that this scheme results in an upper bound on the time-step to ensure numerical stability. This means that simulations may require many thousands of time-steps for completion. The use of the Yee lattice has proven to be very robust in numerical calculations.

FDTD is extremely versatile since the interaction of an electromagnetic wave with matter can be mapped into the space lattice by assigning appropriate values of permittivity to each electric field component, and permeability to each magnetic field component. This can be done without seriously compromising the speed of the method.

The fact that time is observed directly, and the versatility of the method, make it probably the most efficient technique for complex 3D transient problems. An implementation of the FDTD technique dedicated to plasmonics and photonics can be found in the commercial tool Lumerical [15].

Figure 2.

Standard Cartesian Yee cell used in the FDTD technique.

3.1.3. The Finite Integration Technique (FIT) [16]

The Finite Integration Technique was introduced by Weiland in 1977. The word integration does not imply any relation with integral equations. FIT first describes Maxwell’s equations on a grid space. The matrix equations for the electromagnetic integral quantities obtained by FIT possess some inherent properties of Maxwell's equations, for example with respect to charge and energy conservation. This makes them very attractive from a theoretical point of view. FIT can be formulated on different kinds of grids, e.g. Cartesian or general non-orthogonal ones, which is a clear advantage. In the time-domain, the resulting discrete grid equations of FIT are, at least in “some” cases, identical to the discrete equations derived with the classical Finite-Difference Time-Domain (FDTD) method. In contrast to FIT, which is applied to the integral form of the field equations, FDTD (as a subset of the finite integration method) is applied to the differential form of the governing Maxwell curl equations. Theoretical links between FDTD and the FIT approach realized in CST may be found in [17]. A comparative study on Time Domain methods was presented in [18]. In some sense, FIT can thus be considered as a powerful generalization of the FDTD technique. A software tool using FIT and very widely spread is CST Microwave Studio.

3.2. Integral Equation techniques (IE)

Integral equation methods make use of Maxwell’s equations in integral equation form to formulate the electromagnetic problem in terms of unknown currents flowing on the object to be described. These currents are induced by a field incident on the object. This incident field can be a real incident field, travelling in space, or a bounded wave feeding the object for example via a transmission line. An integral equation solution is fundamentally based on the combination of two equations.

In the case of Electric Field Integral Equations(the case of prime importance to the plasmonics community), the solution is based on considering the electric field. The first equation is

$$E^{sca}(r)={\displaystyle \underset{V\text{'}}{\int }G(r,r\text{'})}\cdot J^{ind}(r\text{'})dV\text{'}$$E7

(7)

which gives the scattered electric field $$E^{sca}$$generated in an arbitrary observation point r in space in terms of an induced volumetric electric current distribution $$J^{ind}$$ flowing within the volume $$V\text{'}$$ of the object considered. The kernel $$G(r,r\text{'})$$ is a so-called dyadic Green’s function, which is a tensor, and which relates a Dirac impulse type current flowing in $$r\text{'}$$ to the electric field it generates. The main strength of the integral equation technique is that in many cases the dyadic Green’s function can be calculated either directly analytically, such as in homogeneous media, or as some kind of inverse Fourier transform of an analytic function, as is the case for example in layered media. It is crucial to understand that the dyadic Green’s function is actually a solution of Maxwell’s equations and thus rigorously takes into account the background medium, for example a stack of dielectric layers.This means that unknown currents only have to be assumed on the objects embedded within this background medium. This results in an enormous reduction of unknowns compared to differential equation techniques, which have to model these dielectric layers in just the same way as the objects embedded within them. Basically, the problem formulation automatically covers the entire surrounding space without making any fundamental approximations. As a consequence, the corresponding solution is automatically valid in every point of the background medium. Far field radiation phenomena, surface waves in layered structures, etc., that are vital for efficient and accurate analysis, are analytically included in the solution.

In the case of plasmonics, the second equation is

$$\begin{array}{l}{J}^{ind}(r\text{'})=j\omega \text{ (}\epsilon (r\text{'})-{\epsilon }_0\text{) }{E}^{ind}(r\text{'})\\ =j\omega \text{ (}\epsilon (r\text{'})-{\epsilon }_0\text{) }(E^{sca}(r\text{'})+E^{sca}(r\text{'}))\end{array}$$E8

(8)

It expresses the boundary condition which has to be enforced within the object under consideration. It is a relation between the total field, which is the sum of the incident and the scattered field, and the volumetric electric current, which can be considered as partially being a conduction current (due to the imaginary part of the permittivity), and a polarization current (due to the real part of the permittivity).

Combining (7) and (8) yields the equation from which the currents can be solved. This can be done by chopping up the object considered (only the object, and thus not entire space !!!) in little pieces in which the current variation has a pre-described profile (a so-called basis function). This reduces the problem of currents varying over the object to a discrete (matrix) problem that can be handled on a computer. IE-MoM gives rise to a dense matrix equation, which can be solved using standard matrix algebra technology.

Next to the reduced number of unknowns, there are other theoretical advantages linked to the IE-MoM method. The second advantage is that, if properly formulated, IE-MoM is variationally stable since most of the output parameters are expressed in integral form over the equivalent currents. This means that even if the calculated currents differ considerably from the exact solution, integral parameters over both currents may remain very similar. Further, this will be illustrated by showing that even with a rather rough mesh high quality physical results may be obtained. A third advantage is that MoM does not heavily suffer from field singularities for example near sharp edges, since they are analytically incorporated inside the Green's functions. For differential equation methods special care (= a fine mesh) should be taken to describe correctly these field singularities.

The main disadvantage of IE-MoM is that, although there are many efforts in that direction, there is still a lack of matrix solvers operating on dense matrices which are comparable in efficiency to the solvers used in the differential equation techniques (mainly FEM), which yield sparse matrices. This has precluded the use of the very flexible volumetric IE-MoM technique for modeling complex structures on widespread computer systems.

IE-MoM is normally applied in the frequency domain, i.e. for time-harmonic problems. This means that the solution has to be determined at each frequency of interest.

Much more details can be found in the classic and basic Method of Moments (MoM) book [19]. There are many variants of the method. In general, boundary equations can be enforced also at the boundaries of volumes (utilizing Surface Integral Equations (SIE) [20]), next to inside the entire volumes themselves(applying Volume Integral Equations (VIE) [21] at the inside of the components, as described above). Also, the integral equations can be written down in different forms (dyadic form, mixed-potential form, hybrid forms, etc. [22]), which give rise to specific implementations. Further, the first theoretical developments in the field of computational plasmonics are appearing in literature [23].

In the following sections, several aspects linked with the use of integral equation techniques for plasmonic structures are discussed.

3.2.1. 3D volumetric current

The fact that a 3D volumetric current has to be described, as mentioned already, invokes the need for either a real volumetric MoM implementation or the more common surface approach, where the plasmonic component is described by applying the equivalence principle at its surface. The surface approach has already been used, for example in [24]. The volumetric approach as described in [21], has very recently been introduced in the plasmonics community [25]. For realistic structures, where volumes can be embedded within a layer structure, it offers a very flexible solution technique.

3.2.2. Material characteristics

Although strongly varying material characteristics may seem trivial to implement in an integral equation scheme, in some cases more severe consequences occur. For example, some Green’s function calculation schemes are based on the fact that the material properties do not change over the frequency band of interest. Any particular implementation of IE-MoM needs this essential relaxation into its formalism, which possibly needs to be adapted.

An advantage of a frequency domain technique like IE-MoM compared to a time domain technique like FDTD is that the ellipsometric measurement data can be used directly in the tool, without any further fitting. FDTD techniques developed for the optical range tend to fit the dielectric response of the metal to the experimentally determined dielectric permittivity using a Drude model [26] or more sophisticated models (for example in Lumerical). However, this may create problems (see further).

3.2.3. Occurrence of layer structures

As explained above, IE technique solvers are formulated making use of Green’s functions. These Green’s functions can be formulated for multi-layered structures, such that the background medium of the structure may consist of an arbitrary number of horizontal, infinitely stretched, dielectric and metallic layers, which are taken into account analytically. This environment is particularly interesting for plasmonic structures, because they are traditionally deposited on a flat layered dielectric substrate, for example glass. Thisglass can be contained in the background environment. The only remaining components are local scattering objectswith medium or small dimensions compared to the wavelength, such as the dipole in section 5.2.1, but also dots, rods, monomers, dimmers, rings, discs, etc.. With volumetric integral equations these components are replaced by equivalent volume currents, which appear as the primary unknowns in the resulting integral equations.Since the substrates used in normal circumstances are huge compared to the wavelength, edge effects are extremely small and can be neglected. It may thus be concluded that plasmonic structures have specific features which make a IE-MoM solution quite attractive.

4. Software tools

In the following sections, several solvers are briefly described. Several commercial solvers and one academic solver are considered. It has to be emphasized that the author does not claim that this overview is complete. Since the widespread use of computational electromagnetics in plasmonics is quite recent, it is highly probably that there are more solvers that the author is not aware of.

4.1. Commercial software tools

4.1.1. HFSS [27]: FEM

Since it was one of the first tools in the market, and also due to its generality and flexibility, HFFS is one of the tools heavily used in industrial microwave and millimeter wave design environments. The purpose of HFSS is to extract parasitic parameters (S, Y, Z), visualize 3D electromagnetic fields (near- and far-field), and generate SPICE models, all based on a 3D FEM solution of the electromagnetic topology under consideration. This software is extremely popular and is used for all kinds of purposes. Numerous results for plasmonic topologies can be easily found by googling the words HFSS and plasmonics.

The first step in a new simulation project is to define the geometry of the structure. The geometry is modeled in the GUI or can be imported from another program (AutoCAD, STEP, …). Then the user continues by defining material properties, boundary conditions, and excitations to the different domains and surfaces of the geometry. If desired one can control the meshing process manually. Finally, the solution parameters are defined. The most important are the solution frequency, the order of the base functions and the matrix solver type. The direct matrix solver is default, but an iterative solver can be used as well. Very useful is the automatic adaptive mesh generation and refinement, which in many cases frees the designer of worrying about which mesh/grid to choose. First a solution for an initial mesh is determined and the solution is assessed. If the solution does not qualify, the mesh is refined and a new solution is computed. This procedure is repeated until one of the exit criteria is fulfilled. It is important to note that HFSS offers no curved elements for a better approximation of curved objects. This means that the mesh along a curved surface has to be chosen very fine. This can be achieved with the 'surface approximation' parameter. Finally the solution, i.e. electric and magnetic fields, currents, and S-parameters can be visualized in 1D, 2D and 3D. All solutions can be exported as files. An 'Optimetrics' toolbox is offered for optimizations, parameter, sensitivity and statistical analysis.

4.1.2. COMSOL multiphysics [28]: FEM

Comsol Multiphysics uses the finite element method to solve partial differential equations in 2D and 3D. It is a multiphysics code, which means that it can handle not only electromagnetics, but also acoustics, mechanics, fluid dynamics, heat transfer, etc.. The geometry, material parameters and boundary conditions have to be set up in the graphical user interface (GUI). Then, after defining boundaries and domains the mesh is generated. This includes the selection of the order of the basis functions and the curved mesh elements (order higher than one means curved element). They are used for a better approximation of curved boundaries. The meshing can be steered by specifying a number of nodes and their distribution on each edge of the model. Comsol Multiphysics offers a wide variety of matrix solvers to solve the system of equations. There are direct solvers, more suitable for smaller problems, and iterative solvers suitable for larger problems. Comsol Multiphysics can generate 2D and 3D plots. Quantities derived from the electromagnetic field, like Poynting vector, energy, and energy loss, are available as predefined variables. Comsol Multiphysics does not allow parameter sweeps or optimization. COMSOL is very easily combined with MATLAB. A simulation set up with the GUI can be saved as a Matlab script, which can be easily modified. The Matlab optimization toolbox can be used to optimize a topology.

4.1.3. JCMsuite [29]: FEM

JCMsuite is a software package based on FEM. It contains dedicated tools for certain simulation problems appearing in nano-optics and plasmonics. It incorporates scattering tools, propagation mode tools, and resonance mode tools.The time-harmonic problems can be formulated in 1D, 2D and 3D.The topologies can be isolated or may occur in periodic patterns, or a mixture of both. Interesting is that dedicated tools for problems posed on cylindrically symmetric geometries are available. The electromagnetic fields are discretized withhigher order edge elements. JCMsuite contains an automatic mesh generator (with adaptive features), goal-oriented errorestimators for adaptive grid refinement, domain-decomposition techniques and fast solvers.

4.1.4. CST [30]: FIT

CST Microwave Studio (CST MWS) is based on the finite integration technique (FIT). It allows to choose the time domain as well as the frequency domain approach. Despite the presence of transient, eigenmode, and frequency domain solvers within CST MWS, the transient solver is considered as the flag ship module. The Time Domain Solver calculates the broadband behavior of electromagnetic devices in one simulation run with an arbitrarily fine frequency resolution. The modeling of curved structures using the Perfect Boundary Approximation^® technique and the modeling of thin perfectly electric conducting sheets with the Thin Sheet Technique^® tries to cope with the typical difficulties inherent to FDTD methods for classical structures. Several mesh types can be chosen. The automatic mesh generator detects the important points inside the structure (fixpoints) and locates mesh nodes there. The user can manually add fixpoints on a structure, as well as fully control the number of mesh lines in each coordinate with regards to the specified wavelength. Energy based adaptation of the mesh allows to refine it in a predefined number of passes, providing a mesh refinement of sophisticated design features for the price of longer overall simulation times. Although the FIT in principle can handle material parameters changing over the dielectric volumes defined, this is not implemented yet.

CST, as a general purpose software package being a real competitor for HFSS in the traditional microwave field, has gained a lot of popularity in the last few years. Also for the analysis and design of plasmonic structures, more and more results obtained with CST can be found in literature (just google the words CST and plasmonics). A problem sometimes observed with CST is a ripple in the frequency response in case the tool settings are not appropriate. This is due to the fact that the flagship of CST is inherently a time domain solver.

4.1.5. Lumerical [15]: FDTD

Lumerical is the leading software tool in the plasmonics and photonics community. It is based on the FDTD algorithm, and can be used for 2D and 3D topologies. The topology has to be generated in the GUI or can be imported with GDSII/SEM files. Basic shapes include triangles, rectangular blocks, cylinders, conic surfaces, polygons, rings, user-defined (parametric) surfaces, spheres and pyramids. Several boundary conditions can be used: absorbing (PML), periodic, Bloch, symmetric, asymmetric, and metal boundaries. A non-uniform mesh can be used and automesh algorithms are provided. Since it is a dedicated time domain solver, it has a sophisticated library of Lorentz, Drude, Debye models for the material parameters. Excitation of the structure is possible with waveguide sources, dipoles, plane waves, focused beams and diffraction-limited spots. A scripting language is available to customize simulation. Data can be exported to Matlab or in ASCII.

4.2. Non-commercial software tools

4.2.1. MAGMAS 3D [31]: IE-MoM

MAGMAS 3D is the IE-MoM code developed at the Katholieke Universiteit Leuven, Belgium. It was originally developed in cooperation with the European Space Agency for planar and quasi-planar antenna and scattering structures embedded in a multilayered dielectric background medium and operating in the microwave frequency range. Starting from the topology considered, frequency band needed, and type of excitation, it calculates the network, radiation, and scattering characteristics of the structure under consideration. It is based on a full-wave Mixed-Potential formulation of Electric Field Integral Equations (MP-EFIE), originally applied only to 2D surface currents [22], [32]. New theoretical techniques were developed and implemented within the framework: the Expansion Wave Concept (EWC) [33], [34], the Dipole Modeling Technique (DMT) [35], special de-embedding procedures [36], etc..In 2007, it was extended with the capability to handle volumetric 3D currents with the VIE technique introduced in [21]. This was crucial in view of its application in plasmonics. In 2009, this led to the first verified simulated results in this field [25]. Now, it is extensively used for plasmonics. The MAGMAS mesh is based on a combination of rectangular and triangular mesh cells.Full mesh control is available in manual meshing mode. A Graphical User Interface is available.

To the best knowledge of the author, to date (June 2012), it is the only IE-MoM framework able to handle arbitrary plasmonic structures embedded in a multilayered environment, with the very general and flexible VIE technique. No similar commercial solvers exist, in contrast to the situation at microwave frequencies, where a multitude of IE solvers can be found.

4.2.2.Other. non-commercial tools found in literature

At this moment plasmonics is a mainly experimentally driven research field. Whereas for example in the field of traditional antenna research, numerous papers can be found on modeling as such, this is much less the case in plasmonics, especially in the high impact journals. Most papers are concerned with the description of the physical phenomena occurring and occasionally just mention the numerical tool used.

Most researchers in the plasmonics field use the commercial tools available, and as far as I can see, the main one is Lumerical. In the publications that do treat the numerical modeling of plasmonic structures as such, mainly the FDTD or the FEM technique is applied [37], [38]. Very few papers consider the IE-MoM technique, and if they do, it is the Surface IE technique [39], [40], [41]. In [42] Chremmos uses a magnetic type scalar integral equation to describe surface plasmon scattering by rectangular dielectric channel discontinuities. Even the use of the magnetic current formalism to describe holes, classical at microwave frequencies, already has been used in plasmonics [43]. However, these dedicated developments cannot be categorized under the title „analysis framework“ in the sense that they do not allow to handle a wide range of different topologies.

5. Benchmarking

Nowadays, physicists and engineers rely heavily on highly specialized full wave electromagnetic field solvers to analyze, develop, and optimize their designs. Computer-aided analysis and optimization have replaced the design process of iterative experimental modifications of the initial design. It is evident that the underlying solution method for a software tool may significantly influence the efficiency and accuracy by which certain structure types are analyzed. Nevertheless, the commercial focus increasingly switches from such key theoretical considerations to improvements in the area of layout tools and system-level design tools. Therefore, users may get the wrong impression that a given solver is automatically suited to solve any kind of problem with arbitrary precision. This is of course not true.

This section verifies the plausibility of such expectations by presenting a benchmark study for a few plasmonic structures. The study focuses on the capabilities and limitations of the applied EM modeling techniques that usually remain hidden for the user.

5.1. Benchmarking in literature

In literature, not many comparisons between solvers in the field of plasmonics can be found. In [44] Hoffmann et al. consider a single plasmonic topology, a pair of Au spheres, and analyze this structure with several electromagnetic field solvers: COMSOL Multiphysics (FEM), JCMsuite (FEM), HFSS (FEM), and CST Microwave Studio (the FEM solver available in the Studio is used, not the flagship FDTD solver). The output parameter considered is the electric field strength in between the two spheres. Note that all solvers tested are based on the Finite Elements technique in the frequency domain, so it can be expected that the main issues with these solvers for this topology are the same. Very interesting is the fact that in the paper itself CST is categorized as “inaccurate”. However, afterwards this was corrected. It seems that wrong material settings were used, and after correction accurate results were obtained. This issue led to the fact that this benchmark example is now presented as a reference example on the CST website.

In [45], several numerical methods are tested for 2D plasmonic nanowire structures: not only the Finite Element Method (FEM) and the Finite Difference Time-Domain (FDTD) technique, but also less „commercial“ methods like the Multiple Multipole Program (MMP), the Method of Auxiliary Sources (MAS), and the Mesh-less Boundary Integral Equation (BIE) method are tested. By comparing the results, several conclusions can be drawn about their applicability and accuracy for plasmonic topologies. Differential techniques like FEM and FDTD can reach a high level of accuracy only with a high discretization. In 2D, this is readily affordable on present-day computer systems. There, these techniques have a clear advantage in terms of speed, matrix size, and accuracy. In 3D, due to the rocketing size of the problem, this may not always be that obvious. The advantages may thus disappear in a full 3D analysis of geometrically complicated structures. This happens because matrices become denser or more ill-conditioned. The main conclusion of the paper is that the most efficient method depends on the problem dimension and complexity. This is a similar conclusion as was reached in [46] within the context of the analysis and design of planar antennas.

5.2.Comparison between differential and integral equation techniques

In this section, the differential and integral equation techniques are compared. This is done by choosing a representative solver from each category and using it in the analysis of a basic plasmonic topology. In the category differential techniques, Lumerical is chosen, as it is the most widespread commercial solver in the plasmonic community. Since, as far as we know, there are no commercial integral equation solvers in the plasmonics area, MAGMAS 3D is chosen in this category, the in-house developed solver at Katholieke Universiteit Leuven. The basic topology selected is the plasmonic dipole. It is a structure that can function as a scatterer, but also as a real nano-antenna. The main result of this section is that it draws conclusions on different computational aspects involved in the modeling of nanostructures with the two solvers.

5.2.1. Plasmonic dipole

The plasmonic dipole is depicted in Fig. 3. It is a structure consisting of two nanorods with a gap in between. When the gap is “shortcircuited”, it functions as a single rod, when the gap is open, it functions as a device generating an enhanced electric field there. When it is connected to other nanocircuits, it may function as a nano-antenna, both in transmit (Fig. 3b), and in receive (Fig. 3c). The two rods may be fabricated from metals like Au, Ag, Cu, Al, Cr, etc.. In this section the width W and height H of the dipole are kept constant at 40 nm. This value is well-chosen, since it is a value which can be fabricated with sufficient accuracy using present-day nanofabrication technology. The gap width can have different values, depending on the case considered. Note however that 10 nm is about the minimum that can be fabricated with reasonable accuracy nowadays.

5.2.1.1. Quality of the input

It cannot be over-emphasized that Computational Electromagnetics solvers in general produce correct results consistently for a multitude of different structures only if two conditions are met:

1. The user has sufficient general background in the field, and a more specific knowledge about the solution technique implemented in the solver that he is using. If not, there is a huge danger that the solver is not used in a proper way. For example, although this is not necessary to “operate” the solver, if the user has no basic knowledge and does not realize how the discretization / meshing scheme works, it is impossible to understand and assess the effect of a proper meshing. In the majority of the cases, there will be no problem with that, and automatic meshing schemes will be able to produce good results. However, in unusual and/or challenging cases, often met in scientific research, this issue will be a crucial factor.

2. The user needs to have sufficient knowledge about the actual tool that he is using and its peculiarities, more specific about the way a certain problem has to be imported into the solver. A very good illustration of this has already been addressed in section 5.1 where CST was actually not properly used and produced incorrect results that were published in open literature. Afterwards, this was corrected and the problem disappeared. However, following the proper line of reasoning is not always straightforward. Even experienced researchers sometimes easily can make mistakes. It is the opinion of the author that this is happening too much, even in a lot of peer-reviewed scientific papers.

Figure 3.

The dipole model studied. a) W and H are set equal to 40 nm and the gap G can vary. b) The dipole as transmitting antenna with a model for the feeding structure located in the gap. c) The dipole as receiving antenna excited by a plane wave.

5.2.1.2. Modeling of the feed

The first parameter studied is the impedance that is seen when the nanodevice operates as a transmitting or receiving antenna. In Fig. 4, the impedance simulated with both MAGMAS and Lumerical is given for an Al dipole of 200 nm length and with a gap of 10 nm. The input impedance is calculated as Z = V/I, where Z, V and I are the input impedance, the voltage over the gap G, and the source current, respectively. The exciting source is located in the middle of the gap. MAGMAS uses a physical current filament of finite width w. Lumerical does not have a built-in physically feasible current source model. It uses a (less realistic) pure electrical dipole source, which means that the current I has to be evaluated “manually” by the user based on the magnetic field distribution and Ampere’s law. This also means that the impedance will be somewhat depending on the resolution chosen during meshing. It is seen that there is an excellent agreement in the trend of the curves, i.e. the change of input impedance with frequency for different widths of the source current filament in MAGMAS (5 nm and 13.33 nm) is the same as the one in Lumerical. However, in general there is an off-set. Only for a specific width of the current filament in MAGMAS, nl. 5 nm, the two excitation mechanisms become almost completely equivalent and even the offset between the two solvers disappears. A similar agreement was also observed for both gold and silver dipoles.

The main conclusion of this section is that there is still an issue with the modeling of (localized) feeds, especially in Lumerical. The reason is that the plasmonics community is used to consider objects as scatterers, excited by an incident wave. In this case, the models implemented are satisfactory. However, in view of designing nanocircuits, localized feeds and impedances derived from them become important. The basic rule is that the feeding model has to correspond as closely as possible to the actual feed topology that is going to be used later in practice. In this respect, the localized feeding models available in Lumerical today are unsatisfactory.

Figure 4.

Comparison of impedances obtained with MAGMAS and Lumerical for an aluminum dipole on a glass substrate, L = 200 nm and gap = 10 nm. Two widths w are used in MAGMAS for the source current filament. R and X stand for real and imaginary part.

5.2.1.3. Modeling of losses

The importance of taking into account losses is obvious. Metals at plasmonic frequencies may be very lossy and neglecting this would result in totally erroneous results. However, there is a difference between time domain and frequency domain solvers. Since Lumerical is an FDTD based tool, it works in time domain. This means that it cannot work directly with (complex) permittivities and permeabilities. The material data have to be transformed from functions of frequency into functions of time. This is a very difficult process, involving a convolution, and is subject to various constraints. Local curve fitting of permittivities is very difficult to implement in a time domain solver. Lumerical fits the sampled measurement data with its own multi-coefficient material model, which “provides a superior fit compared to the standard material models (such as Lorentz, Debye, Drude, etc )”. However, this produces artefacts, as is illustrated in Fig. 5. There, the relative error is given between the measured imaginary part of the permittivity (which represents losses) and the value obtained for this imaginary part after fitting, for three metals. It is seen that for gold and aluminum, the error is in the order of 20 %, which is already quite elevated. However, in the case of silver, the error reaches a full 100 %, which is unacceptable. In some cases, it will result in a totally wrong prediction for example of the radiation efficiency of nano-antennas [47]. If one forces the material characteristics to be exactly the same in the two solvers, for example by using the interpolated values also in MAGMAS, instead of the real measured ones, the two solvers produce very similar results, see Fig. 6.

The main conclusion is that a highly accurate prediction of losses in Lumerical is still an issue. The material models have to be further refined, especially for silver.

Figure 5.

Relative error between interpolated and measured imaginary part of permittivity for silver, gold and aluminum.

Figure 6.

Comparison between FDTD and IE-MoM for a gold dipole in free space, L = 250 nm and the gap is 10 nm. Identical material parameters are used in both solvers.

5.2.1.4. Meshing

It is well known that the mesh quality and resolution are key factors in the accuracy of any solver. The electromagnetic coupling between nearby segments may differ considerably due to the specific meshing used, especially in parts of the structure where rapid topological changes occur. The question is whether the meshes used in the solvers are adequate. This issue was investigated by performinga convergence study in terms of the resolution of the meshes [48].

The extinction cross section of a gold monomer (no gap, or in other words G = 0 nm) calculated using MAGMAS 3D with different meshes is plotted in Fig. 7(a). The analysis of these data demonstrates clearly the stability of IE-MoM. The results obtained even with a very rough mesh 6x1x1 (41.7 nm x 40 nm x 40 nm) provide already a very good estimation of the antenna resonance properties. The differences in extinction cross section calculated with the rough and fine meshes are almost negligibly small. This result is obtained thanks to the variational stability of IE-MoM. On the adjacent figure Fig. 7(b) extinction cross sections are calculated using Lumerical. As expected, in general the results obtained with the FDTD method depend considerably on the chosen mesh. The calculated wavelengths as a function of the mesh cell size are plotted in Fig. 7(c). It should be also noted that in contrast to Fig. 7(a) (IE-MoM) in Fig. 7(b) (FDTD) not only peak positions but also their levels depend clearly on the mesh.

The conclusion is that in Lumerical a fine mesh is mandatory for reliable calculations of the monomer resonant wavelength. Thanks to its variational approach, in IE-MoM, a highly dense mesh is not always needed, especially when scattering problems are studied.

Figure 7.

a and b: extinction cross section of a 250 nm x 40 nm x 40 nm gold monomer on a substrate with n=1.5 (a. MAGMAS 3D, b. Lumerical), c. convergence performance of both solvers.

5.2.1.5. Calculation speed

Information on the calculation times for the convergence study of the previous section is given in Table 1. It has to be emphasized that calculation times cannot be directly compared. The computers used were different, and the FDTD tool used 10 processors in parallel. Further, since it is a time domain technique, a calculation time per frequency point cannot be given for FDTD. Nevertheless, combining the data in this table with the data concerning the convergence, given in Fig. 7, it is easily seen that IE-MoM outperforms FDTD for this particular case.

5.2.2. Plasmonic nanojets

In this section it is illustrated what can be reached with computational electromagnetics in the field of plasmonics. The results in this section are extracted from the paper [49], where full details can be found.

IE-MoM		FDTD
Mesh	Calculation time per frequency point (s)	Mesh	Total calculation time (s); (10 processors in parallel)
6 x 6 x 1	5.4	6 x 1 x 1	38
12 x 2 x 2	10.4	12 x 2 x 2	55
25 x 4 x 4	33.6	25 x 4 x 4	225
50 x 8 x 6	238	50 x 8 x 8	672
		100 x 16 x 16	5893
		250 x 40 x 40	20100

Table 1.

Calculation times as a function of the mesh for the gold monomer.

Similarly to a pebble hitting a water surface, the energy deposition by femtosecond laser pulses on a gold surface can produce local melting and back-jet. Contrary to water though, a gold surface can be nanopatterned and the excitation of surface plasmon resonances leads to the appearance of hotspots, literally. This was explicitly proven for the first time in a nanopatterned gold surface, composed of a G shaped periodic structure.It was observed there that laser-induced melting and back-jet occur precisely in the plasmonic hotspots. The aid of computational electromagnetics to rigorously analyze the structure by predicting the hotspots was crucial in explaining this phenomenon. Much more details can be found in [49]. It is clear that this type of insight into the basic interaction mechanisms of light with matter at the nanoscale opens up new possibilities for applications.

Another example of what can be reached by applying computational electromagnetics in plasmonics can be found in [50].

Figure 8.

Illuminating the sample with femtosecond laser pulses produces a polarizationdependentpattern of nanobumps with sharp tips on the sample surface [49].Left: horizontally polarized light, right: vertically polarized light.

Figure 9.

The pattern of melted nanobumps matches that of the plasmonic currents, where Ohmic losses locally increase the temperature [49]. The electric field of light incident on the nanostructures causes the surface charges to oscillate in response to the direction of light polarization: horizontal, in (a), and vertical, in (b). The locations of the resulting electric currents are shown in (c) and (d) respectively. The current maxima are indicated with colored circles matching the pattern of nanobumps in Fig. 8.

6. Conclusions

In this chapter, abrief introductory overview has been given on the use of existing computational techniques and software tools for the analysis and design of plasmonic structures. This field is booming, and many modeling techniques developed at lower frequencies, i.e. in the microwave range, are now being transferred to the plasmonics community, of course with the necessary adaptations. Although the differential equation techniques are the most widespread, both in scientific literature, and in the commercial scene, it is proven that integral equation techniques are a valid alternative. In some cases they are even superior. The author sees a lot of opportunities in this area.

In most cases, the tools are used to analyze structures. Few papers really make dedicated designs. However, this will change in the future. Inevitably, the plasmonics community will follow a similar patch as the traditional antenna community. Whereas 30 years ago, the design of an antenna was based on the accumulated expertise and know-how of many years, nowadays most antennas are designed almost in a single pass by experienced designers using commercial tools. This is possible since the tools in this field have sufficient accuracy and matureness in order to be able to do that. It is my belief that as soon as plasmonics will make the unavoidable shift from the physics to the engineering community, this last community will aim at conceiving applications which will involve designs heavily based on computational electromagnetics.

I would like to end this chapter with the following guideline. The use of two different solvers, based on different theoretical methods (integral and differential) may provide an excellent means to characterize the quality of simulation results. If the two results are in good agreement, it is highly likely that the results are correct. If the two results are in disagreement, a deeper investigation of the structure and its modeling is absolutely necessary.

Acknowledgement

The author gratefully acknowledges the following persons: Zhongkun Ma for the impedance and radiation efficiency comparison between MAGMAS and Lumerical, Dr. Vladimir Volskifor the convergence study with MAGMAS and Lumerical, and V. K. Valev for making available the plasmonic nanojet figures.

We also would like to express our gratitude to the KU Leuven Methusalem project and its beneficiary Prof. V. V. Moshchalkov, whose activities created a real stimulus to become personally active in this field, and the fund for scientific research Flanders (FWO-V) for the financial support.