Understanding the Physical Behavior of Plasmonic Antennas Through Computational Electromagnetics

This chapter focuses on understanding the electromagnetic response of nanoscopic metallic antennas through a classical computational electromagnetic algorithm: volumetric method of moments (V-MoMs). Under the assumption that metals only respond to external electromagnetic disturbancelocally, we rigorously formulate the light-nanoantenna interaction in terms of a volume integral equation (VIE) and solve the equation by using the method of moments algorithm. Modes of a nanoantenna, as the excitation independent solution to the volume integral equation (VIE), are introduced to resolve the antenna ’ s complex optical spectrum. Group representation theory is then employed to reveal how the symmetry of a nanoantenna defines the modes ’ properties and determines the antenna ’ s optical response. Through such a treatment, a set of tools that can systematically treat the interaction of light with a nanoantenna is developed, paving the road for future nanoantenna design.


Introduction
In the past decade, the grand leap in nanotechnology has granted us the very capability to reproduce every single object in the macroscopic world nanoscopically. Since classical antennas have been so optimized for controlling the low-frequency (normally below 1 THz) electromagnetic interaction, their topology naturally becomes the first option to be explored at nanoscale. The nanoscale version of classical antennas, which is always referred to as a "nanoantenna," has been proven to provide an effective route to couple photons in and out of nanoscale volumes. They are considered as excellent tools to probe and control light-matter interaction at the nanoscale (e.g., the interaction between light and molecules) and therefore have become an essential element in the discipline of nanophotonics [1,2]. Till now, the concept of nanoantennas has been applied to many different fields [3][4][5][6][7][8][9][10][11] covering single-molecule detection, magnetic recording, bio-imaging, photochemistry, nanoscale signal processing, optical functional material, and so on.
The general physics of the interaction between electromagnetic fields and the collective oscillations of free electrons in a metal can be described by Maxwell's equations. This physical observation enables the reuse of traditional numerical techniques developed within the microwave community. On the one hand, well packaged in commercial solvers, differential equation-based techniques, the finite difference time domain (FDTD) method [12], and the finite element method (FEM) [13], to name a few, are widely employed by experimental physicists. On the other hand, integral equation-based techniques [14], where the light-nanoantenna interaction is mathematically depicted by surface integral equations (SIEs) or volume integral equations (VIEs), are also of essential importance. Since an integral equation involves a precalculated Green's function (analytically known) and only targets the scatterers, that is, the nanoantennas, it can provide an unbeatable accuracy, efficiency, and physical clarity compared to their differential equation counterpart.
To reach an integral equation formalism, one may start from the equivalence principle, match the electromagnetic field at the nanoantenna's boundary, and obtain surface integral equations. Another possibility is to use the volume equivalence principle. The nanoparticles are replaced by equivalent-induced electric polarization currents. With discretized boundaries or bodies, both the surface integral equation and the volume integral equation can be subsequently solved by the method of moments (MoMs) algorithm [15].
In this chapter, we solely focus on the volume integral equation formalism. A review of the surface integral equation formalism will be presented separately in the future. This chapter is organized as follows. In Section 2, under the local response material assumption, we formulate the interaction between light and a nanoantenna, which is of an arbitrary shape immersed in a generic environment composed of several dielectric regions, in terms of a volume integral equation. Then, the volume integral equation is solved by the method of moments algorithm. To validate our implementation and demonstrate its diversity, simulation results are contrasted with a wide range of experimental data.
In Section 3, the notions of eigenmodes and natural modes of a nanoantenna are presented. A mode is a fixed spatial distribution of some physical quantities (e.g., charge, current, electric field, magnetic field, etc.), and is only determined by the material, the geometry, and the environment of the structure. It is thus independent of the incident field. Two kinds of modes are introduced: eigenmodes and natural modes. While the eigenmodes are still frequency dependent, the natural modes are frequency independent. The latter represent the most fundamental properties of a structure. Different from a closed system (e.g., a cavity), the modes of an open system (i.e., a nanoantenna) are not necessarily orthogonal in an inner product sense. This implies a possible interference between the modes which is the very cause of the formation of an asymmetric line, that is, a Fano resonance, in a nanoantenna's scattering spectrum. Section 3 indicates that the modes of a nanoantenna are not necessarily orthogonal. Nevertheless, the condition under which the modes are orthogonal or not orthogonal is not specified. In Section 4, we supply such a condition by using symmetry arguments. The relation between the symmetries of a nanoantenna and the mode orthogonality is rigorously put into the mathematical framework of the group representation theory. It is found that the modes that belong to different irreducible representations are orthogonal to each other in an inner product sense, and therefore no interference between these modes is allowed.

