## 1. 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–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 pre‐calculated 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.

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

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

Further, it should be realized that the scattered field is induced by the induced current (as in Eq. (2)),

In Eq. (3),

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

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

Approximating the induced current

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

In Eq. (11), *j*th basis function.

Notice that in Eqs. (12) and (13), the integrations are with respect to the volume of the *i*th *basis function*, *j*th *test function,*

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

In Eq. (15), *S*_{i} denotes the base rectangle of the *i*^{th} hexahedral block. ϱ_{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), *i*th layer.*i*th 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,

(19) |

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 comparison—can be found in [26–31].

## 3. Eigenmodes and natural modes of a nanoantenna

### 3.1. 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 *n*th eigenmode with an eigenvalue *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

Apparently, *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, *complex symmetric* operator (i.e., the operator contains complex elements and its transpose is equal to itself

By expanding the induced current in terms of the eigenmodes,

and projecting Eq. (5) onto the eigenmodes in a pseudo inner product sense, another matrix equation can be received,

In Eq. (25), *m*th 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

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. (24)–(27) while using the concept of inner product,

In Eq. (30), the impedance matrix and the vector containing the projection of the incident field on the *m*th eigenmode are primed to emphasize that the inner product is applied here. Notice that now the impedance matrix **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 *n*th“eigen” incident field on the *m*th eigenmode and quantifies the energetic coupling between different eigenmodes. We can give the following two definitions,

*Definition 1: When *

*Definition 2: When *

Quantitatively, by evaluating Eq. (34), a self‐impedance and a mutual impedance are related with the eigenvalues,

In Eq. (35), *m*th and the *n*th 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

### 3.2. 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), *n*th natural mode with a natural frequency

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

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

*i* indicates the degree of degeneracy. In the discussion, it is assumed that the natural mode is not degenerate so that

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

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

Moreover, consider an incoming electromagnetic pulse with a delta time dependence

In Eq. (47),

Accordingly, the “life time” of a natural mode can be defined as *Q*_{mod} can be defined for each natural mode,

Again consider the frequency domain and apply an incident field with an angular frequency

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.

### 3.3. Examples

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

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

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 anti‐symmetric 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)**).

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

## 4. 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 *z*‐axis, the

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,

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.

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.

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 (odd)‐order modes, is allowed and can be evaluated by including the effects of the different eigenmodes’ coupling coefficients,

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.

## 5. 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 equation‐based solver may become an essential complement for commercial solvers which are mostly based on differential equations.