Volume integral equation formulation and volumetric method of moments algorithm for light-nanoantenna interaction
Assume a nanoantenna in space (see Figure 1) shone upon by an incident plane wave. This plane wave is generated by a current source which is very far from the scatterer. It oscillates with a time dependency e jωt with an angular frequency ω. A general field relation, which is the total field, is the sum of the scattered field and the incident field, must be satisfied at every point in the space, Especially, at the nanoantenna the total electric field can be linked with the induced current which includes both the conduction currents due to the motion of free electrons and the displacement currents due to the motion of bound electrons through a dielectric function, In Eq. (2), ε 0 is the vacuum permittivity. However, in general it can be understood as the ambient background permittivity where the nanoantenna is positioned. ε r, ω ð Þis a function of Figure 1. An abstract model for the light interaction with a nanoantenna (scatterer). In the figure, the scatterer is composed of a material characterized by electric and magnetic parameters, that is, ε V r, ω ð Þand μ 0 . Here, the material is assumed to be non-magnetic. For the sake of conciseness, a homogeneous background is assumed and plotted but the background can be inhomogeneous. Lastly, the incident field is generated by a predefined current distribution. observation point and frequency and represents the permittivity of the material constituting the nanoantenna. Here, the locality (i.e., the material property does not depend on the source point and hence no non-local behavior is assumed), the homogeneity (i.e., within the nanoantenna, the permittivity is a constant), and the isotropy (i.e., the dielectric function is a scalar) of the material are assumed. The relation stated in Eq. (2) is the so-called volume equivalence principle.
Further, it should be realized that the scattered field is induced by the induced current (as in Eq. (2)), In Eq. (3), μ 0 is the vacuum permeability and, in the most general case, should be understood as the background permeability. G r, r 0 , ω ð Þis a Dyadic Green's function describing how an elementary-induced current J ind r 0 , ω ð Þ at a source point r 0 generates an electric field at an observation point r: The volume integral takes into account all the induced current spanning the nanoantenna and gives the scattered field at an observation point r.
Focusing on the volume of the nanoantenna, substituting Eqs. (2) and (3) into Eq. (1) and putting all the terms containing the induced current on the left side of the equation, it is readily found that Eq. (4) is the main equation in this chapter. Since the integral is taken with respect to the volume of the nanoantenna, in this manner, we formulate the interaction between light and a nanoscopic scatterer in terms of a volume integral equation. In short, Eq. (4) can be recast in an operator form, The Z operator is named as the impedance operator since this operator links current density with electric field and has the same unit as an impedance. This operator consists of two parts, The first term in Eq. (6) is called the total field operator that links the induced current with the total electric field, The second term in Eq. (6) is called the scattered field operator that links the induced current with the scattered field, By inverting the impedance operator in Eq. (5), in principle the induced current J ind r, ω ð Þ can be calculated. The scattered field in the space other than the region of the nanoantenna can be further calculated by using Eq. (3), and therefore the total field at any point in space can be recovered.
Since analytical solutions for the main equation (4) only exist for a few cases (e.g., the Mie solution in case of a sphere), in general the equation has to be solved numerically. To this end, a method of moments algorithm can be employed [15]. The algorithm starts with the discretization of the antenna with tetrahedral or hexahedral blocks (see Figure 2) and further assumes that the current density flowing from one tetrahedral (hexahedral) element to a neighboring one forms a certain basis function ϕ i r ð Þ, for example, a rooftop shape [16,17].
Approximating the induced current J ind r, ω ð Þby a weighted sum of the basis functions, and inserting Eq. (9) into Eq. (5), one equation with n unknowns is found In order to construct a system of N equations with N unknowns, a testing procedure is applied.
Assuming that a test function is ψ j r ð Þ, it is obtained that In Eq. (11), z ji È É , w i f g, and e j È É represent the impedance matrix, the vector containing the weighting coefficients as in Eq. (9), and the vector containing the projection of the incident field onto the jth basis function. z ji È É and e j È É are Figure 2. Continuum description (a), discretized description (b), and tetrahedral/hexahedral blocks (c) used to mesh the continuous volume.
Notice that in Eqs. (12) and (13), the integrations are with respect to the volume of the ith basis function, V i , and the jth test function, V j . By inverting the impedance matrix, the induced current can be reconstructed and the related physical quantities, for example, the near field and the far field of the nanoantenna, can be calculated.
In our method of moments implementation, Green's function for a planar stratified background is considered. Instead of using the Dyadic form of Green's function as in the above equations, a hybrid mixed-potential form is employed [18][19][20][21][22]. In the spirit of the Fourier transform, Green's function is analytically constructed in the Fourier space (the spectral domain) and then numerically transformed back to the real space (the spatial domain). In the latter step, all the singular behaviors in the Fourier space, for example, the branch point, the pole, the slowly decreasing asymptotic behavior when an observation point is close to a source point, are manually removed and analytically converted back to the real space so that the numerical integration is considerably facilitated.
Besides, if a hexahedral block (see Figure 2(c)) is used, a basis function that varies separately in the transverse direction and in the vertical direction is chosen. A transverse basis function reads, In Eq. (14), the functions Δ i x, y ð Þ and Π i z ð Þ are triangular and rectangular functions, In Eq. (15), S i denotes the base rectangle of the i th hexahedral block. r i gives the distance of a point away from a given edge. l i emphasizes the distance between the chosen edge and its opposite edge. In Eq. (16), the superscripts "+" and "À" delimit the upper limit and the lower limit of the block, respectively. A vertical basis function is In Eq. (17), the triangular and rectangular functions are similarly defined as in Eqs. (15) and (16). Notice that Green's function for a general planar multilayer structure behaves differently in the transverse direction (which is described mathematically by a Bessel function and physically by a cylindrical wave) and in the vertical direction (which is described mathematically by an exponential function and physically by a plane wave). For example, a generic form of Green's function is In Eq. (18), k ρ is the transverse wave number and k iz ¼ is the corresponding vertical wave number in the ith layer. A i k ρ À Á represents the amplitude of a field component in the ith layer.
Based on the above facts, the sevenfold integral (including three from the inner volume, and three from the outer volume, and one from the inverse Fourier transform) in the second term of Eq. (12) can be reduced to a fivefold integral. Take the coupling between a horizontal test function and a horizontal basis function as example, ð In Eq. (19), the twofold integral with respect to the vertical direction can be performed analytically. In this way, an enormous amount of numerical efforts can be saved. Please notice that in order to cope with the recently increasing interest in nanoscale antenna arrays, the calculation of periodic Green's functions has been implemented on the basis of a quasi three-dimensional (3-D) mixed potential integral equation formulation [23,24].
The above MoM implementation is extremely efficient for the so-called extruded structures and can be immediately confirmed by experimental data. Take the Nickel (Ni)-based G shape structure (see Figure 3) as an example. In the simulation, the structure is illuminated by horizontally polarized light. As a result, the induced current is highly localized along the structure's diagonal (see Figure 4(a)). Remembering that the current is directly related with the electric field via Eq. (2), a corresponding intensive electric field concentration can be anticipated (see the second harmonic generation image of the structure's near field in Figure 4(b)). This electromagnetic response induces a local thermal effect, that is, heat accumulation (due to the large Ohmic loss in Ni). The generated heat can go beyond the melting point and finally decorates the optical response of the structure. Similar simulation-experiment comparisoncan be found in [26][27][28][29][30][31].

Eigenmodes of a nanoantenna
The impedance operator in Eq. (5) and accordingly the discretized version in Eq. (11) actually include all the electromagnetic properties of the light-nanoantenna interaction problem. As mentioned in the previous section, for a given incident electric field (e.g., a plane wave or the electromagnetic fields emitted by a fluorescent molecule), the induced polarization current in a nanoantenna can be calculated and the associated scattered fields can be readily derived. This analysis scheme indeed works for every nanoantenna.
Nevertheless, it should be noted that the solution is incident field dependent. In some cases, due to the complexity of a nanoantenna's geometry (e.g., the G-shape structure), an excessive amount of information can be read from a single simulation, for example, a very complex near-field pattern, multiple resonances in the extinction, and scattering spectra. In this sense, the incident field-dependent full solution may not be a good starting point for an analysis. To cope with this problem, in this section we pursue two incident field-independent solutions, that is, based on the eigenmodes and natural modes of a nanoantenna, and reveal their distinct roles in determining the antenna's response.
First, we investigate at the eigenvalue problem derived from Eq. (5) [32,33], In Eq. (20), Z is the impedance operator defined in Eq. (5). For a specific angular frequency ω, J n r, ω ð Þis the nth eigenmode with an eigenvalue λ n ω ð Þ. Some extra attention should be paid to the terminology. Here, we define the eigenmodes of a nanoantenna based on the impedance operator Z. This definition is different from the term "eigenmode" (or "normal mode") in quantum optics. There, eigenmodes are defined based on a differential operator (e.g., the differential operator in the Helmholtz equation), while here eigenmodes are defined linked to the integral operator as in Eq. (20).
Under the assumption that the material constituting a nanoantenna is locally responding, homogeneous and isotropic, the total field operator Z tot (defined in Eq. (7)) is a scalar. On this basis, we can distinguish two different contributions to the eigenvalue, Apparently, λ n, scat ω ð Þ is associated with the scattering field operator (defined in Eq. (8)) and hence is controlled by the geometry (as the volume integral in Eq. (8) is taken with respect to the volume of the nanoantenna) and the environment of the nanoantenna (as the involvement of the Dyadic Green's function in Eq. (8)). After removing the total field operator from both sides of Eq. (20), it is obtained that From Eq. (22), we can readily conclude that When the material constituting a nanoantenna is local, homogeneous, and isotropic, the antenna's eigenmodes are indeed independent from the material of the nanoantenna. They only depend on the antenna's geometry and environment.
Due to the reciprocity of Green's function, Z scat is a complex symmetric operator (i.e., the operator contains complex elements and its transpose is equal to itself Z scat ¼ Z T scat ), which further implies that the eigenmodes are complex functions and orthogonal in a pseudo inner product sense [32], and projecting Eq. (5) onto the eigenmodes in a pseudo inner product sense, another matrix equation can be received, In Eq. (25), e m f g is short for the projection of the incident field onto the mth eigenmode, Eq. (25) is very similar to Eq. (11) in the sense that Eq. (11) uses local basis functions to approximate the induced current, while in Eq. (25) eigenmodes which can be deemed as global basis functions are employed. However, different from Eq. (11), which is in general a full matrix, due to the orthogonality property presented in Eq. (23), the impedance matrix z mn f g is a diagonal one. The weighting coefficients in Eq. (25) can be then evaluated as Thus, combining Eq. (24) with Eq. (27), it can be concluded that The induced current is a weighted sum of eigenmodes. The contribution of an eigenmode is determined by two factors: on the one hand the fact whether the incident field can efficiently excite the eigenmode (the denominator of Eq. (27)); on the other hand the absolute value of the eigenvalue, that is, when the absolute value of the eigenvalue reaches its minimum, the corresponding weight reaches its maximum and hence the eigenmode reaches its resonance.
Interestingly, it should be noted that we may use the concept of inner product as well in the construction of Eq. (25). In contrast to the definition of a pseudo inner product, which is the inner product is defined as Comparing Eq. (28) with Eq. (29) immediately reveals the key difference between the notion of the pseudo inner product and the notion of the inner product, namely the extra complex conjugate operation. In classical electromagnetism, the pseudo inner product is more associated with the reciprocal properties, while the inner product is emphasized in the definition of energy/power-related physical quantities, for example, the work done to a scatterer by an incident field.
In this spirit, we can reformulate Eqs.
In Eq. (30), the impedance matrix and the vector containing the projection of the incident field on the mth eigenmode are primed to emphasize that the inner product is applied here. Notice that now the impedance matrix z 0 mn n o is not a diagonal matrix anymore. This is because the eigenmodes are only orthogonal in a pseudo inner product sense but are not orthogonal in an inner product sense.
Further, we can see the right-hand side of Eq. (20) as an "eigen" incident field, Using Eq. (33), Eq. (31) becomes Hence, physically z 0 mn gives the work done by the nth"eigen" incident field on the mth eigenmode and quantifies the energetic coupling between different eigenmodes. We can give the following two definitions, Definition 1: When n ¼ m, z 0 nn ω ð Þ is defined as the self-impedance of the nth eigenmode.
Definition 2: When n 6 ¼ m, z 0 mn ω ð Þ is defined as the mutual impedance between the mth eigenmode and the n th eigenmode.
Quantitatively, by evaluating Eq. (34), a self-impedance and a mutual impedance are related with the eigenvalues, In Eq. (35), d mn ω ð Þ is the inner product of the mth and the nth eigenmodes and represents a unitless factor taking into account the fact that the eigenmodes are not orthogonal in the inner product sense, that is, not energetically orthogonal Note that in Eqs. (30), (31), and (34), the factor 1 2 has been systematically omitted.

Natural modes of a nanoantenna
As a special case of Eq. (5) where the eigenvalue is at zero, we can define the natural mode for a nanoantenna [34,35] In Eq. (37), J α r, ω α ð Þ is the nth natural mode with a natural frequency ω α . Different from the eigenmodes discussed in the previous section, the natural modes are independent of frequency. Mathematically, the natural modes originate from the impedance operator's null space and mark the non-trivial solutions to Eq. (37). Physically, the natural modes (including the induced current and the associated field distribution) correspond to the resonant modes of an open resonator (in contrast to a closed resonator like a perfect electric conductor (PEC) cavity). Because of the radiation damping and the material loss in a nanoantenna, a natural frequency is a complex number (remember in a PEC cavity with no losses the natural frequency is real).
In our MoM algorithm, we can numerically solve for a natural mode by seeking for the (complex) frequency where the determinant of the impedance matrix is zero, We are especially interested in how a natural mode and a natural frequency affect the optical response to an external perturbation at an angular frequency ω. The induced current is obtained by inverting the impedance operator, As in Eqs. (11), (25), and (30), the impedance operator can be recast for a specific set of basis functions. Eq. (39) can be further written as adj Z ð Þ and det Z ð Þ are, respectively, the adjugate and the determinant of the impedance matrix (operator). In general, the determinant is a polynomial that can be decomposed into a product of ω À ω α ð Þ i where i indicates the degree of degeneracy. In the discussion, it is assumed that the natural mode is not degenerate so that i ¼ 1. The induced current can be rewritten as R α is the system residue matrix at a natural frequency. It can be shown that the residue matrix is a dyadic and can be decomposed as the outer product of a natural mode and its transpose, Substituting Eq. (42) into Eq. (41), the induced current can be expressed as Based on Eq. (43), we can clearly define the weighting coefficient for a natural mode, As the natural frequency is a complex number, For an incident field with an angular frequency ω around ω αr , we can approximate the coupling coefficient by In the above approximation, the imaginary part of the natural frequency is assumed to be small. It is readily seen from Eq. (46) that when ω is around ω αr , the weighting coefficient is maximized.
Moreover, consider an incoming electromagnetic pulse with a delta time dependence δ t ð Þ and apply an Inverse Fourier Transform to Eq. (43) with respect to the frequency. In time domain, Eq. (44) reads, In Eq. (47), k α ω ð Þ describes the coupling between a natural mode and an incident field. Due to the complex nature of the natural frequency, the natural mode oscillates at a frequency ω αr but decays with a rate of ω αi . Since the decaying rate is linked with the dissipated power and the stored energy in a scatterer, Accordingly, the "life time" of a natural mode can be defined as 1 À2ωαi . In this way, a modal quality factor Q mod can be defined for each natural mode, Again consider the frequency domain and apply an incident field with an angular frequency ω around ω αr . Assuming that the αth natural mode is dominant, the induced current is That is, the response near the resonance is largely determined by the excited natural mode. Since a fixed stored energy and dissipated power are associated with the natural mode, the modal quality factor Q mod renders a good estimation for the quality factor Q res at resonance.

Examples
In this section, we take real nanostructures to illustrate the above theoretical discussions on the eigenmodes and natural modes of a nanoantenna.

Eigenmodes of a nanobar
Assume a normal incident field polarized along the longest axis of the nanobar structure ( Figure 5). The real part and the imaginary part of the self-impedance for the first three eigenmodes are demonstrated in Figure 6(a) and (b). Note that the coupling coefficient is inversely proportional to the eigenvalue as in Eq. (27) and the eigenvalue is closely related with the self-impedance as in Eq. (35). When the absolute value of the self-impedance reaches its minimum, the coupling coefficient may maximize. Comparing Figure 6(a) and (b) with Figure 6(c), this is indeed the case for the first mode (the L1 mode) and the third mode (the L2 mode). For the L2 mode, the coupling (the integral in the numerator of Eq. (27)) between the incident light and the mode pattern (see the top surface charge distribution for the first three modes in Figure 6(c)) should be taken into account. Since the L2 mode pattern is symmetric with respect to the short axis of the bar but the normal incident field is anti-symmetric, the coupling (the integral in the numerator of Eq. (27)) is zero, that is, the second mode does not contribute to the optical response at all. Further discussions on the symmetry can be found in Section 4. Lastly, comparing the numerically calculated (by the V-MoM algorithm) and experimentally obtained extinction cross sections with the coupling coefficients, one finds that although the coupling coefficients and the extinction cross sections follow the same trend, there is a shift. This shift is actually from the well-known near-field-far-field shift, which is a direct consequence of the radiation loss. Figure 6. The left column illustrates the real part (a) and the imaginary part (b) of the self-impedance for the L1 -L3 mode, coupling coefficients (c) for the L1 mode and the L3 mode and extinction cross sections (d and e), while the right column demonstrates the real part (f) and the imaginary part (g) of the mutual impedance for the L1 mode and the L3 mode, the selfcoupling power (h) of the L1 mode and the L3 mode, the mutual coupling power (i) between the L1 and L3 modes and the extinction cross sections (j and k). In the inset of (c), the top surface charge distribution for the first three modes is plotted. The colour there is coded from blue to yellow to mark the polarity of the surface charge (in print: grey scale). In (d), the V-MoM simulation is marked by the solid line. In (j), the total dissipated power is marked by the solid line.
In the right column of Figure 6, the role of the mutual impedance in determining the antenna's optical response is studied. As discussed in Eq. (34), the mutual impedance describes the energetic coupling between different modes.
First, since the first mode and the third mode (i.e., the odd-order modes which are antisymmetric with respect to the bar's short axis) have a different symmetry compared to the second mode (the even-order modes are symmetric with respect to the bar's short axis), the mutual coupling between the odd-order modes and the even-order modes is zero. Further, because the second mode is not excited by the incident field at all, the second mode does not contribute in any way to the antenna's response.
Second, the first mode and the third mode are energetically coupled (see the mutual impedance in Figure 6(f) and (g)). As a result, the excitation of one of these two modes would directly lead to the excitation of another. Taking this power transfer in between different modes into account breaks the symmetric line shape in the self-coupling power and generates the so-called asymmetric Fano line shape (see Figure 6(h)-(k)).

Natural modes of nanopatches and nanodimer
In this section, three nanostructures, that is, a gold (Au) nanopatch, a nickel (Ni) nanopatch, and a gold (Au) nanodimer, are considered. The nanopatches are fabricated on top of an SiO 2 (100 nm)/Si layer, while the nanodimer sits on an SiO 2 substrate (see Figure 7). In the numerical simulations, the same mesh (5 Â 4 Â 1) is used for the patch structures, while each monomer in the dimer is discretized by a 3 Â 2 Â 1 mesh. For the sake of simplicity, the multilayer substrate is not explicitly taken into account but is modeled by a homogeneous surrounding environment with an effective refractive index (n ¼ 1.25).
Essentially, two observations can be immediately made from the simulations and experiments shown in Figure 8: First, natural modes are fundamental to the optical response of a nanoantenna. Depending on the incident field, some of the modes are excited, while some of the modes are forbidden. For the patch, due to the comparable dimensions along the long axis and the short axis, it is readily expected that the first two resonant modes are horizontal and vertical dipolar modes. Depending on the excitation, that is, horizontally (vertically) polarized normal incident light, the horizontal (vertical) dipolar mode is excited (see the left column in Figure 8). The same argument applies to the dimer structure (see Figure 8). Especially, the second mode (the L2 mode) can be only excited by an oblique incident field.
Second, natural frequencies control the resonances of a nanoantenna. From Figure 8, it is clear that the real part of a natural frequency gives a "baseline" for the resonance, while the imaginary part finally shifts the resonance to the "right" position. Another role that the imaginary part plays is that it is associated with the loss (both radiation and material dissipation) in the structure. A large imaginary part can "flatten" the resonance (e.g., the gold patch has less material loss than the nickel one, so it has a sharper resonance, or the L2 mode in the dimer has less radiation loss than the other two modes, resulting in a sharper resonance).

Symmetries
In this section, we introduce a set of symmetry arguments to determine under which circumstances the eigenmodes are orthogonal to each other in an inner product sense (as defined in Eq. (29)). Since a natural mode is just a special case of an eigenmode, the developed theory applies to the natural mode case as well. In the following, we shall concentrate solely on the symmetries and the eigenmodes of a nanoantenna.
To make our discussion more pedagogical, we again take the simple bar shape structure (see Figure 9(a)) as an example. The bar structure carries a C 2 symmetry group. There are two symmetry operations in this group: the identity operation E where no transformation is conducted, and a rotation of π about the z-axis, the C 2 symmetry operation.
Since symmetry operations are always applied based on coordinates, we should be able to find a corresponding set of matrices to "represent" these operations. Here, we especially focus on the matrices with the lowest dimensionalities, that is, the irreducible representations. Since the group under discussion is an abelian group, we have two irreducible representations and they are shown in Figure 9(b).
Moreover, in contrast to these transformations operating on coordinates, we follow Wigner's conventions [36,37] and define transformation operators which operate on functions, In Eq. (51), this definition is illustrated for both scalar functions (such as charge, etc.) and vector functions (such as currents, electromagnetic fields, etc.). These transformation operators are commutative with the impedance operator defined in Eq. (5), that is, the sequence of applying the transformation operator and the impedance operator does not affect the final outcome of the calculation, Figure 9.
Combining the group's irreducible representations and its transformation operators, we can further construct a set of projection operators for the group under discussion, In Eq. (53), a projection operator is characterized by the subscript j which marks an irreducible representation. Here, j may run from one to two. l j is the dimensionality of an irreducible representation. Since every irreducible representation has a dimensionality of one, l j is equal to one. Then, the summation is carried out with respect to all the symmetry operations.
Based on the defined projection operator, we can categorize all the eigenmodes according to the irreducible representations (see an illustration of an eigenmode categorization in Figure 9(c)). It is well known that [36,37] the functions (e.g., in the current case, the eigenmodes) that belong to different irreducible representations are orthogonal to each other in an inner product sense. Figure 10. The mutual impedance (a) -(p) between the first four modes (see Figure 9) of the bar structure. This figure is organized as follows. A column represents a "source" mode that generates eigen-incident field as in Eq. (33) while a row represent an "observation" mode that the eigen-incident field is applied onto. Thus, the figure at the cross of the mth row and the n th column denotes the mutual impedance z' mn as defined in Eq. (34). This figure is reproduced from [36]. This is the very condition that controls the modes' orthogonality (in an inner product sense). The above is equivalent to saying that, according to the irreducible representations, an eigenspace (in which one finds all the eigenmodes) can be split into several invariant subspaces (in which one finds all the eigenmodes belonging to an irreducible representation).
To demonstrate this condition, the energetic coupling among the first four eigenmodes for the bar structure is calculated with the V-MoM algorithm discussed in this chapter. Upon inspecting Figure 10, it can be immediately concluded that the even-order modes do not energetically couple with the odd-order modes. This is the same as we have seen in Section 3.3.1, where the vanishing energetic coupling is attributed to the integration of the product of an odd function and an even function. Now, this intuitive idea is formally incorporated in the framework of group theory. That is, since the even and the odd eigenmodes belong to different irreducible representations, they must be orthogonal in an inner product sense.
Following the same line of reasoning as in Section 3.3.1, an excitation is applied (see Figure 11). The incident field is polarized along the bar's long axis but the incident direction is oblique. As already seen in Figure 6, power transfer between modes of the same parity, that is, the even Figure 11. The power transfer (a) -(p) between the first four modes (see Figure 9) of the bar structure. This figure is organized similar to Figure 10. A column represents a "source" mode that generates eigen-incident field as in Eq. (33), while a row represent an "observation" mode that the eigen-incident field is applied onto. However, different from Figure 10, the coupling coefficients of the "source" and the "observation" modes are included. Thus, the figure at the cross of the m th row and the n th column denotes the power transfer p' mn as defined in Eq. (54). This figure is reproduced from [36].
(odd)-order modes, is allowed and can be evaluated by including the effects of the different eigenmodes' coupling coefficients, p 0 mn ¼ c Ã m ω ð Þc n ω ð Þd mn ω ð Þλ n ω ð Þ: ð54Þ Eq. (54) is a direct extension of Eq. (35). However, as stipulated by group theory, the power transfer between modes of different parity is forbidden. Lastly, as a comment, different from Figure 6, due to the obliquely incident field, the even-order modes are now excited, and hence the resonance and the energetic coupling between the second mode and the fourth mode is seen.

Conclusions
This chapter presents a complete set of numerical tools to tackle the nanoscale light-matter interaction problem. By comparing with experimental results, the V-MoM solver and its extensions are proven to be very efficient and accurate. The modal analysis scheme and the group theoretical-based analysis scheme are demonstrated to deliver a deep physical understanding of the physical system under study. Therefore, we believe that the developed integral equationbased solver may become an essential complement for commercial solvers which are mostly based on differential equations.