1. Introduction

Electromagnetic resonators are omnipresent in science and engineering and come in diverse sizes and shapes ranging from microwave resonators, via cavities for gas and semiconductor lasers, to optical microcavities and plasmonic nano-resonators. Common to all of them is the fact that the resonances—i.e., the discrete set of special frequencies that may show up as peaks in scattering spectra—all have associated electromagnetic field distributions, which are often referred to as “modes” of the resonators. In all realistic physical resonators, moreover, there is a certain degree of dissipation of energy to the environment. This effect is responsible for a broadening of the peaks in the spectra and is typically quantified in terms of the so-called quality factor $ Q $. The higher the $ Q $ factor, the longer it takes before an initial excitation of the resonator has dissipated away—either through radiation to the environment or absorption in the material. From the panoply of various physical systems exhibiting electromagnetic resonances, it is not obvious that a common mathematical framework exists, which can be used to precisely capture their physical properties in terms of the resonant and possibly dissipative modes. Nevertheless, such a framework does exist, and it is the goal of this Tutorial to present it in some detail along with a number of relevant applications. Indeed, from a modeling perspective, it is an interesting fact that the dissipative modes of electromagnetic resonators can be calculated as solutions to a specific eigenvalue problem, namely the sourceless Maxwell wave equation subject to a radiation condition that allows only solutions propagating away from the resonator at large distances. These solutions, which have complex resonance frequencies, are known in the literature as resonant states [1,2], morphology-dependent resonances [3,4], or quasinormal modes (QNMs) [5,6].

1.1. Brief Comment on Nomenclature

In the history of physics modeling, the first use of solutions to the wave equation having complex energies is commonly attributed to Gamow [7], and Zel’dovich appears to be the first to have introduced a method for the normalization of such solutions [8]. Both Gamow and Zel’dovich were concerned with problems occurring in quantum mechanics, where the solutions are typically referred to as “resonant states” or “resonance states,” yet the wave nature of the Schrödinger equation implies that mathematical approaches developed in this context are useful also in electromagnetism [9]. “Morphology-dependent resonances” is the name originally given to the resonances observed in various optical microparticles, and which clearly depend on the shape and structure of the particles. For spheres, in particular, the resonances can be immediately related to poles in the Mie expansion coefficients [10] and the associated dissipative modes of the spheres [11]. By comparison to the Schrödinger equation, one can appreciate that the electromagnetic problem of the sphere is very similar to the problem of an electron in a spherical potential well, only with a radially dependent potential, which vanishes at large distances [4]. The terminology “quasinormal modes” appears to have been first used in the context of decaying perturbative solutions to the Einstein equations close to a black hole [12,13]. The similarity between the linearized Einstein equations and Maxwell’s equations [14], as well as the open nature of the universe around a black hole, means that the mathematics are similar to that of electromagnetic resonators in free space. Indeed, some results originally derived for optical resonators have been subsequently applied also to gravitational systems [15]. We shall use the term “quasinormal mode (QNM)” to emphasize the fact that in many ways the modes of interest represent a precise generalization of the well-known modes of closed cavities (with infinite $ Q $ factor) to the case of general dissipative resonators, in which the local electromagnetic energy may be lost in the form of radiation to the environment and possibly absorption in the material.

1.2. Motivation for QNM Models

In many expositions of resonator optics, the dissipation of the cavity mode is treated in a perturbative manner, where the cavity mode is first calculated for the closed cavity with no absorption and subsequently coupled to the environment via mostly phenomenological coupling constants. Such models, for example, have been extremely successful in modeling optical cavities with high $ Q $ factors. Also, when the evolution of electromagnetic resonators has led to smaller and smaller devices (in particular in the broad research area known as nanophotonics), the associated modeling has largely been based on ideas rooted in the perturbative approach. The strive for smaller resonators (often quantified by the effective mode volume $ {V_{\text{eff}}} $) has largely been driven by the associated increase in the attainable electromagnetic field strength, even at the expense of lower $ Q $ factors, as long as the ratio $ Q/{V_{\text{eff}}} $ remains large; a lower $ Q $ factor can even be beneficial, as it relaxes the fabrication tolerances or allows larger bandwidth operation. The limit of low $ Q $ factors, however, is exactly where the perturbative closed-cavity resonator models break down, and where a model of the cavity modes in terms of QNMs may be particularly advantageous. In the case of metallic nano particles, for example, the QNMs describe the localized surface plasmon polaritons that are supported by the nano structures.

A QNM framework is attractive from a conceptual as well as a computational point of view, since in the general case there may be no practical way of calculating or even defining a closed cavity. This is the case, for example, in the second introductory example in Subsection 1.5b. Even in cases where one can make such a definition, the coupling to the environment may lead to substantial frequency shifts and field distortions, which cannot be accurately predicted based on closed-cavity modes. Moreover, typical descriptions based on QNMs are no more complicated than the perturbative closed-cavity models—the main difference is that a QNM framework usually provides an explicit and precise way of calculating the various coupling parameters entering the model. In cases where a single QNM dominates the response, for example, the result of a QNM model of the Purcell factor can be immediately written in exactly the same way as the original formula due to Purcell [16], but with a slightly modified expression for the effective mode volume [17]; this definition of a mode volume for a leaky resonator can even be extended to plasmonic systems by properly accounting for absorption and dispersion in the material [18]. Finally, in systems described by several QNMs, there will generally be a phase difference between the complex QNM amplitudes, and so a formal framework based on QNMs may be very useful.

1.3. Overview of the Existing Literature

Below we attempt to provide an overview of the existing literature on QNMs. As already mentioned in the introduction, the generality of the wave equation means that the usefulness of a framework based on dissipative modes extends through many branches of physics. In keeping with the scope of this Tutorial, however, we shall limit the overview to electromagnetism.

1.3a. Theoretical Developments

Due to the relatively large computational costs of calculating QNMs for general structures, an impressive body of work has been developed for problems in one dimension as well as the analytically tractable cases of cylinders and spheres in two and three dimensions, where the QNMs can be calculated with relative ease. Perturbation theory for QNMs was presented by Lai et al. in Ref. [19], and, in a series of papers, Leung et al. have treated completeness [20–22], perturbation theory [23,24], and dispersive materials [25]. Lee et al. discussed the QNM completeness and electromagnetic Green tensor expansions [26] as well as perturbation theory [27] with applications to dielectric spheres; the case of degenerate perturbation theory was later discussed in Ref. [28]. Recently, the interesting problem of perturbative changes in the surrounding material was addressed by Both et al. [29], and a variant of perturbation theory with deformed resonators was presented in Ref. [30]. In Ref. [31], both one and three-dimensional problems were treated by Muljarov et al. using direct expansions in a subspace of QNMs—a method generally referred to as resonant state expansion [32]—to calculate the effects of material changes; this method was later extended to the case of dispersive media in Ref. [33] and was used for the study of planar photonic crystal structures in Ref. [34], see also Ref. [35] for a recent contribution. The one-dimensional problem was also investigated by Settimi et al. in Refs. [36,37], and Doost et al. have treated slabs [38], cylinders [39], and spheres [40] by a QNM expansion of the Green tensor. For modeling using a biorthogonal basis in one dimension, see also Refs. [41,42]. Recently, the completeness of QNM expansions in spheres made from dispersive materials was investigated by Mansuripur et al. [43]. Armitage et al. used a QNM expansion of the Green tensor to model planar waveguides with oblique incidence of light in Ref. [44], and Lobanov et al. recently suggested the use of a resonant state expansion based on the QNMs of a sphere in combination with the Dyson equation to calculate the scattering properties of general resonators [45].

For treating QNMs in general structures, a variety of numerical methods have been employed and are still under active development for both QNM calculation and normalization. These include the use of volume [17,46] or surface [47–49] integral equation formulations, as well as the Fourier Modal Method (FMM)—also known as rigorous coupled wave analysis—for periodic structures [50–52], or for single resonators by use of so-called perfectly matched layers (PMLs) [18,53]. In two-dimensional coupled cavity-waveguide structures, the QNMs have been calculated by a Dirichlet-to-Neumann technique in Ref. [54], by FMM [55], or by the finite element (FEM) calculations with a nonlocal boundary condition in Ref. [56]. The latter work also discussed normalization of the QNMs via the theory of divergent series. Römer et al. [57] used QNMs obtained as the solutions to the wave equation in FEM calculation with PMLs to study spontaneous emission from emitters in photonic crystal cavities. To this end, the QNMs of interest were normalized by a volume integration, which extended through the PMLs, a technique that was later also used by Sauvan et al. in combination with FMM calculations [18,53]. A FEM formulation leads essentially to one numerical eigenmode per degree of freedom in the problem. The vast majority of these eigenmodes, however, are connected with the PMLs and at first sight do not appear to be useful for modeling. Nevertheless, as shown by Vial et al. for two-dimensional open systems [58], the full set of modes containing the dominant QNMs and the auxiliary PML modes provide a useful basis for expansion of the solutions to the wave equation. These ideas were subsequently further developed by Yan et al. [59] in three dimensions. Muljarov et al. have contributed to the discussion of usefulness of various normalization methods in Ref. [60] (see also Refs. [61–63]). The problem of normalization was recently revisited by Stout et al. [64] by means of another possible regularization of the normalization integral originally due to Zel’dovich [8].

Numerical eigenmode calculations typically require significant computational resources. As an alternative, therefore, direct calculations of the electromagnetic response at complex frequencies may be advantageous, as was pointed out by Bai et al. [65] and Perrin [66]. Such approaches are similar in spirit to the so-called Riesz projection approach [67–69], which, in turn, is closely related to the Green tensor expansion usually employed in the resonant state expansion literature. Although the QNMs are defined in the frequency domain, they have been successfully calculated also with time domain methods—notably the finite difference time domain (FDTD) method—in combination with PMLs and a Fourier transformation [17,70,71].

An inherent exponential divergence of the QNMs at large distances from the resonator means that they cannot be directly used to describe electromagnetic fields in this limit. Instead, the use of a Dyson equation approach based on QNMs has been suggested [70,72]. The far-field problem has also been treated using causality requirements by Abdelrahman et al. in Ref. [73] and Colom et al. in Ref. [74], and in the context of coupled mode theory (CMT) in Ref. [75]. The problem is closely connected to that of QNM-based scattering calculations, which have been treated in a number of ways in Refs. [45,59,65,66,70,76–80]. Similar to the divergence at large distances, the dramatic increase of the electromagnetic feedback close to metal surfaces cannot be captured by a single QNM. In such cases, the response can be conveniently handled by an additional quasi-static contribution to the Green tensor [70], see also Ref. [59].

1.3b. Practical Applications

Even with some fundamental questions still unsolved, there has been a recent bloom in the use of QNMs for practical modeling tasks in nanophotonics. Many of these applications can be seen as refinements of previous calculation methods for resonant structures, where now the QNMs provide explicit and precise definitions of parameters that would normally be inferred by fitting to simulation data or measurements. The close connection between the QNMs and the resonances in scattering matrices of general structures has been clarified [45,76–79], and QNMs have been used for an alternative derivation of the so-called temporal CMT [75], including applications to switching using nonlinear materials [81], see also Ref. [82] for a recent study of nonlinear effects using QNMs. Yang et al. used QNM perturbation theory for sensing applications [83], and QNMs of resonators modeled with a nonlocal material response were presented in Refs. [69,84]. The QNMs in coupled cavity-waveguide systems were used for perturbation theory and Purcell factor calculations [56] as well as broadband local density-of-states (LDOS) calculations in two-dimensional systems [85]. Malhotra et al. subsequently applied the theory in three dimensions to describe on-chip single photon emitters [86]. Also, QNM descriptions have been employed for LDOS calculations in hybrid plasmonic photonic structures [87] and used for theoretically predicting [88] and interpreting [89] electron energy loss spectroscopy results.

Classical laser cavity models can be said to implicitly rely on a QNM picture, namely that of a cavity with a multitude of resonant modes, each of which repeats itself after a full round trip. Partial transmittance at the cavity end facets means that part of the electromagnetic energy escapes, leading to complex resonance frequencies. By compensating the radiative energy loss through a gain medium, the resonance frequency can be shifted towards the real axis, which leads to the characteristic linewidth narrowing at the onset of lasing [90,91], and QNMs were used for calculating the parameters in a Maxwell–Bloch model for a photonic crystal cavity laser in Ref. [92]. In the limit of single energy quanta, microscopic semiclassical theories of light–matter interaction can be set up in a number of ways, often based on the Green tensor. In these limits, a QNM expansion immediately leads to relatively simple and physically appealing models for modified spontaneous emission calculations, in particular the Purcell factor [17,18,60,65–67,71,85,93–95]. As an alternative to methods based on the Green tensor, the QNMs have been used as input to single mode master equations [49,96–98]. The problem of full quantum models based on QNMs has been discussed by Ho et al. [99], Dutra et al. [100], Severini et al. [101], and recently Franke et al. [102].

To round off this overview, we remark that the QNMs of localized electromagnetic resonators are related to the modes of optical fibers and general waveguides, where absorption or radiation will also lead to inherently dissipative modes [9,103–106], and one can conveniently perform projections by use of suitably defined adjoint modes as we do in the derivations to follow. In addition, as mentioned above, QNMs have been studied in areas of physics other than electromagnetic resonators, notably in quantum mechanics [107] and in general relativity [108], but also in acoustics [109]. For a number of reviews on QNMs, see Refs. [6,110,111].

1.4. Scope and Structure

In this Tutorial, we present a variant of the well-established expansion and projection method using biorthogonal modes [112] to model light scattering by localized electromagnetic resonators in terms of QNMs. At first sight, this approach, which is independent of the particular calculation technique used to obtain the QNMs, suffers from the fact that the QNMs in two and three dimensions obey the asymptotic requirement of a radiation condition instead of a boundary condition. For this reason, the traditional approach [112] is not directly applicable in dimensions higher than one, even if variants of it have been successfully used in electromagnetics for the study of lossy periodic structures [113,114]. Nevertheless, as we shall see, it is possible to extend the theory to higher dimensions by ideas originally developed in the literature on resonant states. For definiteness, we limit most of the derivations and illustrations to resonators of finite size embedded in homogeneous materials, where the scattered electromagnetic fields obey the Silver–Müller radiation condition. This does not include all technologically relevant geometries. Indeed, the modeling of infinite planar array structures or membranes at oblique incidence requires the use of periodic boundary conditions in the in-plane directions [51,52,76,115,116], and the electromagnetic field in resonators on top of substrates [59,95] or coupled to optical waveguides [54–56,85,86] obey different radiation conditions, which are not immediately covered by the theory laid out in this article. Such cases must be treated either by modifications of the theory using suitable boundary or radiation conditions, or by theoretical approaches specific to the calculation method of the QNMs. References [18,58,59], for example, describe modal techniques in which the calculation domain truncation by use of PMLs is an integral part of the modeling approach, in particular of the normalization procedure. Similarly, the normalization scheme by use of divergent series in Refs. [56,85,86] arises naturally from the waveguide radiation condition of the QNMs in coupled cavity-waveguide systems. These and other interesting cases are beyond the scope of this Tutorial, and we refer instead to the literature, cf. Subsection 1.3.

In Subsection 1.5, we present two introductory examples of practical QNM modeling applications; transmission through a dielectric barrier in one dimension, and Purcell factor calculations for a plasmonic dimer of nano-spheres in three dimensions. In the course of the Tutorial, we shall repeatedly return to these example structures to exemplify the various calculations. The one-dimensional example is sufficiently simple, so that most of the calculations can be handled analytically, and we much encourage the interested reader to repeat them. As a supplement to this Tutorial, we provide a number of MATLAB files implementing some of the one-dimensional examples (Code 1, Ref. [117]) as well as the code necessary for calculating QNMs of three-dimensional resonators (Code 2, Ref. [118]) using the freeware code MNPBEM [119–121].

The remainder of this Tutorial is organized as follows. In Section 2, we discuss various methods currently in use for practical QNM calculations. Section 3 lays out the basic elements of the theory for general, three-dimensional resonators in homogeneous environments. We define the QNMs as the solutions to the sourceless wave equation subject to the Silver–Müller radiation condition and show how this requirement naturally leads to the definition of adjoint QNMs and a special operator for projection of certain solutions to the wave equation onto the QNMs. Using the projection operator, we construct formal expansions of general electromagnetic fields and the electromagnetic Green tensor in terms of QNMs. In addition, we discuss how the QNMs are directly related to the residues of the Green tensor, as typically exploited in the literature on resonant states, and we show how such an approach leads naturally to a convenient alternative normalization procedure for the QNMs. In Section 4, we discuss the question of convergence of the formal expansions in terms of QNMs, which can be assessed by the limiting behavior of the Green tensor. We discuss how an investigation of the Green tensor can be used to define a region of consistency for the QNM expansions, and how one can subsequently exploit this knowledge to extend the region of consistency. Section 5 is devoted to a number of practical applications of QNM modeling. Subsection 5.1 presents the CMT and scattering calculations, Subsection 5.2 is concerned with hybridization and discusses how one can expand the QNMs of coupled resonators in terms of the QNMs of the individual resonators, Subsection 5.3 presents various applications of perturbation theory, and Subsection 5.4 discusses the use of QNMs for Purcell factor calculations. Finally, Section 6 holds the conclusions. In three appendices, we provide additional information. Appendix A gives an introduction to the use of Green functions and Green tensors, Appendix B discusses the theory and practical application of convergence studies, and Appendix C collected certain lengthy but instructive calculations for which we display only the final results in the main text of the Tutorial.

1.5. Introductory Examples

To illustrate the usefulness of the QNMs in modeling responses of electromagnetic resonators, we consider two examples—a one-dimensional example of a dielectric barrier and a three-dimensional example of a dimer made from two gold nano-spheres. We shall repeatedly return to these introductory examples throughout the article to illustrate the application of the various results.

1.5a. Dielectric Barrier in One Dimension

We consider first the classical example of a dielectric barrier with constant refractive index $ {n_{\text{R}}} = \pi $ in a background with refractive index $ {n_{\text{B}}} = 1 $. For definiteness, we consider propagation along the $ x $ axis, and we take the electric field to point in the $ y $ direction; the barrier has a width of $L$ and is centered around the origin. The top panel in Fig. 1 shows the transmission through the barrier as a function of frequency. It features a number of distinct peaks with unity transmission; these are the well-known Fabry–Perot resonances, each of which can be directly associated with a QNM. The center panel shows a part of the complex frequency spectrum below the real axis, in which the QNM frequencies show up as dark spots. For this particular choice of refractive index, the real part of the QNM frequencies are evidently spaced equally by $ L/\text{c} $, and all QNM frequencies have the same imaginary part. The bottom panel of Fig. 1 shows an example of a QNM field in the vicinity of the resonator. As one might expect, the QNMs represent standing waves inside the resonator, as would be the case also for a closed cavity made from perfect reflectors at each side. Contrary to the closed-cavity case, however, the field is not zero outside, where it propagates away from the resonator at both sides. Notably, the field magnitude clearly increases in the direction away from the resonator, which is a general feature of QNMs and necessitates the use of a normalization that is different from that of the closed-cavity case. Once normalized, however, the QNMs can be used to calculate quantities of typical interest in the description of electromagnetic resonators. As an example, the dashed red curve in the top panel shows the approximate transmission as calculated using the three QNMs with $ 3 \text{c}/L \le \text{Re}\{ \tilde \omega \} \le 5 \text{c}/L $. The approximation is best in the center of the frequency interval, where the relative error is less than 1%. By systematically increasing the number of QNMs, the error can be made arbitrarily small, as we shall see in Subsection 5.1c.

 

Figure 1. Top: transmission through the dielectric barrier as a function of frequency for the case of $ {n_{\text{R}}} = \pi $. The dashed red curve shows the approximate transmission calculated using three QNMs. Center: complex resonance frequency spectrum, in which the logarithmic color coding indicates the magnitude of the left-hand side of Eq. (4). In this way the discrete QNM frequencies show up as dark spots, indicating the zeros of the function in the lower half of the complex plane. Bottom: mode profile of the electric field QNM with $ \text{Re}\{ \tilde \omega \} = 4\text{c}/L $. The red solid and blue dashed curves show the real and imaginary parts of the field, respectively, and the black curve shows the absolute value, which increases (exponentially) as a function of distance from the resonator. The gray shading indicates the extent of the dielectric barrier.

Download Full Size | PDF

1.5b. Plasmonic Dimer of Gold Spheres

As a second example, we consider a dimer made from two gold spheres of radius $ R = 50\,\text{nm} $ embedded in air and separated by a distance $ d = 50\,\text{nm} $. The electromagnetic response of gold is modeled with a Drude model for the relative permittivity of the form

 

Figure 2. Top: Purcell factor of a dimer made from gold nano-spheres (see inset) as a function of frequency. Red curves show the approximate responses calculated using a single QNM corresponding to $ {\tilde \omega _1} $ (dashed) or the two QNMs corresponding to $ {\tilde \omega _1} $ and $ - \tilde \omega _1^* $ (solid). Center: complex QNM spectrum, in which the logarithmic color coding indicates the smallest magnitude eigenvalue of the operator $ {{\bf M}_{\bf G}}(\omega ) - 1 $, cf. Eq. (8). In this way, the QNM resonance frequencies show up as dark blue spots indicating the position of the zeros of this function. The lowest-order eigenvalues of interest are indicated with black circles. Bottom: field profiles showing the absolute values of the three QNM fields corresponding to the frequencies indicated in the spectrum above. Only $ \tilde{{\bf f}}_1({\bf r}) $ contributes to the Purcell factor in the gap center.

Download Full Size | PDF

2. QNM Calculation Methods

Although the dissipative nature of physical electromagnetic resonators is easy to accept, the introduction of associated dissipative modes of the electromagnetic field appears to be less intuitive. Indeed, scientists will routinely refer to the $ Q $ factor of a resonator, but rarely to the $ Q $ factor of the dissipative mode of the electromagnetic field associated with the resonator. One conceptual difficulty, it appears, is that the word “mode” is sometimes reserved for solutions to Hermitian eigenvalue problems, and any mode associated with dissipation would fall outside this category. It is conceptually very fruitful, however, to broaden the scope of the word “mode” to mean simply the eigenfunctions of a given differential equation problem. In electromagnetism, for example, this could be the eigenfunctions of the sourceless Maxwell’s equations. A central point in the theory of differential equations, however, is that a differential equation in itself does not constitute a well-defined problem; only by introducing a set of boundary or radiation conditions do we obtain a sufficiently rigorous statement of the mathematical problem, the solutions to which we can then refer to as “modes.” From this point of view, the modes of a closed cavity belong to a different class of modes than the modes of a leaky cavity, because they fulfill different boundary conditions. Mathematically, therefore, the conceptual as well as computational difficulties associated with QNMs can be seen as a consequence of the fact that the QNMs obey a radiation condition instead of a boundary condition, and radiation conditions are comparably difficult to handle, even if important progress has been made over the past decades.

In practical QNM calculations, one typically does not enforce the radiation condition explicitly, but rather imposes the correct radiation behavior by one of several possible strategies. One option is to directly look for solutions to the wave equation in terms of analytic continuation of functions obeying the correct radiation condition at real frequencies. For spheres in free space, in particular, this provides an analytically tractable approach for calculating the QNMs in terms of spherical wave functions [19,31], and we use this approach also for the case of the dielectric barrier in one dimension below. A numerical variant of this approach is QNM calculations via a FMM framework [18,51–53,55]. Another option is to calculate the QNMs via an integral equation formulation [17,46,48,49], where the radiation condition is built in via analytical continuation of the Green tensor, which manifestly fulfills the correct radiation condition. A third option, which is likely the most widely adopted, relies on PML-type truncations of the calculation domain to eliminate as much as possible reflections at the domain boundaries and thereby effectively emulate radiation in free space [17,18,65,116]. This calculation method has traditionally been the workhorse for many optical cavity calculations, but the lack of an explicit correspondence with the radiation condition meant that many properties of the modes were not broadly appreciated—in particular the divergent nature of the QNMs at large distances. The explicit correspondence between QNM calculations with an integral equation formulation and with PML truncations was illustrated in Ref. [17], and Maes et al. [122], de Lasson et al. [123], and Lalanne et al. [124] recently compared a number of distinct calculation methods and implementations commonly in use in the literature. In addition to direct numerical solutions, numerous approximate methods exist, which may be used to devise simple analytical descriptions of the QNM fields as well as insight to the physical mechanisms responsible for the partial trapping of the electromagnetic field. Models of linear defect cavities have been presented in Ref. [125] and plasmonic nano-rods and resonators are treated in Refs. [126,127], respectively. For resonators with high $ Q $ factors, such as microtoroids, for example, one can also obtain approximations to the QNMs by calculating the modes of the resonator embedded in a closed cavity and subsequently estimate the radiative energy loss by integration of the field on the resonator surface [128].

In Subsections 2.1 and 2.2 below, we present details of the calculations for the QNMs of the two example material systems from Subsection 1.5. We note that the dielectric barrier is a classical QNM example, which has been presented in numerous publications, including Refs. [31,111]. The QNMs in a variant of the plasmonic dimer of gold nano-spheres were previously presented in Ref. [46].

2.1. QNMs of the Dielectric Barrier

A general electromagnetic field in a one-dimensional system with piecewise homogeneous materials can be expanded in forward and backward propagating plane waves. For the dielectric barrier of width $ L $ and refractive index $ {n_{\text{R}}} $ embedded in a background of refractive index $ {n_{\text{B}}} $, we find that the $ m^{\prime}$th electric field solution to the wave equation with the requirement of purely outgoing waves in the regions left and right of the resonator can be written in the form

As an alternative to the analytical approach taken here, we could have solved the problem numerically with algebraic boundary conditions connecting the electric and magnetic QNM fields,

This is always possible in one dimension, where the Silver–Müller radiation condition turns into a regular boundary condition, which can be applied at any distance from the resonator.

Supplementary Code. As a supplement to this Tutorial, we provide MATLAB code to conveniently calculate and plot the QNMs of the dielectric barrier [117], and we encourage interested readers to investigate the variation in the QNM resonance frequencies and field distributions as a function of the mode index $ m $ or the permittivity difference.

2.1a. Limit of Vanishing Permittivity Difference

It is instructive to consider the limit of vanishing permittivity difference between the resonator and the surrounding medium, for which one might expect to recover the plane wave solutions of the wave equation in homogeneous media. From Eq. (5), it follows immediately that the spacing of the real parts of the complex resonance frequencies tend to $ \Delta \omega = \pi /{n_{\text{B}}} $, whereas the imaginary parts diverge along paths parallel to the negative imaginary axis. The logarithmic divergence, however, is extremely slow. For a 10% difference in refractive index, setting $ {n_{\text{R}}} = 1.1 $ and $ {n_{\text{B}}} = 1 $, we find $ {\gamma _m}L/\text{c} \approx 2.8 $. Reducing the refractive index difference by a factor 10 by setting $ {n_{\text{R}}} = 1.01 $ and $ {n_{\text{B}}} = 1 $, we find $ {\gamma _m}L/\text{c} \approx 5.3 $.

As the refractive index contrast between the dielectric barrier and the surrounding medium vanishes, the QNMs change in a continuous fashion and become more and more leaky as the imaginary part of the resonance frequency tends to minus infinity (albeit extremely slowly). In this limit, the real-frequency electromagnetic response due to a given QNM blends with that of the other QNMs, in a subtle way mimicking the continuous spectrum of the free one-dimensional problem. The individual QNMs, however, retain the form in Eq. (2), which is manifestly different from the plane waves because the QNMs obey different boundary conditions than the plane waves. Also, whereas the plane waves can be used as a basis along the entire real line, the QNM expansions converge to the correct solution only inside or close to the resonator, as we shall discuss further in Section 4.

2.2. QNMs of the Plasmonic Dimer

Even though the plasmonic dimer consists of only two spheres, the geometry is too complicated for any analytically tractable approach, and we must resort to numerical calculations. All results in Fig. 2 were calculated with the volume integral equation (VIE) method described in Ref. [46]. This method is specialized to the problem of collections of spherical scatterers, for which it provides relatively precise results. To illustrate the generality of the results, and to comment on the issues of convergence and consistency, we also carry out nominally identical calculations with the use of the freely available code MNPBEM by Hohenester [119–121], which is an implementation of a boundary element method (BEM) formulation by de Abajo and Howie [129]. For the QNM calculations using MNPBEM, we follow the approach in Ref. [49], and we provide supplementary code to enable the reader to repeat several of the calculations [118]. Both integral equation methods benefit from the fact that they are defined by use of the electromagnetic Green tensor. Consequently, they manifestly respect the correct radiation condition, and the numerical errors are therefore expected to derive primarily from the spatial discretization.

2.2a. Calculations Using Spherical Wave Functions

In the VIE formulation of Ref. [46], a general electric field QNM is expanded in spherical wave functions inside each sphere $ i $ as

 

Figure 3. Complex QNM spectrum for the dimer of gold nano-spheres showing, at each position $ \omega = {\omega _{\text{R}}} + \text{i}{\omega _{\text{I}}} $ in the lower half of the complex plane, the logarithm of the smallest magnitude eigenvalue of the operator $ {{\bf M}_{\bf G}}(\omega ) - 1 $, cf. Eq. (8). The QNM frequencies show up as distinct dark blue spots indicating the minima in the modal landscape. The black frame in the upper right corner indicates the extent of the frequency map displayed in Fig. 2. Evidently, the entire spectrum extends far below this region.

Download Full Size | PDF

The spectrum is clearly mirror symmetric with respect to the imaginary axis, which is a general feature of the QNM spectrum, and we shall return to it in Subsection 3.1. The possible degeneracy of the QNMs is fully captured by Eq. (8), but is not immediately clear from Fig. 3.

The accuracy is governed by the cutoff parameter $ {l_{{\max}}} $. In practice, therefore, the complex resonance frequency can be thought of as a (discrete) function of $ {l_{{\max}}} $, and to estimate the accuracy of a given calculation, we consider the difference between results obtained using $ {l_{{\max}}} $ and $ {l_{{\max}}} + 1 $,

For the case of the complex QNM resonance frequency $ {\tilde \omega _1} $, Fig. 4 shows the logarithm of $ {D_{ + 1}} $ as a function of $ {l_{\text{max}}} $. To a very good approximation, the points corresponding to large values of $ {l_{\text{max}}} $ fall on a straight line, indicating an exponential convergence of the form

 

Figure 4. Convergence analysis for the calculation of $ {\tilde \omega _1} $ using spherical wave functions showing, as a function of $ {l_{\text{max}}} $, the difference $ {D_{1 + }}({l_{\text{max}}}) $ between results obtained using $ {l_{\text{max}}} $ and $ {l_{\text{max}}} + 1 $. Red circles and blue squares correspond to real and imaginary parts of the difference, respectively. Solid and dashed lines show the corresponding fits to the last two data points of each set.

Download Full Size | PDF

With this assumption, we can use fits to the data points in Fig. 4 to estimate the true value of $ {\tilde \omega _1} $ to an accuracy of approximately one part in a billion,

2.2b. Calculations Using a Boundary Element Method

In the BEM approach of Refs. [119–121,129], the electric and magnetic vector potentials $ \phi ({\bf r}) $ and $ {\bf A}({\bf r}) $ are cast in terms of surface charges and currents, $ \sigma ({\bf r}) $ and $ {\bf h}({\bf r}) $, inside and outside the material boundaries. The material boundary itself is then discretized in surface elements, and we can formally write the corresponding discretization of the surface charges and currents as

For a given calculation mesh approximating the resonator geometry, we can set up Eq. (13) to map out the complex frequency plane as in Fig. 3 and find the QNMs at the positions where the eigenvalues of the operator $ {\bf \Sigma }(\omega ) $ vanish. The discretization introduces numerical errors via the discrete approximation to the surface charges and currents and via the approximation of the spherical surface through piecewise flat triangles. As in the case of the VIE method, we can investigate the change in resonance frequency as a function of the fineness of the mesh to assess the convergence properties of the BEM method. Because of the high estimated accuracy of the VIE method, we can use the result in Eq. (11) as a reference value for calculating the error. Figure 5 shows the relative error as a function of average side length in a double logarithmic plot. The calculated relative errors in Fig. 5 are consistent with a polynomial convergence in which the calculated value depends on the triangle side length $ h $ as

 

Figure 5. Convergence analysis for the BEM calculation of $ {\tilde \omega _1} $ showing the relative error $ {{\cal E}_{\text{BEM}}} = |{\tilde \omega _{\text{BEM}}} - {\tilde \omega _{\text{ref}}}|/|{\tilde \omega _{\text{ref}}}| $ as a function of average side length $ h $ of the triangular mesh. The reference value $ {\tilde \omega _{\text{ref}}} $ is given in Eq. (11), while $ {\tilde \omega _{\text{BEM}}} $ is the value obtained with the BEM approach. Circles indicate the calculated values, and the black line shows a fit to the last four data points. Inset shows an example of a relatively coarse mesh with ${\log }_{10} \{ h/d\} = - 0.5 $.

Download Full Size | PDF

Supplementary code. As a supplement to this Tutorial, we provide a number of MATLAB files that enable the calculation of QNMs using MNPBEM (Code 2, Ref. [118]). In particular, we provide the files necessary to reproduce the Purcell factor spectrum and the mode profiles in Fig. 2 as well as the convergence analysis in Fig. 5 (see also Appendix C.2a). Even though the examples focus on the plasmonic dimer of gold nano-spheres, they may also serve as a convenient template for analysis of QNMs in other geometries.

3. Theoretical Framework

In this section, we lay out the basic elements of QNM modeling theory for general, three-dimensional resonators in homogeneous environments and discuss how the resulting expressions relate to various formulations in the literature.

The main emphasis is on a biorthogonal approach, in which the differential equation and the radiation condition are used in a systematic way to define an adjoint problem along with an associated projection operation. In Subsection 3.1, we start by precisely defining the differential equation problem of interest, namely Maxwell’s equations with no free charges subject to the Silver–Müller radiation condition; the QNMs are defined as the solutions to this problem in the absence of sources. Subsection 3.2 shows how one can use the differential equation problem and the radiation condition as inspiration to define the adjoint QNMs and a useful projection operator with the interesting properties that the projection of one QNM onto another is zero, whereas the projection of a QNM onto itself provides a well-known normalization formula for QNMs. The projection operator can be used to expand certain electromagnetic fields, such as the electromagnetic Green tensor, in the vicinity of electromagnetic resonators, as we discuss in Subsection 3.3. A conceptually distinct but complementary approach to QNM modeling is offered by the Green tensor approach presented in Subsection 3.4, which is inherently very generally applicable because the Green tensor by definition obeys the correct radiation condition, whatever that may be. In practical calculations along such a route, one can use the normalization procedure in Subsection 3.4a, provided one can calculate the Green tensor at complex frequencies.

We note that the definition in Subsection 3.1 for the differential equation problem and its solutions is consistent with the general literature, cf. Subsection 1.3. In particular, the formulation in terms of six-component field vectors and the matrix Green tensor in Subsection 3.1a is similar to that of Refs. [29,130]. The use of a two-component formalism and a biorthogonal approach to QNM modeling in one dimension were shown in Refs. [22,24,41], respectively, see also Ref. [42]. In the presentation below, the introduction of the adjoint QNMs and the normalization procedure in Subsection 3.2 arise by extending the classical biorthogonal approach [112] by insights gained from the requirement that the QNMs obey the Silver–Müller radiation condition. The introduction of a projection operator is arguably similar in spirit to the modeling approach in Refs. [18,53], but the resulting expression is different and features an additional surface integral; the associated normalization integral was previously presented in Ref. [130]. The results for the normalization in the case of dispersive materials in Subsection 3.2b have been derived by different methods in Refs. [18,25,59,60,130], see also Ref. [66], and we note that variants of the formal expansions in Subsection 3.3, as well as the complementary analysis in Subsection 3.4, have been shown to arise by a number of different approaches and a variety of material systems [18,20,21,26,27,31,38–40,42,53,60,66,130].

3.1. Differential Equation Problem and Its Solutions

Assuming a time dependence of the form $ \exp \{ - \text{i}\omega t\} $, we shall be interested in electric and magnetic fields, $ {\bf E}({\bf r},\omega ) $ and $ {\bf H}({\bf r},\omega ) $, which solve Maxwell’s equations with electric and magnetic source currents, $ {{\bf J}_{\text{s}}}({\bf r},\omega ) $ and $ {{\bf M}_{\text{s}}}({\bf r},\omega ) $, respectively, and no free charges. In particular, we require the fields to obey the divergence equations $ \nabla \cdot [{ \epsilon _0}{ \epsilon _{\text{R}}}({\bf r}){\bf E}({\bf r},\omega )] = \nabla \cdot [{\mu _0}{\mu _{\text{R}}}({\bf r}){\bf H}({\bf r},\omega )] = 0 $, where $ { \epsilon _0} $ and $ {\mu _0} $ denote the permittivity and permeability of free space, and $ { \epsilon _{\text{R}}}({\bf r}) $ and $ {\mu _{\text{R}}}({\bf r}) $ denote the relative permittivity and permeability, which we assume, initially, to be dispersionless and absorptionless, i.e., $ { \epsilon _{\text{R}}}({\bf r}) $ and $ {\mu _{\text{R}}}({\bf r}) $ are assumed to be real. The case of dispersive and absorptive materials is treated in Subsection 3.2b, and in general we shall limit the discussion to piecewise homogeneous and passive materials without gain. At sufficiently large distances, moreover, we shall assume that $ { \epsilon _{\text{R}}}({\bf r}) $ and $ {\mu _{\text{R}}}({\bf r}) $ take on the constant values $ { \epsilon _{\text{B}}} = n_{\text{B}}^2 $ and $ {\mu _{\text{B}}} = 1 $. In addition to the divergence equations, we require the fields to obey the curl equations, which we write for the six-component electromagnetic field vector $ \underline {\bf F} ({\bf r},\omega ) = [{\bf E}({\bf r},\omega ),{\bf H}({\bf r},\omega {)]^{\text{T}}} $ in the form

$ \underline {\underline {\bf W} } = \text{diag}\{ { \epsilon _0}{ \epsilon _{\text{R}}}({\bf r}),{\mu _0}{\mu _{\text{R}}}({\bf r})\} $, and $ \underline {\bf J} (\omega ) = [{{\bf J}_{\text{s}}}({\bf r},\omega ),{{\bf M}_{\text{s}}}({\bf r},\omega {)]^{\text{T}}} $. Last, but not least, we require the fields to obey the Silver–Müller radiation condition [131,132] in the form

In general, at positions sufficiently far away, the solutions to Eqs. (15), (17), and (18) can be written in terms of outgoing waves with exponential factors of the form

In piecewise homogeneous materials with no free charges, the electric and magnetic fields are manifestly transverse at all positions except the material interfaces, where the divergence equations lead to the familiar requirement of continuity of the normal components of $ { \epsilon _{\text{R}}}({\bf r}){\bf E}({\bf r},\omega ) $ and $ {\mu _{\text{R}}}({\bf r}){\bf H}({\bf r},\omega ) $. In many practical approaches to the calculation of electromagnetic fields in problems with piecewise homogeneous material distributions, one therefore expands the unknown fields in transverse solutions to the curl equations in each of the material domains individually, and subsequently uses the boundary conditions to connect them. In a similar spirit, we shall in general focus mostly on the solutions to the curl equations and the Silver–Müller radiation condition, but with the implicit assumption that the electric and magnetic fields obey also the divergence equations and therefore the correct boundary conditions at all interfaces.

3.1a. Matrix Green Tensor

To solve the problem defined by Eqs. (15), (17), and (18), we can use a matrix Green tensor formalism by generalizing the ideas of Green functions and Green tensors, which we review in Appendix A. In this way, we find that the general solution can be calculated via integration as

3.1b. QNMs

In the absence of sources, Eqs. (15), (17), and (18) have no solutions at real frequencies, except possibly the case $ \omega = 0 $, which leads to nontrivial solutions only for longitudinal fields. Nevertheless, it is possible to find solutions at complex frequencies, which behave as outgoing waves and obey Eqs. (17) and (18) in the weaker sense $ A \to B $ if $ A/B \to 1 $. These solutions can be expanded in outgoing spherical wave functions of the form in Eq. (21), whose analytical continuations onto the real frequency axis obey Eqs. (17) and (18) in the stronger sense $ A \to B $ if $ A - B \to 0 $. These solutions are what we refer to as the QNMs. To explicitly distinguish them from general fields, we use the notation $ {\tilde {\underline {\bf F} }_n}({\bf r}) = [\tilde{\mathbf{f}}_n({\bf r}),\tilde{\mathbf{g}}_n({\bf r}{)]^{\text{T}}} $ for the six-component QNMs, in which case we can write the defining equation for the QNMs in compact form as

Since $ \underline {\underline {\bf D} } $ is real, it is clear from Eq. (25) that if $ {\tilde {\underline {\bf F} }_n}({\bf r}) $ is a QNM with resonance frequency $ {\tilde \omega _n} $, then $ \tilde {\underline {\bf F} }_n^*({\bf r}) $ is another QNM with resonance frequency $ - \tilde \omega _n^* $.

From the chosen time dependence of the form $ \exp \{ - \text{i}\omega t\} $ and the assumption of no gain in the materials, we can infer that the imaginary part of the complex QNM resonance frequency $ {\tilde \omega _n} = {\omega _n} - \text{i}{\gamma _n} $ must be negative ($ {\gamma _n} \gt 0 $) in order to correspond to an overall temporal decay of the local electromagnetic field as energy radiates away from the resonator or is absorbed in the material. This is consistent with the example calculations in Subsections 2.1 and 2.2. Indeed, as we shall see, the QNM resonance frequencies can be associated with the poles of the electromagnetic Green tensor, which must all be located in the lower half of the complex plane in order to ensure causality. Also, as we shall see, from the definition of the $ Q $ factor of a resonance as the ratio of the angular resonance frequency to the full width at half-maximum (FWHM) of the spectral response, we can calculate the $ Q $ factor pertaining to a single QNM resonance frequency as

The temporal decay of the electromagnetic fields associated with each of the QNMs has an interesting effect on the spatial variation of the QNMs. From Eq. (20), it is clear that if the frequency is complex with a negative imaginary part, then $ {\tilde {\underline {\bf F} }_n}({\bf r}) $ must increase exponentially at sufficiently large distances outside the resonators. This increase is visible, for example, in the QNM of the dielectric barrier in Fig. 1. The exponential increase of the QNMs is a real effect, in the sense that if one excites a resonator close to the resonance frequency of a QNM, the electromagnetic field in the resonator will subsequently decay exponentially in time as energy leaks to the environment. From the exponential nature of the decay, the field just outside the resonator is directly proportional to the field in the resonator; this can be seen clearly in the CMT calculations in Subsection 5.1, for example. Therefore, the exponential temporal decay for a time-harmonic field dependence is correspondingly mapped onto what appears to be an exponential growth of the spatial dependence of the field propagating away. As time goes on, the field in the resonator decays to zero, and the field outside the resonator propagates further away. In this way, the field at the wavefront becomes exponentially large relative to the field in the resonator [9], as illustrated in Fig. 6.

 

Figure 6. Sketch of the exponential decay of a field in a resonator due to radiation. At $ t = 0 $, the field is approximately confined to the resonator, illustrated by gray shading. At all times, the field is continuous across the resonator boundary, and as the field inside the resonator decays in time, the field outside propagates away (at the speed of light). In this way, the exponential temporal decay of the field inside the resonator is mapped inversely onto the spatial distribution outside, and the field at the wavefront becomes exponentially large compared to the field in the resonator.

Download Full Size | PDF

For practical problems of interest, the calculations can often, if not always, be carried out in regions close to the resonators, wherefore the exponential growth of the field does not lead to unphysical results. Indeed, consistency of the QNM approximations is restricted to regions close to the resonator, as we discuss in Subsection 4.3. As a result, the exponential growth of the fields outside the resonator should under no circumstances be interpreted in such a way that the coupling strength of, say, an electric dipole to the field in the resonator increases with increasing separation. On a related note, we remark that there is no immediate connection between the magnitude of the individual QNM fields and the electromagnetic energy density. In particular, this means that the exponential growth of the QNM fields outside the resonators cannot in any way lead to the conclusion that the fields carry an infinite amount of energy.

The complex resonance frequencies, and the associated divergent nature of the QNMs, are consequences of the fact that they appear as solutions to a non-Hermitian eigenvalue problem. For this reason, they cannot immediately be used for modeling by the well-known mathematical framework developed for Hermitian eigenvalue problems—in particular, they cannot be normalized by the norm that is often used for such problems. Instead, a slightly more general approach based on the adjoint eigenvalue problem is needed, as we discuss in the next section.

3.2. Adjoint QNMs and Normalization

As in usual approaches to field expansions, we now seek an appropriate operator for projection of a given field onto the QNMs. As a starting point, we consider the familiar inner product $ \langle {\underline {\bf F} _A}({\bf r})|{\underline {\bf F} _B}({\bf r})\rangle $ between two arbitrary electromagnetic fields $ {\underline {\bf F} _A}({\bf r}) = [{{\bf E}_A}({\bf r}),{{\bf H}_A}({\bf r}{)]^{\text{T}}} $ and $ {\underline {\bf F} _B}({\bf r}) = [{{\bf E}_B}({\bf r}),{{\bf H}_B}({\bf r}{)]^{\text{T}}} $ defined with the weight function $ \underline {\underline {\bf W} } = \text{diag}\{ { \epsilon _0}{ \epsilon _{\text{R}}}({\bf r}),{\mu _0}{\mu _{\text{R}}}({\bf r})\} $ as

Using the relation $ \nabla \cdot [{\bf A} \times {\bf B}] = {\bf B} \cdot [\nabla \times {\bf A}] - {\bf A} \cdot [\nabla \times {\bf B}] $ and the divergence theorem, one can rewrite the expression as

To determine the adjoint QNMs, we now seek to define both the associated differential operator and the corresponding radiation condition in order to ideally satisfy the relation $ \langle \tilde {\underline {\bf F} }_m^{\ddagger }({\bf r})|\underline {\underline {\bf D} }\; \underline {\bf F} ({\bf r},\omega )\rangle = \langle {\underline {\underline {\bf D} } ^{\ddagger }}\tilde {\underline {\bf F} }_m^{\ddagger }({\bf r})| \underline {\bf F} ({\bf r},\omega )\rangle $ [112]. From Eq. (30), it follows that we should define the adjoint operator as $ {\underline {\underline {\bf D} } ^{\ddagger }} = - \underline {\underline {\bf D} } $. Moreover, in the limit of large volume (and real frequencies), the fields $ {\bf E}({\bf r},\omega ) $ and $ {\bf H}({\bf r},\omega ) $ fulfill the Silver–Müller radiation condition on the boundary $ \partial V $. In this limit, therefore, one can rewrite the surface integral as

Compared to Eqs. (17) and (18), the signs have flipped, so that the analytical continuation of the Poynting vector $ \tilde{{\bf S}}_m^{\ddagger }({\bf r}) = \tilde{{\bf f}}_m^{\ddagger }({\bf r}) \times \tilde{{\bf g}}_m^{\ddagger }({\bf r}) $ points inward. Writing out the matrix equation, one can verify that if $ {\tilde {\underline {\bf F} }_n}({\bf r}) = [\tilde{\bf f}_n({\bf r}),\tilde{\bf g}_n({\bf r}{)]^{\text{T}}} $ is a solution to Eq. (25) obeying the Silver–Müller condition in Eqs. (17) and (18), then $ \tilde {\underline {\bf F} }_m^{\ddagger }({\bf r}) = [\tilde{\bf f}_m({\bf r}), - \tilde{\bf g}_m ({\bf r}{)]^{\text{T}}} $ is a solution to the equation

The adjoint QNMs were chosen so as to make the surface integral vanish in the limit of large volume (and real frequencies). At finite distances (and complex frequencies), the surface integral is nonzero in general, and we must keep it. For the special case of $ \omega = {\tilde \omega _m} $, it follows from Eq. (36) that the surface integral vanishes identically, but so does the denominator in the second term of Eq. (39). To investigate this limit, we can expand the electromagnetic field $ \underline {\bf F} ({\bf r},\omega ) $ to first order as

This relation, although valid as written, is slightly deceptive, in that it might suggest that one can directly insert any sum of QNMs in Eq. (39) to get the projection onto the QNM of interest. Such an approach will not work in general, however, because of the frequency dependence of the second term in Eq. (39). In practice, therefore, we find meaningful expressions only when working with solutions to the wave equation (which is not the case for general superpositions of QNMs, because of the explicit frequency dependence in Eq. (25)).

For the special class of electromagnetic fields $ \underline {\bf F} ({\bf r},\omega ) $, which solve Eqs. (15), (17), and (18), the operation $ \langle \langle \tilde {\underline {\bf F} }_m^{\ddagger }({\bf r})|\underline {\bf F} ({\bf r},\omega )\rangle \rangle $ does provide the projection onto the QNM $ {\tilde {\underline {\bf F} }_m}({\bf r}) $, as we show in Subsection 3.3. The necessity of the fields $ \underline {\bf F} ({\bf r},\omega ) $ to obey the Silver–Müller radiation condition lies in the fact that only if this is true does the surface integral vanish in the limit $ \omega \to {\tilde \omega _m} $. In particular, Eq. (37) is valid for any $ \underline {\bf F} ({\bf r},\omega ) $ that solves Eq. (15), but only if $ \underline {\bf F} ({\bf r},\omega ) $ behaves as the proper analytical continuation of the QNM is it true that $ \underline {\bf F} ({\bf r},\omega ) \to {\tilde {\underline {\bf F} }_m}({\bf r}) $ as $ \omega \to {\tilde \omega _m} $. Moreover, as shown in Appendix C.4, the sum of the two integrals in Eq. (39) is independent of size and shape of the volume $ V $, as long as it contains the electromagnetic resonator. If the integral is extended to infinity using a complex coordinate transformation, for example, by the use of PMLs, the surface integral vanishes, and the projection reduces to the exact expression for the projection due to Sauvan et al. [18].

In one dimension, the radiation condition turns into an algebraic boundary condition, and the second term in Eq. (39) vanishes identically. The first term then defines a projection of an arbitrary sum of QNMs onto $ {\tilde {\underline {\bf F} }_m}({\bf r}) $, as suggested from Eq. (43).

3.2a. Degenerate Modes

For degenerate QNMs with $ {\tilde \omega _m} = {\tilde \omega _n} $, but $ {\tilde {\underline {\bf F} }_m}({\bf r}) \ne {\tilde {\underline {\bf F} }_n}({\bf r}) $, it is possible that the modes are not immediately orthogonal under the operation in Eq. (39). In such cases, as long as the individual QNMs are normalizable, one can always define a new set of QNMs via the Gram–Schmidt orthogonalization process. Moreover, in certain cases of degenerate QNMs, Eq. (42) may not be immediately useful for normalization. In the case of a spherical resonator, for example, the QNMs can be chosen to have an azimuthal dependence of the form $ \exp \{ \pm \text{i}\varphi \} $ for which the integral in Eq. (42) vanishes. Any small perturbation of the circular shape, however, will break the degeneracy and result in QNMs with azimuthal dependencies of the form $ \exp \{ \text{i}\varphi \} \pm \exp \{ - \text{i}\varphi \} $. In this case, therefore, one can define the QNMs by these linear combinations, for which Eq. (39) is immediately applicable [61]. In the particular case of spheres, for which the QNMs can always be written as products of the form $ {\underline {\bf F} _m}({\bf r}) = {{\bf R}_m}(r){{\bf P}_{lm}}(\theta ,\varphi ) $, where $ {{\bf P}_{lm}}(\theta ,\varphi ) $ is any of the vector spherical harmonics of order $ (l,m) $, one can also define the adjoint modes via complex conjugation of the angular dependence only, as was done in Refs. [26,27]. Throughout this work, however, we shall exclusively use the definition $ \tilde {\underline {\bf F} }_m^{\ddagger }({\bf r}) = [\tilde{\bf f}_m({\bf r}), - \tilde{\bf g}_m({\bf r}{)]^{\text{T}}} $ and assume that any problems arising from degeneracies can be handled by defining the QNMs as a suitable linear combination, as described above.

Under very special circumstances, the complex spectrum of QNM frequencies may feature a so-called exceptional point, where two or more QNMs coalesce so that $ {\tilde \omega _m} = {\tilde \omega _n} $ and $ {\tilde {\underline {\bf F} }_m}({\bf r}) = {\tilde {\underline {\bf F} }_n}({\bf r}) $, for $ m \ne n $. This very interesting phenomenon has been found to occur by careful symmetry breaking of ring resonators, for example [135]. In principle, the theory laid out in this Tutorial is not immediately suited to handle such cases. In most practical applications, however, the system is never exactly at the exceptional point, and therefore one can in principle apply the theory as is, provided the calculations are sufficiently accurate to distinguish the (small) differences in the QNMs.

3.2b. Dispersive and Absorptive Materials

The biorthogonal framework can be immediately extended to the technologically interesting case of resonators made from dispersive and absorptive materials by use of auxiliary fields governing the material response [59,136]. In the case of a Drude material response as in Eq. (1), we can write the general relative permittivity distribution as

As described in detail in Appendix C.3, one can follow an approach completely analogous to the one used for the dispersionless case above, and define generalized QNM fields $ {\tilde {\underline {\bf F} }_m}({\bf r}) = [\tilde{{\bf f}}_m({\bf r}),\tilde{{\bf g}}_m({\bf r}),\tilde{{\bf j}}_m({\bf r}{)]^{\text{T}}} $ as the solutions to an equation of the same form as Eq. (25) along with a corresponding projection operator

In the limit $ \omega \to {\tilde \omega _n} $, we find that the normalization can be written as in Eq. (42) with the substitution $ { \epsilon _{\text{R}}}({\bf r}) \to \eta ({\bf r},{\tilde \omega _m}) $, where

3.2c. Adjoint QNMs and Normalization for the Dielectric Barrier

In one dimension, the radiation condition turns into an algebraic boundary condition, cf. Eq. (6). Therefore, if we choose the integration volume for the normalization of the QNM of the dielectric barrier to be $ |x| \le L/2 $, the second term in Eq. (39) vanishes identically, provided the adjoint QNMs obey the boundary conditions

As in the general case, we can now see that if $ \tilde{\text{f}}_m(x) $ and $ \tilde{\text{g}}_m(x) $ are electric and magnetic field QNMs, then $ \tilde{\text{f}}_m^{\ddagger }(x) = \tilde{\text{f}}_m(x) $ and $ \tilde{\text{g}}_m^{\ddagger }(x) = - \tilde{\text{g}}_m(x) $ are the associated adjoint electric and magnetic field QNMs. In the absence of the second term in Eq. (39), we can immediately take the limit $ \omega \to {\tilde \omega _m} $. In this way, we can write the integral for the QNM normalization as

Note that even though we chose the integration to be $ |x| \le L/2 $ for convenience, the value is independent of this choice, as long as the integration boundaries are beyond the extent of the barrier. In Eq. (50), this property is a result of the fact that the integrand vanishes identically outside the resonator. Inserting the explicit form of $ \tilde{\text{f}}_m (x) $ and $ \tilde{\text{g}}_m(x) $, the integral simplifies substantially, and we find

3.3. Formal Expansions in Terms of QNMs

We now first take the point of view that a given field can be expanded in QNMs within a volume $ V $ enclosing the resonator, in which case the projection operator arises naturally and provides us with the expansion coefficients. Subsequently, we will worry about the validity of the expansion, i.e., the convergence and consistency of the formal expansion.

3.3a. Formal Expansion of a General Electromagnetic Field

We consider a general electromagnetic field $ \underline {\bf F} ({\bf r},\omega ) $, which solves Eqs. (15), (17), and (18). Multiplying Eq. (15) from the left with $ \tilde {\underline {\bf F} }_m^{\ddagger }({\bf r})\underline {\underline {\bf W} } $ and integrating over a volume $ V $ containing the resonator and all sources, we find that

The peaks in Fig. 2 suggest that the QNMs are associated with poles in the electromagnetic response of the resonator. We will therefore assume that the electromagnetic fields of interest can be represented by the series

 

Figure 7. Example integration contours in the complex frequency plane for the projection of a general electromagnetic field onto the QNMs of the dielectric barrier in one dimension of width $ L $. An initial integration contour around the frequency $ \omega $ (top panel) is enlarged into a contour $ {\Gamma _N} $ while in the process excluding $ N $ poles at all QNM resonance frequencies inside by integration around them, as illustrated for the case of $ N = 4 $ (center panel). Arrows indicate the direction of integration. All contributions from the vertical dashed lines connecting $ {\Gamma _N} $ with the small circles around the poles at the QNM resonance frequencies vanish, because they are integrated in both directions. If the integral along $ {\Gamma _N} $ vanishes in the limit of infinitely large radius, the projection is given in terms of contributions from the QNMs only (bottom panel).

Download Full Size | PDF

Next, we express $ \underline {\bf F} ({\bf r},z) $ in the numerator of Eq. (54) using Eq. (53), but we explicitly evaluate $ {a_n}(z) $ at $ z = \omega $, which is warranted since for analytic functions $ \text{f}(\omega ) $ and $ g(\omega ) $, we have

With this procedure, the left-hand side of Eq. (52) can be rewritten as

The additional contribution from the integral along the outer contour $ {\Gamma _N} $ vanishes because of the functional form of the integrand. Using the orthogonality in Eq. (43), we conclude that

Finally, returning to Eq. (52), we find that if the electromagnetic field can be expanded as in Eq. (53), then the expansion coefficients are given as

From the expansion in Eq. (53) or the projection in Eq. (58), we can now appreciate that for a constant source, the electromagnetic response due to each of the QNMs is a Lorentzian centered on $ \omega = {\tilde \omega _m} $ and with a FWHM of $ 2{\gamma _m} $, thus justifying the definition in Eq. (26). Note, however, that in general the contributions to the spectrum from different QNMs interfere, and depending on the particular resonator and the excitation conditions, this may lead to substantial deviations from Lorentzian line shapes in practice.

It is instructive to consider the problem of the source-free electromagnetic field. Setting $ \underline {\bf J} ({\bf r},\omega ) = 0 $, it follows from Eq. (59) that there are no nontrivial solutions to the source-free wave equation of the form in Eq. (53) with $ {a_n}(\omega ) \ne 0 $. Instead, the solutions are of the form $ \underline {\bf F} ({\bf r},\omega ) = {\tilde {\underline {\bf F} }_n}({\bf r}) $, i.e., they are QNMs. In this case, analytic continuation of the QNMs implies that we must have $ {a_n}(\omega ) = (\omega - {\tilde \omega _n})(1 + {\cal O}(\omega )) $.

3.3b. Formal Expansion of the Matrix Green Tensor

As a special case, the QNMs can be used for expansion of the matrix Green tensor. Because of the transverse nature of the QNMs in regions with constant permittivity, the expansions are limited to the transverse part of the Green tensor, as further discussed in Subsection 4.1. To proceed, we start from Eq. (23) and assume that the matrix Green tensor can be expanded as

Inserting in Eq. (60), the matrix Green tensor expansion then takes the form

We note that because of the symmetric property of the QNM resonance frequency spectrum, as illustrated in Fig. 3, it follows from Eq. (63) that the electric field Green tensor satisfies the so-called crossing relation,

 

Figure 8. Illustration of the QNM approximation to the electric field Green function for the dielectric barrier in one dimension from Subsection 1.5a. Top: QNM Green function approximation $ {G_{500}}(x,x^{\prime},\omega ) $ as a function of $ x $, with $x^{\prime} = 0 $, $ {\omega _0} = 4L/\text{c} $, and in units of $ {G_0} = \text{c}/2{n_{\text{R}}}\omega $. Red solid and blue dashed curves show the real and imaginary parts, respectively, and the black curve shows the absolute value. The gray shading indicates the extent of the dielectric barrier. Bottom: relative error $ {{\cal E}_{500}}(x,x^{\prime},{\omega _0}) $.

Download Full Size | PDF

3.3c. Green Function Expansion for the Dielectric Barrier

As a first test of the expansion in Eq. (63), we return to the example of the dielectric barrier from Subsection 2.1, and consider a QNM approximation to the electric field Green function of the form

 

Figure 9. Convergence analysis for the QNM Green function expansion inside the dielectric barrier. Red circles show the relative integrated error $ {{\cal I}_N}(0,{\omega _0}) $ with fixed $ x^{\prime} = 0 $ in the region $ - L/4 \lt x \lt L/4 $ around the center of the barrier, and blue squares show the relative error $ {{\cal E}_N}(0,0,{\omega _0}) $ at $ x = x^{\prime} = 0 $ in the center of the barrier. For both data sets, $ {\omega _0} = 4L/\text{c} $, and the solid and dashed lines show linear fits to the last three data points.

Download Full Size | PDF

As a first glimpse at the convergence properties of the QNM Green function approximation in Eq. (65), Fig. 9 shows a double logarithmic plot of the relative error in Eq. (66), as well as the integrated error

Supplementary code. The supplementary code Code 1, Ref. [117] provides the files necessary to reproduce Fig. 9, and we encourage interested readers to investigate the convergence of the QNM Green function approximation when changing the number of terms in Eq. (65) or study the nontrivial divergence of the QNM Green function approximation as $ x \to \infty $.

3.4. QNMs as the Residues of the Matrix Green Tensor

To lay the foundation for the convergence analysis in Section 4, and as a complementary QNM modeling approach to that taken in Subsection 3.3, it is illustrative to see how the QNMs appear as the residues of the matrix Green tensor. The starting point is that the Green tensor has a number of poles, all of which are located in the lower half of the complex plane. Close to a pole $ {\Omega _n} $, we assume that the Green tensor can in general be approximated as

It follows from Eq. (69) and the inherited radiation condition that each column of $ {\underline {\underline {\boldsymbol \rho } } _n}({\bf r},{{\bf r}^\prime}) $ solves the defining equations for the QNMs in Eqs. (17), (18), and (25). Therefore, the complex QNM resonance frequencies are the poles of the Green tensor, $ {\Omega _n} = {\tilde \omega _n} $, and $ {\underline {\underline {\boldsymbol \rho } } _n}({\bf r},{{\bf r}^\prime}) $ is proportional to the QNM $ {\tilde {\underline {\bf F} }_n}({\bf r}) $, as was tacitly assumed in Eq. (60). In the vicinity of $ \omega \approx {\tilde \omega _n} $, we write this as

Relaxing now the condition $ {\bf r} \ne {{\bf r}^\prime} $, we use Eqs. (68) and (70) in Eq. (23), multiply from the left with $ \tilde {\underline {\bf F} }_n^{\ddagger }({\bf r})\underline {\underline {\bf W} } $, and integrate over the volume $ V $ to find that

Last, by taking the limit $ \omega \to {\tilde \omega _n} $, we find that the residues of the matrix Green tensor can be expressed as

3.4a. Alternative Normalization Procedure

From a numerical point of view, especially with integral equation techniques, the use of a volume integral for normalization is not ideal, since one must necessarily evaluate the field at numerous positions within the integration domain. Bai et al. [65] suggested a combined calculation and normalization approach for QNMs, which is useful when one has access to any of the components of the matrix Green tensor at complex frequencies. In a slightly different formulation, but exploiting the same idea, we can derive a variant of this normalization evaluation by noting that if Eq. (72) holds, then the inverse of the normalization integral can be inferred from the residue of any of the components of $ \underline {\underline {\bf G} } ({\bf r},{{\bf r}^\prime},\omega ) $ at $ \omega = {\tilde \omega _m} $. In particular, for $ \alpha ,\beta \in \{ x,y,z\} $, one can use the $ \alpha ,\beta $ component of the electric field Green tensor along with the $ \alpha $ and $ \beta $ components of the electric field QNM at the positions $ {\bf r} $ and $ {{\bf r}^\prime} $, respectively, to calculate the inverse normalization integral as

A particularly interesting property of Eq. (74) is that it can be handled effectively by the trapezoidal rule because of the periodic nature of the integral [137]. In practice, we can always find a circle of radius $ R $ centered on the QNM resonance frequency $ {\tilde \omega _m} $ and surrounding only this one pole. Parameterizing the curve in terms of the angle $ \theta $ as $ z(\theta ) = {\tilde \omega _m} + R\exp \{ \text{i}\theta \} $, we write the integral as

3.4b. Alternative Normalization for the Dielectric Barrier

Using the analytical expression for the electric field Green function in Appendix C.1b, we can evaluate the right-hand side of the (inverse) normalization integral in Eq. (74) by explicitly calculating the residue of the Green function at $ z = {\tilde \omega _m} $. Using the eigenfunctions in Eq. (2), it takes the value

 

Figure 10. Convergence analysis for the numerical calculation of the inverse normalization integral for the QNM of the dielectric barrier with $ m = 4 $. The figure shows the relative error $ {\cal E}(N) = |\langle \langle {\tilde {\underline {\bf F} }_4}(x)|{\tilde {\underline {\bf F} }_4}(x)\rangle \rangle _N^{ - 1} - \langle \langle {\tilde {\underline {\bf F} }_4}(x)|{\tilde {\underline {\bf F} }_4}(x)\rangle \rangle _{\text{ref}}^{ - 1}|/|\langle \langle {\tilde {\underline {\bf F} }_4}(x)|{\tilde {\underline {\bf F} }_4}(x)\rangle \rangle _{\text{ref}}^{ - 1}| $ (red circles) and the difference $ {D_{ + 2}}(N) $ (blue squares) between results obtained using $ N $ and $ N + 2 $ points in the trapezoidal sum. Solid and dashed lines show the corresponding fits to the last two data points of each set.

Download Full Size | PDF

The logarithm of $ {D_{ + 2}}(N) $ is also shown in Fig. 11, and the points fall on a straight line with a slope similar to that of the data for the relative error. Using the analysis in Appendix B, we can estimate the value of the inverse normalization integral to be

 

Figure 11. Convergence analysis for the calculation of the inverse normalization integral in the case of the plasmonic dimer. The figure shows, as a function of $ {l_{\text{max}}} $, the difference $ {D_{1 + }}({l_{\text{max}}}) $ between results obtained using $ {l_{\text{max}}} $ and $ {l_{\text{max}}} + 1 $. Red circles and blue squares correspond to real and imaginary parts of the difference, respectively. Solid and dashed lines show the corresponding fits to the last two data points of each set.

Download Full Size | PDF

3.4c. Normalization for the Plasmonic Dimer

We have no reference calculations of the normalization integral for the QNMs of the plasmonic dimer. Moreover, we have no analytical expressions for the QNMs or the electric field Green tensor, so both will have to be calculated numerically, and any numerical error will ultimately limit the accuracy of the calculated inverse norm.

The convergence properties of the trapezoidal rule approximation to the inverse normalization integral in Eq. (76) is almost independent of the accuracy of the integrand. Therefore, we can perform calculations with a relatively low accuracy of the integrand using $ {l_{\text{max}}} = 2 $ to find the number of elements $ N $ necessary for a given accuracy of the integral. Subsequently, using this value of $ N $, we then vary $ {l_{\text{max}}} $ to estimate the true value of the inverse normalization integral as done also for the complex QNM frequency in Subsection 2.2. To this end, we define

4. Convergence and Consistency

To answer the question of convergence and consistency of the formal expansions, we follow a line of argumentation similar to that of Ref. [2] for solutions to the Schrödinger equation and Refs. [20,21] for optical cavities. We use the word “convergence” to describe the property of certain series, that they tend to a finite value as the number of terms is increased. The more restrictive requirement of “consistency” is reserved for the cases where the finite limiting value is also the correct function value of the underlying continuous problem that one hopes to approximate, typically via some sort of discretization. Appendix B.1 discusses these properties in the context of general numerical calculations, and within the present framework, the use of $ N $ QNMs can be understood simply as a specific choice of basis function set, with the parameter $ h = 1/N $ controlling the accuracy.

In practice, of course, we are mostly interested in the cases where the QNM expansions are consistent. As we have already seen in Fig. 8, this is not always the case, and the QNM expansion may even be convergent at some points without being consistent. Even though the answer to the general question of convergence and consistency of QNM expansions appears to be out of reach at this point, we can make some useful progress by carefully analyzing the structure of the appropriate Green tensor as in Refs. [2,20,21] and combine this insight with numerical investigations. In this way, we can introduce so-called regions of consistency, which generally depend on the materials as well as the geometries of the individual resonators, and by exploiting the results of the analyses, we can device a prescription to formally extend these regions.

For the one-dimensional resonator example, we can analyze the region of consistency exactly and show that the QNM expansion of the electric field Green function and other expansions that can be derived from the Green function are consistent only for positions $ x $ and $ x^{\prime}$ strictly within the resonator. Nevertheless, for positions $ x^{\prime} $ within the resonator, the region of consistency for $ x $ extends beyond the boundary, in agreement with Fig. 8. To handle problems with general resonators, we describe a numerical estimate of the region of consistency and illustrate the practical use of the estimation method by full three-dimensional calculations for the plasmonic dimer of gold nanospheres.

As noted above, the argumentation in Subsections 4.1 and 4.2 for analyzing the convergence and consistency of the formal QNM expansions is similar to that of Refs. [2,20,21] (see also Refs. [32,39]). A practically identical integral approach is used in the so-called Riesz projection in Refs. [67–69], but with the difference that the extent of the outer contour is kept at a finite radius. The regions of consistency to be introduced in Subsection 4.3 arguably represent a natural generalization of the results in Refs. [20,21].

4.1. Convergence and Consistency of the QNM Green Tensor Expansions

To analyze the question of convergence and consistency of the QNM Green tensor expansions, we will employ the same calculation method of complex contour integration that we used in deriving the formal expansion of a general electromagnetic field in Subsection 3.3a. The starting point for the analysis is the fact that for any $ \omega \ne 0 $ at which $ \underline {\underline {\bf G} } ({\bf r},{{\bf r}^\prime},\omega ) $ is analytic, we can write the product of $ \omega $ and the matrix Green tensor as

Assuming, for now, that $ z\;\underline{\underline {\bf G} } ({\bf r},{{\bf r}^\prime},z) $ has no pole at $ z = 0 $, we can use the expression for the residue of the Green tensor in Eq. (72) to find that

If the integral along the contour $ {\Gamma _N} $ does not vanish, but rather tends to a finite value, we conclude that the QNM expansion is convergent, but not consistent, since the finite value of the integral is exactly the difference between the QNM expansion and the correct value of the Green tensor on the left-hand side of Eq. (84). If the Green tensor has one or more branch cuts, one must deform the contour by carefully excluding the lines in the complex plane at which the integrand is not analytic [34,39,44,52,76]. This will lead to additional integral terms in Eq. (84), and the convergence and possible consistency of the QNM expansion will then depend also on whether or not these integrals tend to a finite value or maybe even vanish. The possible additional integral terms will in general interfere with the part of the spectrum due to the QNMs, and this may lead to non-Lorentzian line shapes [76].

4.1a. Alternative Expression for the Green Tensor

If, instead of Eq. (82), we start with the expression

If we now assume (as is sometimes the case) that $ \underline {\underline {\bf G} } ({\bf r},{{\bf r}^\prime},\omega ) $ has a simple pole at $ \omega = 0 $, we can use Eq. (72) along with the same arguments as above to conclude that if $ \underline {\underline {\bf G} } ({\bf r},{{\bf r}^\prime},z) $ vanishes sufficiently fast in the limit $ |z| \to \infty $, then the matrix Green tensor can be expressed as

4.1b. Electric Field Green Function in One Dimension

In one dimension, the electric field Green function in general is the sum of a background Green function and a scattering term, as discussed in Appendix A.1. In cases where the scattering term has no pole at $ \omega = 0 $, the residue of the total Green function in this point is given by the residue of the background Green function, which is $ {a_0} = \text{ic}/2{n_{\text{B}}} $, cf. Eq. (A10) in Appendix A.1a. By direct application of Eqs. (84) and (87), and assuming that the Green function $ G(x,x^{\prime},z) $ vanishes sufficiently fast in the limit $ |z| \to \infty $, we find

It is illustrative to consider what happens if one substitutes the scattered part of the Green function, $ {G_{\text{scat}}}(x,x^{\prime},\omega ) $, for $ \underline {\underline {\bf G} } ({\bf r},{{\bf r}^\prime},\omega ) $ in Eq. (85). In this case, one can carry through the same arguments to find that if $ {G_{\text{scat}}}(x,x^{\prime},z) $ vanishes sufficiently fast in the limit $ z \to \infty $, then

4.1c. Electric Field Green Tensor in Three Dimensions

In three dimensions, the electric field Green tensor can be split in transverse and longitudinal parts as discussed in Appendix A.2. By definition of the QNMs as solutions to the source-free Maxwell equations, they are manifestly transverse at all positions except the interfaces between the piecewise homogeneous materials making up the resonators. By directly applying the divergence to either of the QNM approximations in Eqs. (82) or (85) at positions other than the interfaces, all of the terms in the expansions vanish, and we conclude that the QNMs therefore can lead to an expansion for the transverse part of the response only. To obtain the full Green tensor, one must add the longitudinal part.

Assuming that the transverse part of the Green tensor has no pole at $ \omega = 0 $, we can find a QNM expansion for it by substituting $ {\bf G}_ \bot ^{\text{EE}}({\bf r},{{\bf r}^\prime},z) $ for $ \underline {\underline {\bf G} } ({\bf r},{{\bf r}^\prime},z) $ in either Eq. (82) or (85) to find that if the integral along the contour $ {\Gamma _N} $ vanishes, the transverse part of the electric field Green tensor can be written in terms of the QNMs as

The fact that the QNM expansion will provide only the transverse part of the Green tensor was explicitly pointed out in Refs. [26,27], as well as Ref. [40], which discusses the need to include also the longitudinal modes. The second equality was also discussed in Ref. [33].

4.2. Consistency of the General Field

Any solution to Eqs. (15), (17), and (18) can be written in terms of the matrix Green tensor as in Eq. (22). Starting from this equation, and rewriting the function $ \omega \underline {\underline {\bf G} } ({\bf r},{{\bf r}^\prime},\omega ) $ as in Eq. (82), we find

In this way, the expansion of the general field inherits the convergence properties of the Green tensor, and we conclude that the formal expansion in Eq. (53) is consistent if $ z\underline {\underline {\bf G} } ({\bf r},{{\bf r}^\prime},z) \to 0 $ for $ |z| \to \infty $.

4.3. Regions of Consistency

From the discussion in Subsections 4.1 and 4.2, it is clear that a sufficient condition for consistency of the QNM expansion is for the matrix Green tensor $ \underline {\underline {\bf G} } ({\bf r},{{\bf r}^\prime},z) $ to vanish sufficiently fast in the limit $ |z| \to \infty $ along all directions in the complex frequency plane. It follows that by investigating this limit for all pairs of positions $ {\bf r} $ and $ {{\bf r}^\prime} $, one can in principle trace out the border of a region in space where the expansion is consistent. We shall refer to this as the “region of consistency,” and in practical calculations, we are often interested in the electric field Green tensor, wherefore we shall focus on this in the following analysis. The behavior of the Green tensor is almost always exponential when the frequency is varied along lines parallel to the negative imaginary axis. Indeed, in the homogeneous background, it is proportional to $ \exp \{ \text{i}z|{\bf r} - {{\bf r}^\prime}|/\text{c}\} $ and thus tends to zero exponentially as $ |z| \to \infty $ in the upper half of the complex plane. The scattered part of the Green tensor represents the scattering of the electromagnetic field due to the presence of material in an otherwise homogeneous background. From a physical point of view, therefore, it is clear that also the scattered part will vanish as $ |z| \to \infty $ in the upper part of the complex plane, and we shall focus only on the frequencies in the lower half, where the behavior is nontrivial.

For problems in one dimension [20], and for the case of spherically symmetric resonators [21], Leung et al. have used the WKB approximation to assess the limiting behavior of the Green tensor. The conclusion is that the QNM expansions are consistent when both positions $ {\bf r} $ and $ {{\bf r}^\prime} $ are inside a boundary set by the outermost spatial discontinuity in the function (or even the derivative of the function) that describes the resonator material. These findings are fully in line with very convincing direct evaluations of the associated sum closure relations [20,21] and higher-order advanced perturbation theory [21,23,27] as well as resonant state expansions [31,33,40].

With certain numerical calculation methods, such as the integral equation method in Refs. [46,138], it is possible to investigate this limit numerically by directly calculating the Green tensor at various positions in the complex plane along the negative imaginary axis. While such an analysis cannot constitute a regular proof of consistency, it can provide a convincing picture of the region of consistency for general resonators. In particular, it can be used in cases where there is no clear definition of the resonator boundary, such as in the case of the plasmonic dimer in Subsection 1.5. As we shall see, for certain choices of $ {{\bf r}^\prime} $, the region of consistency for $ {\bf r} $ may extend beyond the volume of the resonator material. What we have found consistently, however, is that the QNM expansion of the Green tensor seems to be consistent only when at least one of the spatial arguments is within the resonator material.

4.3a. Formally Extending the Region of Consistency

In cases where the region of consistency is bounded by the outermost discontinuity in the permittivity distribution, at a distance $ {x_A} $ from the origin, say, one can immediately devise a new geometry with an arbitrarily small, but discontinuous, perturbation $ \Delta \epsilon $ of the permittivity distribution at a distance of $ {x_B} \gt {x_A} $, as illustrated in Fig. 12. By identical arguments as those leading to the conclusion that the QNM expansion is consistent for $ |x| \lt {x_A} $ in the original geometry, one can appreciate that the QNM expansion will be consistent for $ - {x_A} \lt x \lt {x_B} $ in the new geometry, cf. Fig. 12.

 

Figure 12. Sketch of the material distribution for a dielectric barrier with permittivity $ { \epsilon _{\text{R}}} $ for $ - {x_A} \lt x \lt {x_A} $ and a small additional permittivity perturbation $ \Delta \varepsilon $. For any $ \Delta \epsilon \gt 0 $, the region of consistency extends to the outermost boundary of the perturbation at $ x = {x_B} $.

Download Full Size | PDF

For a given QNM approximation with a fixed number of terms $ N $, one can choose a $ \Delta \epsilon \gt 0 $ sufficiently small so that the changes to the QNM approximation are arbitrarily small. It would then seem that the same QNM approximation, which was consistent only for $ |x| \lt {x_A} $ in the original geometry, is now consistent also for $ {x_A} \lt x \lt {x_B} $. The resolution to this puzzle is the fact that even if the difference between the two QNM approximations, each with $ N $ terms, can be made arbitrarily small, the difference will not remain small if more terms are added. Indeed, only one of them will tend to the correct value of the field as the number of terms is increased. This observation, therefore, is mostly of formal relevance. In practical calculations, for example, one is often interested in single- or few-QNM approximations at positions close to the resonator, yet possibly outside the formal region of consistency. Even if a direct comparison to reference calculations shows a convincing approximation, the truncation of a series that formally does not converge to the correct value is arguably problematic. By extending the region of consistency as suggested above, one can immediately remove this concern by formally considering a different material system, the response of which can be chosen to be arbitrarily close to that of the original system over the bandwidth of interest.

4.3b. Region of Consistency for the Dielectric Barrier

The electric field Green function for the dielectric slab can be written in closed form, as detailed in Appendix C.1b, and one can therefore analytically investigate the behavior of $ G(x,x^{\prime},z) $ as $ |z| \to \infty $. In the general case, however, it is difficult to investigate the limit analytically. To illustrate the viability of a numerical approach, we use a one-dimensional version of the integral equation method in Refs. [46,138] to investigate the limiting behavior of the Green function along the negative imaginary frequency axis, $ \omega = - \text{i}\xi $, for increasing values of $ \xi $. Figure 13 shows the magnitude of the Green function for fixed $ x^{\prime} = 0 $ and three different values of $ x $ to the left of the resonator. From the logarithmic scale, it is evident that the magnitude of the Green function varies exponentially on the curve oriented downwards in the complex plane. Moreover, for $ x $ sufficiently close to the resonator, there is a qualitative change in behavior so that instead of growing exponentially, the magnitude falls off exponentially.

 

Figure 13. Top: magnitude of the electric field Green function, in units of $ {G_0}( - \text{i}\xi ) = \text{c}/2\xi $, for $ \omega = - \text{i}\xi $ plotted as a function of $ \xi $. The position $ x^{\prime} = 0 $ is held fixed at the center of the resonator of length $ L $, and $ x $ is set equal to $ x = - 3L $ (red), $ x = - 2L $ (green), and $ x = - L $ (blue). Bottom: geometry of the dielectric barrier in one dimension showing the permittivity as a function of position (thin black curve). The gray and hatched shadings indicate the regions of consistency for $ x^{\prime} = 0 $ (gray) and $ x^{\prime} = L/4 $ (hatched), as indicated by the vertical lines.

Download Full Size | PDF

The change in behavior occurs at positions slightly to the left of $ x = - 2L $. From a detailed analysis of the Green function, we find that it vanishes in the limit $ |z| \to \infty $ if

To set up a practical calculation method for estimating the region of consistency, we note that the slope of the curves in Fig. 13 is proportional to the distance of $ x $ from the boundary. Therefore, we can estimate the region of consistency by plotting the logarithm of the magnitude of the Green function when varying the real space position and for different values of $ \omega = - \text{i}\xi $. This results in a number of straight lines with different slopes, and the crossing points will approximate the boundary of the region of consistency, as illustrated in Fig. 14. The approximation becomes better when using lines pertaining to larger negative imaginary frequencies $ \xi $. Indeed, calculating the intersection between lines calculated using $ \xi L/2\pi \text{c} $ and $ \xi L/2\pi \text{c} + 0.25 $, the relative error becomes exponentially smaller with increasing $ \xi $, as seen in the bottom panel of Fig. 14. From a practical point of view, the fact that the result of this purely numerical investigation of the Green function agrees with the analytical result gives us confidence that we can apply a similar approach for general structures where no analytical results are available. In particular, we will apply this method to estimate the region of consistency of the plasmonic dimer in Subection 4.3c.

 

Figure 14. Top: magnitude of the electric field Green function as a function of position for fixed frequencies along the negative imaginary axis, $ \omega = - \text{i}\xi $, with $ \xi L/2\pi \text{c} = 0.5 $ (red), $ \xi L/2\pi \text{c} = 1 $ (green), and $ \xi L/2\pi \text{c} = 1.5 $ (blue). Bottom: relative error $ {{\cal E}_{\text{RoC}}}(\gamma ) $ when estimating the boundary of the region of consistency as the intersection of lines as in the top panel calculated using $ \xi L/2\pi \text{c} $ and $ \xi L/2\pi \text{c} + 0.25 $. The black line shows a linear fit to the data points.

Download Full Size | PDF

Extending the region of consistency for the dielectric barrier. To illustrate the ideas put forward in Subsection 4.3a, we consider now the addition of a small (constant and real valued) perturbation to the background permittivity of width $ \Delta X $ and centered on $ {x_0} = 3L $. Using the Dyson equation, we can write the total electric field Green function of the dielectric barrier and the additional perturbation as discussed in Appendix A.1b. To appreciate that the total Green function may vanish for certain choices of $ x $ and $ x^{\prime} $ and increasing values of $ \xi $, we consider the case of a sufficiently narrow barrier, so that the integrand is approximately constant. Choosing $ x = {x_0} $, we can then approximate Eq. (A16) as

In this way, we can see how the asymptotic properties of $ {G_{\text{tot}}}({x_0},x^{\prime}, - \text{i}\xi ) $ are related to the asymptotic properties of the Green functions $ G({x_0},x^{\prime}, - \text{i}\xi ) $ and $ G({x_0},{x_0}, - \text{i}\xi ) $. The magnitudes of all three Green functions are shown in Fig. 15 along with a direct reference calculation of $ {G_{\text{tot}}}({x_0},x^{\prime}, - \text{i}\xi ) $.

 

Figure 15. Magnitude of different electric field Green functions along the negative imaginary axis $ \omega = - \text{i}\xi $ in units of $ {G_0}( - \text{i}\xi ) = \text{c}/2\xi $. The positions $ x^{\prime} = L $ and $ {x_0} = 3L $ are both outside the original resonator and outside the region of consistency of the single dielectric barrier. Red and blue solid lines show characteristics of the Green functions for the single dielectric barrier as indicated, and the gray dashed line shows behavior of the approximate total Green function with circles indicating the corresponding reference calculation.

Download Full Size | PDF

From the above analysis, we conclude that the faster divergence of the denominator in Eq. (95) results in an overall exponential decay in the asymptotic form of $ {G_{\text{tot}}}({x_0},x^{\prime}, - \text{i}\xi ) $. Notably, the asymptotic form is independent of the actual values chosen for $ \Delta \varepsilon $ and $ \Delta X $. Therefore, we can choose them arbitrarily small so as to not influence a given QNM approximation with fixed number of terms, as argued in Subsection 4.3a.

4.3c. Region of Consistency for the Plasmonic Dimer

Following a similar approach as for the dielectric barrier in one dimension, we can now investigate the region of consistency for the dimer of gold nano-spheres from Fig. 2. To this end, we must rely on numerical calculations of the three-dimensional electric field Green tensor at complex frequencies, which is not as well-behaved as in the one-dimensional case. Therefore, we first verify that the method works by investigating the case of a single gold nano-sphere, for which we can calculate the region of consistency analytically, as shown in Appendix C.5. From the analysis, we know that the region of consistency for the single sphere is itself a sphere with the same center. Therefore, placing the center of the single gold nano-sphere at the origin, we fix the position of $ {{\bf r}^\prime} $ inside the sphere by setting $ x^{\prime} = R/2 $ and $ y^{\prime} = z^{\prime} = 0 $, and calculate the radius of the region of consistency for 100 different angles equally distributed on a circle in the plane through the center. Figure 16 shows the average relative error, as a function of the negative imaginary frequency $ \xi $, when estimating the boundary of the region of consistency as the intersection of lines corresponding to various complex frequencies $ \omega - \text{i}\xi $ along a line parallel to the imaginary axis with fixed $ \omega d/2\pi \text{c} = 0.05 $.

 

Figure 16. Relative error $ {{\cal E}_{\text{RoC}}}(\xi ) $ as a function of imaginary frequency $ \xi $ when calculating the radius of the region of consistency for a single gold sphere based on the electric field Green tensor at complex frequencies $ \omega - \text{i}\xi $ using $ \xi d/2\pi \text{c} $ and $ \xi d/2\pi \text{c} + 1 $. The black line shows a linear fit to the last three data points.

Download Full Size | PDF

 

Figure 17. Regions of consistency for the dimer of gold spheres from Fig. 2 showing, for different values of $ {{\bf r}^\prime} $ inside one of the spheres, the positions $ {\bf r} $ for which the absolute value of the electric field Green tensor tends to zero in the limit $ |z| \to \infty $. Dark gray shading corresponds to $ {{\bf r}^\prime} = {{\bf R}_1} $ in the center of the first sphere (indicated by a $ + $), and light gray corresponds to $ {{\bf r}^\prime} = {{\bf R}_1} - 0.9R{\bf x} $ close to the boundary of the first sphere (indicated by a $ \times $). Note, the dark shading extends under the light shading.

Download Full Size | PDF

At sufficiently large values of $ \xi $, the matrix problem becomes poorly conditioned, which limits the attainable relative error. This issue becomes particularly severe when using large values of the cutoff parameter $ {l_{\text{max}}} $ governing the expansion in terms of spherical wave functions, cf. Subsection 2.2a. From the analysis in Appendix C.5, however, we know that the limiting behavior of the Green tensor is governed by the components with $ l = 0 $, so we can carry out the analysis using only these wave functions, which also significantly speeds up the calculations. Despite the limited number of calculation points available, the average relative errors appear to tend to zero as a function of $ \xi $ in a polynomial manner, with the lowest attainable relative error on the order of a few parts in a thousand. Although the convergence and the accuracy appear to be markedly different from that of the exponential convergence found in the one-dimensional case, we still consider this a useful approach, since in practice we are mostly interested in the overall shape of the region of consistency.

Turning to the case of the plasmonic dimer, we start by fixing $ {{\bf r}^\prime} $ inside one of the two nano-spheres and vary the position $ {\bf r} $ in the plane through $ {{\bf r}^\prime} $ and the centers of both spheres. As in the one-dimensional case and the case of a single sphere, we find that the Green tensor behaves in an exponential manner and diverges for positions relatively far away, whereas it tends to zero at positions within a distinct curve in the plane defining the region of consistency, as shown in Fig. 17 for two different positions of $ {{\bf r}^\prime} $ inside the sphere.

As expected, the region of consistency depends on the choice of $ {{\bf r}^\prime} $, but in both cases, the region of consistency appears to be a sphere centered on the other gold nano-sphere. As in the previous cases studied, the boundary of the region of consistency appears to tend to the boundary of the sphere when $ {{\bf r}^\prime} $ is moved toward the boundary from the inside. For $ {\bf r} $ and $ {{\bf r}^\prime} $ both outside the spheres, we have found no numerically reliable examples where the Green tensor does not appear to diverge in the limit $ \xi \to \infty $.

5. Applications

In this section, we apply the QNM modeling framework to a number of problems of interest in nanophotonics. Subsection 5.1 discusses the derivation of the CMT equations using either the projection operator in Eq. (39) or the so-called field equivalence principle [139]. As an alternative to the CMT—which provides the QNM expansion of the total field—we discuss also how one can calculate the scattered field by means of the QNM expansion of the Green function inside the resonator. In Subsection 5.2, we apply the field equivalence principle to investigate QNM hybridization by calculating the QNMs of coupled systems based on the QNMs of the individual resonators. Subsection 5.3 presents the use of QNMs for perturbation theory calculations, and lastly, in Subsection 5.4, we discuss the use of QNMs for Purcell factor calculations. In all cases, we show practical applications of the results using one or both of the example material systems from Subsection 1.5.

For the scattering calculations and CMT, we note that variants of these results were derived for one-dimensional systems in Ref. [37], and for coupled cavity-waveguide systems in Ref. [75]; the CMT equations for the scattered field were previously presented in Ref [59]. For the coupled resonators in Subsection 5.2, we remark that a variant of this problem, which includes as a special case the example presented in Subsection 5.2a., was treated in Refs. [140,141]. Variants of the problem of perturbation theory with QNMs in Subsection 5.3 have been investigated in Refs. [19,21,23,27–29,51,83]. The Purcell factor derivations in Subsection 5.4 follow the approach of Ref. [17]. An identical result was presented via a different approach and for resonators made from dispersive and absorptive materials in Ref. [18], see also Refs. [60–63]. In addition, we note that variants of Purcell factor calculations with QNMs in different environments have been investigated by a number of authors and methods in Refs [56,65–67,70,85,93–95].

5.1. Scattering Calculations and CMT

As we have seen in Subsection 3.4, the QNMs are intimately related to the poles of the Green tensor. Therefore, in general, it is possible to use the QNMs to calculate the scattered field resulting from a given input field, an exercise generally referred to as the construction of the scattering matrix [45,78,79], and which is directly related to calculations of experimentally relevant quantities such as scattering and extinction cross sections [59,80]. Instead of treating the problem in a scattering framework, one can also take the point of view that the electromagnetic resonator can act as a temporal energy storage when excited by an incoming pulse. In this case, the problem is conceptually identical to the so-called (temporal) CMT, which represents a physically appealing modeling tool for integrated optical circuitry based on cavities coupled through waveguides [142–144]. Indeed, a QNM framework was used to derive the (temporal) CMT equations of coupled cavity-waveguide systems in Ref. [75], and very similar ideas were recently applied to excitation of plasmonic resonators in Ref. [59].

In classical scattering calculations, the total electromagnetic field is split in two parts corresponding to the incoming field and the scattered field as $ {\underline {\bf F} _{\text{tot}}}({\bf r},\omega ) = {\underline {\bf F} _{\text{in}}}({\bf r},\omega ) + {\underline {\bf F} _{\text{scat}}}({\bf r},\omega ) $. The incoming field is taken to be a solution to Maxwell’s equations in the background material without the electromagnetic resonator, and in general it does not obey the Silver–Müller radiation condition. The scattered field represents the change in electromagnetic field profile caused by the resonator and does obey the Silver–Müller radiation condition. Application of the projection operator in Eq. (39) requires the field $ \underline {\bf F} ({\bf r},\omega ) $ to obey the Maxwell curl equations inside the volume $ V $, as well as the Silver–Müller radiation condition. For a rather general approach to scattering calculations with QNMs, we now artificially change the total field by subtracting off the incoming field at positions outside and on the border of the volume $ V $. In essence, this approach is similar to the use of a total field/scattered field technique in numerical electromagnetism [145–147]. In this way, we define the field $ {\underline {\bf F} _{\text{TFSF}}}({\bf r},\omega ) $ as

To set up equations for the driving of the field in an electromagnetic resonator by an incoming field, we start from the defining equation for the total field,

Last, assuming that $ {\underline {\bf F} _{\text{TFSF}}}({\bf r},\omega ) $ can be expanded as in Eq. (53), we can use the same arguments as in Subsection 3.3a to arrive at a CMT type expression for the expansion coefficients of the form

To derive the temporal CMT equations, we can follow the exact same procedure as in Ref. [75] by first Fourier transforming Eq. (53) and defining

Differentiating with respect to $ t $ and rearranging the resulting integrand, we find the temporal CMT equation

5.1a. Derivation Using the Field Equivalence Principle

The same expression can be derived using the field equivalence principle and the QNM expansion of the Green tensor in Eq. (63). This approach has previously been applied to derive the CMT equations for coupled cavity-waveguide structures in Ref. [75]. The starting point is the equivalent surface currents

By inserting the relevant components of the QNM expansion for the matrix Green tensor in Eq. (62) and rearranging the terms in the form of Eq. (53), we find

From the derivation, it appears that the integration surface in Eqs. (100) or (107) must be strictly outside the material defining the resonator in order to ensure that the scattered field can be written in terms of purely outward propagating fields, or to ensure that the equivalent surface currents can be written in terms of the incoming field only. In practice, however, the integration can be performed on the inside of the boundary, since changes in the normal components of $ \tilde{\bf f}_m({{\bf r}^\prime}) $ and $ \tilde{\bf g}_m({{\bf r}^\prime}) $ across the boundary do not affect the value of the integral.

5.1b. Scattered Field Calculations

While Eq. (100) provides the QNM expansion of the total field inside the volume $ V $, it may also be interesting to calculate the expansion coefficients of the scattered field explicitly. Assuming for simplicity $ {\mu _{\text{R}}} = 1 $, we consider the scattered part of the electric field, which can be calculated based on the electric field Green function and the incoming field $ {{\bf E}_{\text{in}}}({\bf r},\omega ) $ as

In cases where the resonator is made from a homogeneous material, we can use the vector Green theorem of the second kind and the wave equation to rewrite the expression as a surface integral,

Even though Eqs. (100) and (111) are very similar in shape, they are not identical. Whereas the former provides the QNM expansion of the total field inside the volume $ V $, the latter provides only the scattered field. Consequently, the difference $ {c_n}(\omega ) = {a_n}(\omega ) - {b_n}(\omega ) $ gives the QNM expansion of the incoming field inside the volume $ V $.

5.1c. Scattering Calculations and CMT for the Dielectric Barrier

The calculations in Subsection 5.1 include a surface integral coupling the incoming field to each of the QNMs at the boundary between the total field and the scattered field regions. For the scheme to work at all positions inside the dielectric barrier, we must choose the boundary of $ V $ to coincide with the physical boundaries at $ x = \pm L/2 $. To calculate the transmitted and reflected field at positions outside the resonator, we apply the field equivalence principle in Eq. (106) with equivalent surface currents at $ x = L/2 $ and $ x = - L/2 $ given in terms of the QNM expansion of the scattered field in Eq. (109).

The QNM approximation to the total field inside the volume $ V $ can be written as

Instead of approximating the total field directly using Eq. (112), we can calculate the QNM approximation to the scattered field as

To calculate the transmitted electric field beyond the barrier at $ x \gt L/2 $, one can treat the QNM approximation to the total field at the boundary of the total field region as the input field in the field equivalence principle. The corresponding calculation of the total field at positions $ x \lt - L/2 $ is much more delicate. In principle, one can calculate the reflected field at the position $ x = - L/2 $ by subtracting off the incoming field from the total field, but the expansion in Eq. (112) fails at the position $ x = - L/2 $ because of the manifest outward propagating nature of the QNMs. One option, then, is to use the approximate field value at a position just inside the resonator where the QNM approximation converges to the correct total field. Another, much more efficient, method is to calculate the scattered field at the left boundary directly from the QNM approximation to the scattered field in Eq. (114) with $ x = - L/2 $.

 

Figure 18. Transmission through the one-dimensional dielectric barrier from Subsection 1.5a. Top: QNM approximation to the total electric field in the case of a scatting type calculation with an incoming field from the left of the form $ {E_{\text{in}}}(x,{\omega _0}) = {E_0}\exp \{ i{\omega _0}x/\text{c}\} $ with $ {\omega _0} = 4\text{c}/L $ and using 1001 QNMs with $ - 500 \le n \le 500 $. Red solid and blue dashed curves show the real and imaginary parts, respectively, and the black curve shows the absolute value. Gray shading indicates the extent of the dielectric barrier. Bottom: relative error as a function of position.

Download Full Size | PDF

Considering now the explicit case of an incoming plane wave from the left of the form $ {E_{\text{in}}}(x,{\omega _0}) = {E_0}\exp \{ i{\omega _0}x/\text{c}\} $ with $ {\omega _0} = 4\text{c}/L $, Fig. 18 shows the calculated approximation to the total electric field using $ N = 500 $ along with the relative error. At positions outside the resonator, the relative error is constant, because the approximate reflected and transmitted fields are calculated using the field equivalence principle. The failure of the QNM approximation to the total field at $ x = - L/2 $ is evident as a peak in the relative error.

To assess the convergence properties of the QNM approximation to the total field in Eq. (112), Fig. 19 shows a double logarithmic plot of the relative integrated error

 

Figure 19. Convergence analysis for the CMT approximation to the total field inside the dielectric barrier. Red circles show the relative integrated error $ {{\cal I}_N}({\omega _0}) $ in the region $ - L/4 \lt x \lt L/4 $ around the center of the barrier, and blue squares show the relative error $ {{\cal E}_N}(x,{\omega _0}) $ at $ x = 0 $ in the center of the barrier. Solid and dashed lines show linear fits to the last three data points of each set.

Download Full Size | PDF

As a curious consequence of the region of consistency, we show in Fig. 20 an example of a CMT calculation in which the total field region is defined by $ - (\pi + 1)L/2 \lt x \lt L/2 $. In this case, it follows from the expression for the region of consistency in Eq. (93) that the QNM expansion will be convergent only for $ x^{\prime} \gt 0 $. This is exactly what we observe as an abrupt jump in the field and in the associated relative errors at $ x = 0 $. At positions left of the total field region, the electric field is the sum of the incident field and the reflected field as calculated at the boundary $ x = - (\pi + 1)L/2 $, where the QNM expansion of the scattered field is not convergent. The transmitted field is calculated by use of the field equivalence principle at the right boundary $ x = L/2 $, where the QNM expansion is convergent. In comparison to Fig. 19, the relative error is larger, because the rate of convergence of the total field is not as good as when coupling the field at $ x = - L/2 $.

 

Figure 20. Top: QNM approximation to the total electric field similar to Fig. 18 but calculated for a total field region defined by $ - (\pi + 1)L/2 \lt x \lt L/2 $. Bottom: relative error as a function of position.

Download Full Size | PDF

Transmission spectrum calculations for the dielectric barrier. We return now to the transmission calculations in Fig. 1 and the claim that we can make the error in the QNM approximation to the transmission arbitrarily small. We define the transmission through the dielectric barrier as the ratio between the outgoing electric field at $ x = L/2 $ to the incoming electric field at $ x = - L/2 $. There is no scattering back at positions beyond the right-most boundary, and since the electric field is continuous at the boundary, we can write the transmission as

In Fig. 21, we show the real and imaginary parts of the transmission spectrum from Fig. 1 along with the CMT approximation as calculated using a sum similar to that in Eq. (112), but using only the three QNMs with indices $ n \in \{ 3,4,5]\} $. The bottom panel shows the associated relative error as well as the relative error in the single QNM approximation using only $ n = 4 $; the minima in the relative errors are at approximately 2% and 7%. As more QNMs are included in the sum, the minimum error decreases, and the bandwidth of the approximation increases in a symmetric way around the center point.

 

Figure 21. Top: real (red) and imaginary (blue) parts of the transmission through the dielectric barrier as a function of frequency. Dashed curves indicate the CMT approximation using the three QNMs with $ n \in \{ 3,4,5\} $. Bottom: relative error $ {\cal E}(\omega ) = |{T_{\text{CMT}}}(\omega ) - {T_{\text{ref}}}(\omega )|/{T_{\text{ref}}}(\omega ) $. Dark shading shows the error in the QNM approximation in the top panel, and light shading shows the error in the single QNM approximation using only $ n = 4 $ (not shown in the top panel).

Download Full Size | PDF

When using a modest number of QNMs, as in Fig. 21, it is useful to choose the QNMs symmetrically around the frequency of interest. Indeed, while the relative error is as low as a few percent at $ \omega L/\text{c} = 4 $, it is on the order of unity at $ \omega L/\text{c} = 2 $. The rate of convergence, however, is expected to be independent of the center point for the summation, so we can generally approximate $ {E_{\text{tot}}}(x,\omega ) $ either directly using Eq. (112) or as the sum of the incoming field and the scattered field as in Eq. (115). The associated relative errors in the transmission become simply the relative errors in the total field at $ x = L/2 $ as calculated using either Eq. (113) or Eq. (116), respectively. Figure 22 shows the relative errors in the transmission as a function of the number of QNMs in the sum. The error in the direct total field QNM approximation appears to tend to zero in a second-order polynomial fashion, which we attribute to the special point at the resonator boundary. The scattered field approach, which is inherently based on integrating across the barrier, appears to inherit the first-order polynomial convergence, as was the case also for the relative integrated error, cf. Fig. 19.

 

Figure 22. Convergence analysis for the CMT approximation to the transmission through the dielectric barrier at the frequency $ {\omega _0} = 4\text{c}/L $. Red circles and blue squares show the relative error $ {\cal E}_N^{\text{s}}(L/2,{\omega _0}) $ in calculations based on the scattered field in Eq. (115) or the relative error $ {{\cal E}_N}(L/2,{\omega _0}) $ in calculations based on the total field in Eq. (112), respectively. Solid and dashed lines show linear fits to the last three data points of each set.

Download Full Size | PDF

Supplementary code. With the supplementary code [117], we provide the files necessary to reproduce the scattering calculation in Fig. 18. We encourage interested readers to vary the number of QNMs in the sum or change the left boundary of the total field region to reproduce the curious example in Fig. 20.

5.1d. Coupled Mode Theory for the Plasmonic Dimer

We now turn to the plasmonic dimer and consider the electromagnetic response of the dimer when illuminated by a plane wave of magnitude $ {E_0} $ polarized along the dimer axis. Figure 23 shows the broadband electric field at the position $ {{\bf r}_0} $ directly in the middle between the two spheres. As in the case of the Purcell factor in Fig. 2, a number of distinct peaks are visible, each of which can be attributed primarily to a single QNM. The maximum of the highest peak, for example, is found at $ \omega d/2\pi \text{c} \approx 0.11 $, which corresponds roughly to the real part of $ {\tilde \omega _1} $, cf. Eq. (11).

 

Figure 23. Magnitude of the local electric field at the point $ {{\bf r}_0} $ in the middle between the two gold spheres of the plasmonic dimer when illuminated by a plane wave of magnitude $ {E_0} $ and polarization along the dimer axis. Red dashed and solid curves show the single- and two-QNM CMT approximations when using $ {\tilde {\underline {\bf F} }_1}({\bf r}) $ or $ {\tilde {\underline {\bf F} }_1}({\bf r}) $ and $ \tilde {\underline {\bf F} }_1^*({\bf r}) $, respectively, and gray shading shows the reference calculation.

Download Full Size | PDF

The QNM $ {\tilde {\underline {\bf F} }_1}({\bf r}) $ is commonly referred to as the “dipolar mode” of the dimer, because of the field pattern, which resembles the field between two particles with opposite charges. From a physical point of view, this mode is therefore expected to couple effectively to the incoming plane wave, wherefore it is also known as a “bright mode.” Other modes are so-called “dark modes” [148], because they cannot be effectively excited from the far-field. Typically, the distinction between bright and dark modes is inferred from symmetry considerations of the exciting field and the QNMs as in Ref. [149], for example. The expression for the expansion coefficient in Eq. (100) exactly captures the excitation of the different QNMs by an incoming field, and hence provides a precise mathematical distinction between bright and dark modes in general resonators.

Even if the QNM expansion is not formally convergent at the position of interest between the two spheres, we can extend the region of consistency by the procedure discussed in Subsection 4.3a. In this case, for example, we can imagine including a sphere of finite but vanishingly small permittivity difference between the gold spheres making up the plasmonic dimer. Alternatively, we can embed the entire dimer in a large sphere of vanishingly small permittivity difference with respect to the background. In both cases, we expect then the QNM expansion to be convergent at the position $ {{\bf r}_0} $. At the same time, the change in the QNMs of interest can be made arbitrarily small (see also Subection 5.3d), so the expansion coefficients $ {a_n}(\omega ) $ of any truncated series can be calculated simply by use of the QNMs of the original plasmonic dimer. The red dashed and solid curves in Fig. 23 show the magnitude of the CMT approximation to the electric field when using either $ {\tilde {\underline {\bf F} }_1}({\bf r}) $ or $ {\tilde {\underline {\bf F} }_1}({\bf r}) $ and $ \tilde {\underline {\bf F} }_1^*({\bf r}) $. As in the case of the dielectric barrier, the CMT approximation captures not only the magnitude, but also the phase of the electromagnetic field. Focusing on the low-frequency region, Fig. 24 shows the real and imaginary parts of the electric field reference calculations, as well as the relative error. The relative error of the two-QNM CMT approximation is below 0.05 over a large bandwidth and as low as a few parts in a thousand at best.

 

Figure 24. Top: real (red) and imaginary (blue) parts of the electric field at the point $ {{\bf r}_0} $ in the center between the two spheres in the plasmonic dimer, when illuminated by a plane wave of magnitude $ {E_0} $ with polarization along the dimer axis. Dashed curves indicate the CMT approximation using the QNMs $ {\tilde {\underline {\bf F} }_1}({\bf r}) $ and $ \tilde {\underline {\bf F} }_1^*({\bf r}) $. Bottom: relative error $ {\cal E}({\bf r},\omega ) = |{E_{\text{ref}}}({\bf r},\omega ) - {E_{\text{CMT}}}({\bf r},\omega )|/|{E_{\text{ref}}}({\bf r},\omega )| $ as a function of frequency.

Download Full Size | PDF

5.2. Coupled Resonators

When two electromagnetic resonators are placed in close proximity, we expect them to couple to each other and produce new, hybridized resonances. Indeed, this is a fundamental phenomenon, which can be observed in all areas of physics. In the so-called tight-binding approach, or the linear combination of atomic orbitals (LCAO), the hybridized resonances of coupled, localized wave functions can be calculated through the use of an overlap integral [150]. In this way, one can immediately appreciate the origin of the coupling, and one can use this methodology to enhance or suppress the coupling, or to calculate approximate band structures, for example. The QNMs, however, are not localized, so we cannot immediately apply the LCAO framework. The exponential divergence, in particular, would lead to larger and larger overlap integrals for increasing distance between the resonators. For QNMs in layered systems, this problem was addressed in Refs. [140,141] by setting up and solving a self-consistent matrix equation built from the QNMs of individual layers in order to find the poles of the scattering matrix. Below, we take a slightly different approach and use the QNMs of the individual resonators in combination with the field equivalence principle to set up a set of equations for the expansion coefficients of the scattered field in the compound system.

On the boundary of a given resonator, we assume that the scattered field at frequency $ \omega $ can be expanded as in Eq. (109). The scattered field is radiated away and may itself lead to scattering off another resonator. In this way, and quite analogous to multiple scattering theory, we can set up a self-consistent set of equations for the scattered field expansion coefficients $ b_m^{(i)} $, where $ i $ now counts the scatterers. By use of Eq. (106) with equivalent surface currents given in terms of the QNM expansion in Eq. (109), the scattered field from resonator $ j $ at position $ {\bf r} $ is given as

5.2a. Coupled Dielectric Barriers

We now consider the system made from two copies of the dielectric barrier in Subsection 1.5a that are positioned a distance of $ D = L $ apart. This system does not allow a simple, closed-form expression for the QNMs or the resonance frequencies, and so we turn to numerical methods. In practice, we use a one-dimensional formulation of the VIE formulation from Subsection 2.2a, but we note that these calculations can be performed with practically all the methods discussed in Section 2. Scattering between the two barriers leads to a rather complicated spectrum of complex QNM resonance frequencies, which in general cannot be understood from a simple hybridization picture. For the particular separation of choice, however, we find two QNM resonance frequencies close to $ {\tilde \omega _1} $ of the single barrier system,

 

Figure 25. Top: complex spectrum of the double barrier system showing the two QNM frequencies near the position of $ {\tilde \omega _1} $ of the single barrier system (indicated by a black circle). Center and bottom: mode profiles showing the absolute values (black solid lines) as well as the real (red solid) and imaginary (blue dashed) parts of the electric field QNMs corresponding to $ {\tilde \omega _{\text{even}}} $ (center) and $ {\tilde \omega _{\text{odd}}} $ (bottom). Gray shading indicates the extent of the dielectric barriers.

Download Full Size | PDF

To set up Eq. (120) for the double barrier system, we note that the boundaries of the resonators are simply the four points $ x_1^ - $, $ x_1^ + $, $ x_2^ - $, and $ x_2^ + $, indicating the left ($ - $) and right ($ + $) boundaries of resonators 1 and 2. Starting from Eq. (111), we can use the wave equations to formulate $ b_m^{(1)} $ and $ b_m^{(2)} $ in terms of the electric field only as

Expanding the scattered fields as in Eq. (109) symmetrically around $ n = 1 $,

 

Figure 26. Convergence analysis for the calculation of $ {\tilde \omega _{\text{even}}} $ (red circles) $ {\tilde \omega _{\text{odd}}} $ (blue squares) of the double barrier system by means of a QNM expansion showing the relative error $ {{\cal E}_{\text{QNM}}} = |{\tilde \omega _{\text{QNM}}} - {\tilde \omega _{\text{ref}}}|/|{\tilde \omega _{\text{ref}}}| $ as a function of number of QNMs used in the expansions (in each barrier). Solid and dashed lines show the corresponding fits to the last four data points.

Download Full Size | PDF

5.3. First-Order Perturbation Theory

If the material defining the electromagnetic resonator is slightly perturbed by a local change $ \Delta \epsilon ({\bf r}) $ or $ \Delta\mu({\bf r}) $, we expect to be able to calculate the dominant change in the QNM resonance frequency using first-order perturbation theory [19]. Indeed, we may in principle calculate the change in the QNM resonance frequency—as well as changes in the QNM field profile—to arbitrary accuracy by going to perturbation theory of sufficiently high order. In practice, the use of higher-order perturbation theory quickly becomes cumbersome, especially in the so-called Rayleigh–Schrödinger formulation [151]. An alternative to the Rayleigh–Schrödinger formulation, which is known as the Brillouin–Wigner formulation, expresses the perturbed quantity as an implicit series expansion [151]. As yet another alternative, it was illustrated in Ref. [31] that one can use the QNMs of the unperturbed system to formulate a linear eigenvalue problem for the perturbed system, a method referred to as resonant state expansion [32].

5.3a. Perturbations in the Permittivity

To see how the normalization integral in Eq. (42) arises naturally in this process, we write the defining equation for the $ m^{\prime}$th QNM of the perturbed system as

Inserting in Eq. (129) and expanding to first order, we find that the perturbation $ \Delta \tilde {\underline {\bf F} }({\bf r},\omega ^{\prime}) $ solves the equation

Multiplying from the left with $ \tilde {\underline {\bf F} }_m^{\ddagger }({\bf r})\underline {\underline {\bf W} } $ and integrating, we can use Eq. (30) to find the expression

Last, writing out the expression for the surface integral and making use of Eqs. (41) and (43), we can write the first-order correction in the form

In cases with $ \Delta \mu({\bf r}) = 0 $, this reduces to the well-known form

5.3b. Shifting Boundaries

Perturbation theory with shifting boundaries represents a nontrivial extension of the theory, since the local change in material properties close to the boundary of the resonator material will generally change in a nonperturbative way. The problem has been solved by Lai et al. [19] and by Johnson et al. [152] for Hermitian eigenvalue problems, but for which the central arguments immediately carry over to the case of first-order perturbation theory of QNMs. For the case of a resonator made from a nonmagnetic material in a background with permittivity $ { \epsilon _{\text{B}}} $, we can, therefore, immediately infer that [19,152]

5.3c. Perturbation Theory for the Dielectric Barrier

To illustrate the practical use of Eq. (136), we use the one-dimensional example of the dielectric barrier, for which the QNM resonance frequencies are known analytically, cf. Eq. (5). Assuming the permittivity of the dielectric barrier is changed by a small and constant amount $ \Delta \epsilon $, we find from Eq. (136) that the first-order change in resonance frequency should be

To assess the validity of this result, we can compare it directly with the explicit expression in Eq. (5), which is valid for all values of $ {n_{\text{R}}} $. Substituting $ {n_{\text{R}}} + \Delta n $ in place of $ {n_{\text{R}}} $, we can make a series expansion around the point $ \Delta n = 0 $ to find that the first-order change in the resonance frequency is indeed given by Eq. (140). Figure 27 shows the real and imaginary parts of $ {\tilde \omega _4} $ when the permittivity is changed by adding or subtracting as much as $ \Delta \epsilon = \pm 5 $. The first-order approximation, as calculated from Eq. (140), provides the tangents to the curves for the exact results in the point $ \Delta \epsilon = 0 $, as expected.

 

Figure 27. Top: real (red circles) and imaginary (blue squares) parts of the shift in the complex QNM resonance frequency $ {\tilde \omega _4} $ as a function of permittivity change $ \Delta \epsilon $ of the dielectric barrier. Solid and dashed lines show the real and imaginary parts of the corresponding first-order perturbation theory result. Bottom: absolute value of the difference between the exact result and the first-order perturbation theory result.

Download Full Size | PDF

5.3d. Perturbation Theory for the Plasmonic Dimer

To illustrate the influence of the vector nature of the QNMs, we consider the case where a small sphere with radius $ {R_\Delta } $ and local constant permittivity change $ \Delta \epsilon $ are inserted at the point $ {{\bf r}_0} $ directly between the gold spheres of the plasmonic dimer.

For simplicity, we approximate the integral by assuming the QNMs are constant throughout the volume $ {V_\Delta } $ of the small sphere, in which case Eq. (136) can be written as

The assumption that the electric field QNM $ \tilde{\bf f}_1({\bf r}) $ is constant throughout the volume $ \Delta V $ is a remarkably good approximation, even for relatively large spheres, as measured by the difference

 

Figure 28. Top: real (red circles) and imaginary (blue squares) parts of the shift in the complex QNM resonance frequency $ {\tilde \omega _1} $ as a function of permittivity change $ \Delta \epsilon $ of a sphere with radius $ {R_\Delta } = R/4 $ in the center between the gold spheres of the plasmonic dimer. Dashed lines indicate the first-order perturbation theory result, and solid lines show the result of more properly accounting for the vector nature of the QNM by use of the known polarizability of a small sphere. Bottom: relative difference $ {\cal E}(\Delta \epsilon ) = |\Delta \tilde \omega - \Delta {\tilde \omega _1}|/|\Delta \tilde \omega | $ between the solid lines and the reference calculations in the top panel.

Download Full Size | PDF

The factor $ \Delta \epsilon /(1 + \Delta \epsilon /3) $ evidently plays the role of a linear susceptibility, so we can immediately substitute it in place of $ \Delta \epsilon $ in Eq. (141). In this way, we obtain an approximation to first order in the constant field, $ \tilde{\bf f}_1({\bf r}) \approx \tilde{\bf f}_1({{\bf r}_0}) $, which treats the change in polarization to all orders. The solid black curves in Fig. 28 show the results of this approach, which clearly captures the correct resonance frequency change over a much larger range of permittivity changes than the direct application of Eq. (141).

Shifting boundaries. To illustrate the practical use of Eq. (137), we consider the shift in complex resonance frequency $ {\tilde \omega _1} $ from the example in Fig. 2 resulting from increasing or decreasing the size of the gold spheres. At each point on the surface, the (normalized) electric field QNM $ \tilde{{\bf f}}_1({\bf r}) $ is split in a parallel and a normal component to enable the evaluation of the integral in Eq. (137). Figure 29 shows the real part of the differential contribution to the frequency shift plotted on the surface of the dimer. Compared to the field profile in Fig. 2, the shape of the QNM is partly recognized by the relatively large contribution from the parts of the spheres closest to the dimer center. By standard numerical integration, we can now evaluate the surface integrals to find

 

Figure 29. Real part of the differential contribution to the complex frequency shift of $ \Delta {\tilde \omega _1}/\Delta R $ in units of $ 2\pi \text{c}/{d^2} $ plotted on the surface of the dimer.

Download Full Size | PDF

As a reference, we compare this to high-accuracy calculations of the resonance frequency by use of the VIE method as detailed in Subsection 2.2a. Figure 30 shows the movement of $ {\tilde \omega _1} $ in the complex plane as a function of varying radius of the spheres. The first-order perturbation result in Eq. (147) provides the tangent to the curve in the point $ z = {\tilde \omega _1} $ as expected.

 

Figure 30. Zoom of the Complex QNM spectrum in Fig. 2 in the vicinity of $ {\tilde \omega _1} $. The full curve shows the movement of the complex resonance frequency when changing the radius of the gold spheres continuously between 40 nm and 60 nm as indicated. The dashed curve shows the tangent to the curve as given by first-order perturbation theory with shifting boundaries.

Download Full Size | PDF

5.4. Purcell Factor Calculations

The Purcell factor [16] is a common figure of merit for nanophotonic resonators and provides a measure for the enhanced emission rate of an electric dipole source inside the resonator relative to the emission in a homogeneous material. In the so-called weak coupling regime, it is well known that the rate of spontaneous emission of a quantum emitter, with resonance frequency $ \Omega $ and dipole moment $ {\bf d} = d {\hat{\bf e}_{\text{d}}} $, is proportional to the imaginary part of the transverse electric field Green tensor [153,154],

Upon dividing the general result in Eq. (148) by the corresponding emission in a homogeneous medium of refractive index $ {n_{\text{R}}} $, for which we have the relation

If the quantum emitter is placed in a resonator, we can expand the transverse Green tensor by use of the QNMs. For instance, using the first of the expressions in Eq. (90), we find

When a single QNM $ \tilde{\bf f}_m({{\bf r}_0}) $ is sufficient to adequately approximate the Green tensor, and when the emitter frequency is tuned to coincide with the real part of the QNM resonance frequency, $ \Omega = {\omega _m} $, we find the approximate relation

Last, using $ Q = {\omega _m}/2{\gamma _m} $, as defined in Eq. (26), and $ { \epsilon _{\text{R}}} = n_{\text{R}}^2 $, we can write the expression in the exact form due to Purcell [16],

If the emitter is not tuned into resonance with the QNM, one can write the approximate Purcell factor by multiplying a correction factor onto the expression in Eq. (153), but it is typically much easier to work directly with the QNM Green tensor expansion in Eq. (151), as we do below.

5.4a. Purcell Factor in Gap Center of the Plasmonic Dimer

The plasmonic dimer supports a large number of QNMs, some of which show up as distinct peaks in scattering spectra or Purcell factor calculations, as noted already in a somewhat qualitative manner in connection with Figs. 2 and 3. In this section, we discuss the practical application of the QNM approximations in detail. The red dashed curve in Fig. 2 shows the result of approximating the electric field Green tensor by a single term as

The approximate expression arguably captures the main qualitative features of the spectrum around the resonance, and the relative error at the peak is as low as a few percent. The agreement can be dramatically improved, however, by using the symmetry property of the QNMs to ensure the correct crossing behavior of the Green tensor. By adding to Eq. (155) the corresponding term from the mode $ \tilde{{\bf f}}_1^*({\bf r}) $ at the complex resonance frequency $ - \tilde \omega _1^* $, we find

Large bandwidth approximation. To obtain a deeper understanding of the convergence properties of the QNM Green tensor expansion, it is illustrative to calculate it over a wider frequency range that covers several resonances. The high-symmetry point between the two nano-spheres is well suited for this purpose, because a large fraction of the QNMs vanishes at this point. Tables 1 and 2 list the complex resonance frequencies as well as the generalized effective mode volumes of the first five QNMs of interest with nonvanishing electric field components along the dimer axis at the point in the center of the gap between the two spheres. These QNMs are located in the region $ 0 \lt {\omega _n} \lt {\omega _{\text{sp}}} $ and $ {\gamma _1} \le {\gamma _n} $, where $ {\omega _{\text{sp}}} = {\omega _{\text{p}}}/\sqrt 2 $ is the frequency at which the in-plane wave vector $ {k_\parallel }(\omega ) $ of a surface plasmon polariton on a planar Drude metal surface in air tends to infinity [154]. In this limit, the associated in-plane wavelength $ {\lambda _\parallel }(\omega ) = 2\pi /{k_\parallel }(\omega ) $, which describes propagation along the interface, tends to zero, wherefore the finite curvature of the spheres becomes irrelevant and the spheres look locally like plane surfaces. Nevertheless, the spherical shapes lead to a periodicity requirement on the QNMs, which can only be met at certain in-plane wavelengths. As the spacing between the in-plane wavelengths becomes smaller, one can fit an ever increasing number of oscillations around the sphere, and this effect shows up as an accumulation point in the spectrum at $ \omega \approx {\omega _{\text{sp}}} $. Figure 31 shows the Purcell factor from Fig. 2 along with Purcell factor calculations based on Green tensor approximations of the form

Table 1. Complex Resonance Frequencies

View Table | View all tables in this article

Table 2. Generalized Effective Mode Volumes

View Table | View all tables in this article

 

Figure 31. Top: Purcell factor, Eq. (150), for a dipole oriented along the dimer axis at the position $ {{\bf r}_0} $ in the center between the two gold spheres as in Fig. 2 but showing a larger bandwidth. Gray shading shows the reference calculation $ F_{\text{P}}^{\text{ref}} $. The red dashed curve shows the approximate Purcell factor resulting from Eq. (155), and the solid colored lines show the approximations $ {F_{\text{P},N}} $ resulting from Eq. (150) by substituting the QNM approximations to the Green tensor in Eq. (158) for increasing values of $ N $. Red, green, blue, purple, and dark gray curves show the result $ {F_{\text{P},N}} $ with $ N = 1,2,3,4,5 $, respectively. The dashed vertical line indicates the position of $ {\omega _{\text{sp}}} $. Bottom, relative error $ |\Delta {F_{\text{P}}}(\omega )|/|F_{\text{P}}^{\text{ref}}(\omega )| $ corresponding to the calculations in the top panel.

Download Full Size | PDF

6. Conclusion

In this Tutorial, we have presented a biorthogonal approach to modeling electromagnetic resonators using QNMs and used this approach, in combination with examples of a dielectric barrier and a plasmonic dimer in one and three dimensions, respectively, to analyze and illustrate the usefulness and limitations of resonator models based on QNMs.

Starting from the QNMs as solutions to the Maxwell curl equation with the Silver–Müller radiation condition, we have illustrated how a biorthogonal approach leads to a useful definition of adjoint QNMs and an associated projection operator. When the operator is applied to project a QNM onto itself, it provides a well-known formulation of the QNM norm. By applying the operator to certain classes of electromagnetic field problems, including Green tensor calculations, one can derive formal expansions of the fields in terms of QNMs. For the case of the Green tensor, we have discussed how this approach compares to an alternative approach in which the QNMs are directly associated with the residues of the analytical continuation of the Green tensor.

The question of if and when the formal solutions are convergent and consistent was discussed in some detail, building on ideas from the literature. Much of the convergence analysis is based on the central enabling fact that QNM expansions of the electromagnetic Green tensor $ \underline {\underline {\bf G} } ({\bf r},{{\bf r}^\prime},\omega ) $ are possible for positions within the resonators and possibly beyond, depending on the parameters of the resonator and position $ {{\bf r}^\prime} $, as discussed in Subsection 4.3. In general, for a given choice of $ {{\bf r}^\prime} $, one can numerically trace out a region of consistency for which the Green tensor expansion is expected to converge to the correct result in the limit when infinitely many QNMs are included in the expansion. Even then, in certain cases, there might be additional contributions stemming from branch cuts in the analytical continuation of the Green tensor. At positions inside the region of consistency, the actual convergence is generally a relatively slow function of the number of QNMs $ N $; for the dielectric barrier, where we can investigate the convergence explicitly, we find it to be of the order $ 1/N $ in general, cf. Figs. 9 and 26. For certain special points in the CMT calculations, such as the center of the barrier and the right-most boundary, the QNM approximation was found to sometimes converge as $ 1/{N^2} $ as seen in Figs. 19 and 22.

Depending on the calculation method and level of accuracy, the computational costs for individual QNMs are relatively large. In combination with the modest convergence properties, this means that the appeal of a QNM decomposition is largest when the response of the system is governed mainly by a few QNMs. In such cases, a QNM approach provides a physically appealing and analytically tractable mode decomposition that is often very accurate. In the transmission spectrum for the dielectric barrier in Fig. 1, the relative error of the one and three QNM approximations is as low as 7% and 2%, respectively, as analyzed in Subsection 5.1c.. Similarly, in the Purcell factor approximation in Fig. 2, the relative errors of the one and two QNM approximations are as low as a few percent and one part in a thousand, as discussed in Subsection 5.4a.

The final part of the Tutorial was devoted to four examples of QNM modeling: scattering calculations and CMT, coupled resonators, perturbation theory, and Purcell factor calculations. The derivation of the scattering calculations in Subsection 5.1 provided an interesting display of the importance of the radiation conditions. Indeed, the total field in a scattering calculation does not fulfill the radiation condition, and therefore cannot be expanded in terms of the QNMs by use of the projection operator derived in Subsection 3.2. The scattered field, however, does fulfill the radiation condition, and by splitting the total field in an incoming and a scattered field part, the projection onto the QNMs leads to a physically appealing expression in which the incoming field acts as a source for the total field inside the resonator. To derive expressions for the scattered field, one can instead work with an integral equation and the well-known QNM expansion of the Green tensor. The same integral equation can be used to set up a linear equation system for calculating the hybridized QNMs in coupled resonances, as illustrated in Subsection 5.2. Conceptually, such an approach is well known from tight-binding type theories, but the exponential divergence of the QNMs outside the individual resonators complicates the formulation in practice. Subsection 5.3 presented a number of applications of first-order perturbation theory in order to illustrate how the QNM normalization integral arises naturally in this process, and to show practical limitations and extensions, such as the case of shifting boundaries and the use of the polarizability of small scatterers to improve the accuracy. Lastly, in Subsection 5.4, we discussed the use of QNMs in the calculation of Purcell factors, and showed how the original formula due to Purcell arises naturally from a single QNM approximation to the Green tensor.

We hope that this Tutorial can serve to illustrate how the use of QNMs provides a mathematically rigorous framework for the modeling of a multitude of different phenomena associated with electromagnetic resonators. As argued in the introduction, many models of resonator phenomena—such as laser models or models for light propagation through coupled cavity-waveguide structures—have been implicitly relying on a mode decomposition, although almost exclusively by treating the cavity as a closed system and rarely by explicitly defining the cavity modes as QNMs obeying a radiation condition. By fully exploiting a QNM approach, one can usually obtain models of the same complexity, but with explicit and precise definitions of the various model parameters, such as mode volumes and coupling constants. In turn, such QNM-based models facilitate both a qualitative understanding of the underlying physics and a high level of quantitative accuracy in calculations of observables for the direct interpretation of experiments. In combination, these characteristics offer unique predictive capabilities, which may enable the generation of optimized designs with considerably reduced effort relative to optimization schemes based on direct numerical solution of Maxwell’s equations.

Appendix A: Green Functions and Green Tensors

The Green tensor plays a prominent role throughout this Tutorial and in many areas of theoretical physics in general. In this appendix, therefore, we provide an introduction to the concept of Green functions and Green tensors with an emphasis on their use in electromagnetics. Starting from the scalar wave equation, we discuss how one can construct an equation for the associated Green function and why this function is particularly useful for the solution of the inhomogeneous problem. The electromagnetic field is a three-dimensional vector field in general, however, and therefore one must generalize the concept of the Green function to a Green tensor. Although we cannot cover the entire subject, we hope to provide a sufficient foundation to at least make this Tutorial self-contained. For a general and more thorough exposition of the subject, we refer to general textbooks on applied mathematics, such as those by Morse and Feshbach [163] or Arfken and Weber [164]. For a classical and very thorough presentation focusing on electromagnetics, see the book by Tai [165].

A.1. Electric Field Green Function in One Dimension

As a starting point, we consider the scalar wave equation in one dimension for the electric field in an environment described by the permittivity $ { \epsilon _{\text{R}}}(x,\omega ) $, subject to a localized and spatially varying current density $ j(x,\omega ) $,

(A1)$$\partial _{xx}^2E(x,\omega ) + { \epsilon _{\text{R}}}(x,\omega )\frac{{{\omega ^2}}}{{{\text{c}^2}}}E(x,\omega ) = - \text{i}\omega {\mu _0}j(x,\omega ),$$

and with the additional requirement that the electromagnetic field obeys the Silver–Müller radiation condition at $ x \to \pm \infty $, where the permittivity is assumed to be constant. The usefulness of the Green function stems from the fact that the solution to Eq. (A1) can be written as

(A2)$$E(x,\omega ) = \text{i}\omega {\mu _0}\int_{ - \infty }^\infty G(x,x^{\prime},\omega )j(x^{\prime},\omega )\text{d}x^{\prime},$$

where the Green function $ G(x,x^{\prime},\omega ) $ is the solution to the differential equation

(A3)$$\partial _{xx}^2G(x,x^{\prime},\omega ) + { \epsilon _{\text{R}}}(x,\omega )\frac{{{\omega ^2}}}{{{\text{c}^2}}}G(x,x^{\prime},\omega ) = - \delta (x - x^{\prime}),$$

subject to the Silver–Müller radiation condition. Indeed, inserting Eq. (A2) in Eq. (A1), dividing by $ \text{i}\omega {\mu _0} $, and moving the differentiation with respect to $x$ inside the integral over $x^{\prime}$, we find the relation

(A4)$$\int_{ - \infty }^\infty \left[ {\partial _{xx}^2G(x,x^{\prime},\omega ) + { \epsilon _{\text{R}}}(x,\omega )\frac{{{\omega ^2}}}{{{\text{c}^2}}}G(x,x^{\prime},\omega )} \right]j(x^{\prime},\omega )\text{d}x^{\prime} = - j(x,\omega ),$$

and by use of Eq. (A3) and the defining property of the $ \delta $-function,

(A5)$$\int_{ - \infty }^\infty \delta (x - x^{\prime})f(x^{\prime})\text{d}x^{\prime} = f(x),$$

we find that $ E(x,\omega ) $, as given by Eq. (A2), is indeed a solution to the original problem. From the form of the solution, the Green function can be seen as a propagator from the source position at $ x^{\prime} $ to the observation point $ x $. Due to the linearity of the differential equation, the solution can be calculated as the superposition of contributions stemming from the source at different points $ x^{\prime} $. The solution obeys the correct radiation condition because the Green function is required to obey this condition as well.

A.1a. Background Green Function in One Dimension

For the special case of a homogeneous permittivity $ { \epsilon _{\text{R}}}(x,\omega ) = { \epsilon _{\text{B}}} = n_{\text{B}}^2 $, there are a number of different ways for calculating the explicit solution to Eq. (A3), which we refer to as the “background Green function” and denote by $ {G_{\text{B}}}(x,x^{\prime},\omega ) $. One relatively general solution method, which we will use also to calculate the Green function for the dielectric barrier in Appendix C.1b, consists of expanding the unknown Green function in basis functions obeying the homogeneous differential equation for $ x \ne x^{\prime} $ as well as the radiation condition. In the case of a homogeneous permittivity, we shall therefore use the ansatz,

(A6)$${G_{\text{B}}}(x,x^{\prime},\omega ) = \left\{{\begin{array}{*{20}{c}}{A(x^{\prime},\omega ){\text{e}^{ - \text{i}{k_{\text{B}}}x}}}&{x \lt x^{\prime}}\\{B(x^{\prime},\omega ){\text{e}^{\text{i}{k_{\text{B}}}x}}}&{x \gt x^{\prime}}\end{array}} \right.,$$

where $ {k_{\text{B}}} = {n_{\text{B}}}\omega /\text{c} $. Demanding continuity at $ x = x^{\prime} $, we find

(A7)$$A(x^{\prime},\omega ) = B(x^{\prime},\omega ){\text{e}^{\text{i}2{k_{\text{B}}}x^{\prime}}}.$$

To find an expression for $ B(x^{\prime},\omega ) $, we integrate Eq. (A3) from $ x = x^{\prime} - \epsilon $ to $ x = x^{\prime} + \epsilon $ to find the condition

(A8)$$\mathop {\lim }\limits_{ \epsilon \to 0} \{ {G_{\text{B}}}^{\prime}(x^{\prime} + \epsilon ,x^{\prime},\omega ) - {G_{\text{B}}}^{\prime}(x^{\prime} - \epsilon ,x^{\prime},\omega )\} = - 1,$$

where $ {G_{\text{B}}}^{\prime}(x,x^{\prime},\omega ) $ denotes the derivative of $ {G_{\text{B}}}(x,x^{\prime},\omega ) $ with respect to $ x $. Inserting Eq. (A6), and making use of the relation in Eq. (A7), we find the condition

(A9)$$B(x^{\prime},\omega ) = \frac{\text{i}}{{2{k_{\text{B}}}}}{\text{e}^{ - \text{i}{k_{\text{B}}}x^{\prime}}}.$$

Putting it all together, we find that the background Green function in one dimension can be written in terms of the distance $ X = |x - x^{\prime}| $ as

(A10)$${G_{\text{B}}}(X,\omega ) = \frac{\text{i}}{{2{k_{\text{B}}}}}{\text{e}^{\text{i}{k_{\text{B}}}X}}.$$

A.1b. Scattering Calculations Based on the Green Function

In typical scattering calculations, we consider an incoming field $ {E_{\text{in}}}(x,\omega ) $, which is a solution to Eq. (A1) in the homogeneous background defined by constant permittivity $ { \epsilon _{\text{R}}}(x,\omega ) = { \epsilon _{\text{B}}} $,

(A11)$$\partial _{xx}^2{E_{\text{in}}}(x,\omega ) + { \epsilon _{\text{B}}}{k^2}{E_{\text{in}}}(x,\omega ) = - \text{i}\omega {\mu _0}j(x,\omega ),$$

and we seek to calculate the effects of a scattering object defined by the permittivity difference $ \Delta \epsilon (x,\omega ) = { \epsilon _{\text{R}}}(x,\omega ) - { \epsilon _{\text{B}}} $. In many cases, the source $ j(x,\omega ) $ is taken to be far away from the scattering object, although this is not a requirement. As a special case, the use of a plane wave as the incoming field corresponds to moving the source infinitely far away. To proceed, we define the total field as a sum of the incoming field and a scattered field, $ {E_{\text{tot}}}(x,\omega ) = {E_{\text{in}}}(x,\omega ) + {E_{\text{scat}}}(x,\omega ) $, and insert this definition in the wave equation. Making use of the relation $ { \epsilon _{\text{R}}}(x,\omega ) = \Delta \epsilon (x,\omega ) + { \epsilon _{\text{B}}} $ and collecting terms corresponding to the incoming field and the scattered field, respectively, we find that the scattered field obeys the equation

(A12)$$\partial _{xx}^2{E_{\text{scat}}}(x,\omega ) + { \epsilon _{\text{R}}}(x,\omega ){k^2}{E_{\text{scat}}}(x,x^{\prime},\omega ) = - {k^2}\Delta { \epsilon _{\text{R}}}(x,\omega ){E_{\text{in}}}(x,\omega ).$$

In comparison to Eq. (A1), the term with the incoming field evidently takes the place of the source current term on the right-hand side of the equation. Therefore, if we know the Green function for the full problem $ { \epsilon _{\text{R}}}(x,\omega ) = \Delta \epsilon (x,\omega ) + { \epsilon _{\text{B}}} $, we can use it to calculate the scattered field by use of Eq. (A2) with $ {k^2}\Delta { \epsilon _{\text{R}}}(x,\omega ){E_{\text{in}}}(x,\omega ) $ in place of $ \text{i}\omega {\mu _0}j(x,\omega ) $. By adding the incoming field, we can express the total field as

(A13)$${E_{\text{tot}}}(x,\omega ) = {E_{\text{in}}}(x,\omega ) + {k^2}\int_{ - \infty }^\infty G(x,x^{\prime},\omega )\Delta \epsilon (x^{\prime},\omega ){E_{\text{in}}}(x^{\prime},\omega )\text{d}x^{\prime}.$$

The use of Eq. (A13) relies on the Green function for the full problem, which may be nontrivial to calculate, especially in higher dimensions. In such cases, one can often benefit from an approximation to the Green function based on a limited number of QNMs.

In cases where one does not have access to the Green function for the full problem, the total field can be formulated in an alternative way by use of the background Green function as

(A14)$${E_{\text{tot}}}(x,\omega ) = {E_{\text{in}}}(x,\omega ) + {k^2}\int_{ - \infty }^\infty {G_{\text{B}}}(x,x^{\prime},\omega )\Delta \epsilon (x^{\prime},\omega ){E_{\text{tot}}}(x^{\prime},\omega )\text{d}x^{\prime},$$

where $ {E_{\text{tot}}}(x,\omega ) $ is now the solution of an implicit integral equation, since it appears on both sides of the equality. Such types of equations are known in the mathematical literature as Fredholm equations, and in the physics literature on scattering phenomena, Eq. (A14) is often referred to as the Lippmann–Schwinger equation [166]. The equivalence between Eqs. (A14) and (A1) can be seen by first writing the latter as

(A15)$$\partial _{xx}^2{E_{\text{tot}}}(x,\omega ) + { \epsilon _{\text{B}}}{k^2}{E_{\text{tot}}}(x,\omega ) = - \text{i}\omega {\mu _0}j(x,\omega ) - \Delta \epsilon (k,\omega ){E_{\text{tot}}}(x,\omega ).$$

Next, inserting the expression in Eq. (A14) for the total field and making use of Eq. (A11), we can follow identical steps as those taken in connection with Eq. (A4) to see that the equality is indeed fulfilled.

Comparing Eqs. (A1) and (A3), it follows that the Green function can be interpreted as the electric field at the position $ x $ due to a point source at $ x^{\prime} $, and this provides a powerful way of calculating the Green function in general structures. By substituting $ {G_{\text{B}}}(x,x^{\prime},\omega ) $ for $ {E_{\text{in}}}(x,\omega ) $ in Eq. (A14) and interpreting $ {E_{\text{tot}}}(x,\omega ) $ as a Green function, we find the relation

(A16)$$G(x,x^{\prime},\omega ) = {G_{\text{B}}}(x,x^{\prime},\omega ) + \frac{{{\omega ^2}}}{{{\text{c}^2}}}\int {G_{\text{B}}}(x,y,\omega )\Delta \epsilon (y,\omega )G(y,x^{\prime},\omega )\text{d}y.$$

In the physics literature on scattering phenomena, this is known as the Dyson equation [167].

A.2. Electric Field Green Tensor in Three Dimensions

The basic principle of the Green function is very general, and the Green function approach can be applied also to vector fields in higher dimensions. In particular, it can be applied to the electric field which, in the absence of magnetic source currents, solves the equation

(A17)$$\nabla \times \nabla \times {\bf E}({\bf r},\omega ) - { \epsilon _{\text{R}}}({\bf r},\omega )\frac{{{\omega ^2}}}{{{\text{c}^2}}}{\bf E}({\bf r},\omega ) = \text{i}\omega {\mu _0}{{\bf J}_{\text{s}}}({\bf r},\omega ).$$

Following a similar argumentation as that leading to Eq. (A2), we find that the solution to Eq. (A17) with the Silver–Müller radiation condition is

(A18)$${\bf E}({\bf r},\omega ) = \text{i}\omega {\mu _0}\int {{\bf G}^{\text{EE}}}({\bf r},{{\bf r}^\prime},\omega ){{\bf J}_{\text{s}}}({{\bf r}^\prime},\omega )\text{d}V^{\prime},$$

where the integral is over all space, provided the electric field Green tensor solves the equation

(A19)$$\nabla \times \nabla \times {{\bf G}^{\text{EE}}}({\bf r},{{\bf r}^\prime},\omega ) - { \epsilon _{\text{R}}}({\bf r},\omega )\frac{{{\omega ^2}}}{{{\text{c}^2}}}{{\bf G}^{\text{EE}}}({\bf r},{{\bf r}^\prime},\omega ) = \delta ({\bf r} - {{\bf r}^\prime}){\bf I},$$

where $ {\bf I} $ is the $ 3 \times 3 $ unit matrix, and each of the columns is required to satisfy the Silver–Müller radiation condition for the electric field in Eq. (19). Owing to the vector nature of the electric field in three dimensions, the electric field Green tensor must be a $ 3 \times 3 $ matrix in order to enable sufficient freedom to correctly describe the field resulting from a source with arbitrary and possibly position-dependent polarization through the superposition integral in Eq. (A18). Each of the three columns corresponds to one of the three principal orientations of a vectorial source in three dimensions. In Cartesian coordinates, for example, it follows from Eq. (A18) by setting $ {{\bf J}_{\text{s}}}({\bf r},\omega ) = \delta ({\bf r} - {{\bf r}^\prime}){J_{\text{s}}}(\omega ){{\hat{\bf e}}_{\alpha}}$ that each column ɑ of the electric field Green tensor is proportional to the electric field at $ {\bf r} $ due to a point source of magnitude $ {J_{\text{s}}}(\omega ) $ located at $ {{\bf r}^\prime} $ and pointing in the direction of the unit vector $ {\hat{\bf e}}_\alpha\in\{{\hat{\bf x}},{\hat{\bf y}},{\hat{\bf z}} $. We note that Eq. (A19) is consistent with the equation for the scalar Green function in Eq. (A3) for problems with sufficient symmetry. In the one-dimensional case, we may limit the discussion to the $ \hat{{\bf y}}\hat{{\bf y}} $ component of the tensor and define $ G(x,x^{\prime},\omega ) = \hat {{\bf y}} \cdot {{\bf G}^{\text{EE}}}(x,x^{\prime},\omega ) \cdot \hat{{\bf y}} $. Writing out the curl operator, the equation for this one-dimensional Green function reduces to the exact form of Eq. (A3).

Although arguably correct as written, Eq. (A18) has an ambiguity of great practical importance in cases where the electric field is evaluated in the source region, because the Green tensor diverges at this point. To explicitly address this problem in an unambiguous way, the equation is sometimes written in the form [168,169]

(A20)$${\bf E}({\bf r},\omega ) = \text{i}\omega {\mu _0}\int_{V - \delta V} {{\bf G}^{\text{EE}}}({\bf r},{{\bf r}^\prime},\omega ){{\bf J}_{\text{s}}}({{\bf r}^\prime},\omega )\text{d}V^{\prime} + \frac{{{{\bf J}_{\text{s}}}({\bf r},\omega )}}{{\text{i}\omega { \epsilon _0}}}{\bf L},$$

where the integral is taken over the extent of the source current except the possible point $ {\bf r} = {{\bf r}^\prime} $, which is to be avoided by use of a suitable exclusion volume $ \delta V $ and subsequently taking the limit $ \delta V \to 0 $. The choice of exclusion volume affects the explicit form of the source dyadic $ {\bf L} $, as discussed by Yaghjian [169]. In a large number of practical problems, the exclusion volume can be taken to be spherical, in which case $ {\bf L} = {\bf I}/3 $. Like any other vector field, the electric field can be decomposed into a transverse (divergence free) part with $ \nabla \cdot {{\bf E}_ \bot }({\bf r},\omega ) = 0 $ and a longitudinal (curl free) part with $ \nabla \times {{\bf E}_{||}}({\bf r},\omega ) = 0 $. The electric field Green tensor inherits this decomposition, so that we may in general write

(A21)$${{\bf G}^{\text{EE}}}({\bf r},{{\bf r}^\prime},\omega ) = {\bf G}_ \bot ^{\text{EE}}({\bf r},{{\bf r}^\prime},\omega ) + {\bf G}_\parallel ^{\text{EE}}({\bf r},{{\bf r}^\prime},\omega ).$$

A.2a. Electric Field Green Tensor in Homogeneous Media

Like the scalar analogue, the electric field Green tensor in homogeneous media can be calculated in different ways [165,170,171], and because of the translational invariance of the system, we expect that it depends only on the difference between the positions $ {\bf r} $ and $ {{\bf r}^\prime} $. Note, however, that because of the vector nature of the electric field, the direction between the two points becomes important. Defining $ {\bf R} = {\bf r} - {{\bf r}^\prime} $ and $ R = |{\bf R}| $, it can be written as [171]

(A22)$${\bf G}_{\text{B}}^{\text{EE}}({\bf R}) = - \frac{1}{{3k_{\text{B}}^2}}\delta ({\bf R}) + \left( {1 + \frac{{\nabla \nabla }}{{k_{\text{B}}^2}}} \right){g_{\text{B}}}(R,\omega ),$$

in which $ {k_{\text{B}}} = {n_{\text{B}}}\omega /\text{c} $, and

(A23)$${g_{\text{B}}}(R,\omega ) = \frac{{{\text{e}^{\text{i}{k_{\text{B}}}R}}}}{{4\pi R}}$$

is the scalar electric field Green function in the homogeneous material. This, in turn, can also be written in terms of the spherical Hankel function of the first kind and order zero as

(A24)$${g_{\text{B}}}(R,\omega ) = \frac{{\text{i}{k_{\text{B}}}}}{{4\pi }}{h_0}({k_{\text{B}}}R).$$

The delta-function term in Eq. (A22) is connected with the ambiguity in Eq. (A18) when $ {\bf r} = {{\bf r}^\prime} $. Although it appears naturally from the field expansion technique in Ref. [171], it does not show up in all derivations. Therefore, it is sometimes added by hand in order to ensure the validity of certain sum rules [170] based on integrals similar to Eq. (A18). The inclusion of the delta-function term thus provides an illustrative, if not entirely rigorous, connection between the two expressions for the electric field in Eqs. (A18) and (A20) in the case of a spherical exclusion volume. For completeness, we note that the transverse and longitudinal parts of the background Green tensor are  [171]

(A25)$${\bf G}_{\text{B} \bot }^{\text{EE}}({\bf r},{{\bf r}^\prime},\omega ) = \frac{{{\text{c}^2}}}{{4\pi { \epsilon _{\text{B}}}{\omega ^2}{R^3}}}\left[ {{\bf I} - 3\frac{{{\bf RR}}}{{{R^2}}}} \right] + \left( {1 + \frac{{\nabla \nabla }}{{k_{\text{B}}^2}}} \right){g_{\text{B}}}(R,\omega ),$$

(A26)$${\bf G}_{\text{B}\parallel }^{\text{EE}}({\bf r},{{\bf r}^\prime},\omega ) = - \frac{{{\text{c}^2}}}{{3{ \epsilon _{\text{B}}}{\omega ^2}}}\delta ({\bf R}) - \frac{{{\text{c}^2}}}{{4\pi { \epsilon _{\text{B}}}{\omega ^2}{R^3}}}\left[ {{\bf I} - 3\frac{{{\bf RR}}}{{{R^2}}}} \right],$$

respectively. Note that the defining equation for the electric field Green tensor in Ref. [171] is different from Eq. (A19) by a factor of $ - 1 $.

A.2b. Scattering Calculations Using the Green Tensor

In the vectorial case, the scattering integral in Eq. (A13) and the Lippmann–Schwinger equation in Eq. (A14) generalize to [169,172]

(A27)$$\begin{split}{{\bf E}_{\text{tot}}}({\bf r},\omega ) &= {{\bf E}_{\text{in}}}({\bf r},\omega ) + \frac{{{\omega ^2}}}{{{\text{c}^2}}}\int_{V - \delta V} {{\bf G}^{\text{EE}}}({\bf r},{{\bf r}^\prime},\omega )\Delta \epsilon ({{\bf r}^\prime},\omega ){{\bf E}_{\text{in}}}({{\bf r}^\prime},\omega )\text{d}V \\ &\quad- {\bf L}\frac{{\Delta \epsilon ({\bf r},\omega )}}{{{ \epsilon _{\text{B}}}}}{{\bf E}_{\text{in}}}({\bf r},\omega ),\end{split}$$

(A28)$$\begin{split}{{\bf E}_{\text{tot}}}({\bf r},\omega ) &= {{\bf E}_{\text{in}}}({\bf r},\omega ) + \frac{{{\omega ^2}}}{{{\text{c}^2}}}\int_{V - \delta V} {\bf G}_{\text{B}}^{\text{EE}}({\bf r},{{\bf r}^\prime},\omega )\Delta \epsilon ({{\bf r}^\prime},\omega ){{\bf E}_{\text{tot}}}({{\bf r}^\prime},\omega )\text{d}V \\&\quad- {\bf L}\frac{{\Delta \epsilon ({\bf r},\omega )}}{{{ \epsilon _{\text{B}}}}}{{\bf E}_{\text{tot}}}({\bf r},\omega ),\end{split}$$

respectively, where the volume integrals are taken over the volume of $ \Delta \epsilon ({\bf r}) \ne 0 $ except the possible point $ {\bf r} = {{\bf r}^\prime} $, which is to be avoided by use of a suitable exclusion volume $ \delta V $ as discussed in connection with Eq. (A20) above.

As noted above, each column of the electric field Green tensor is proportional to the electric field at $ {\bf r} $ due to a point source at $ {{\bf r}^\prime} $ pointing in one of the three principal directions. Therefore, one can also use an approach similar to Eq. (A28) to calculate the Green tensor of a scattering geometry defined by $ \Delta \epsilon ({\bf r},\omega ) $ by interpreting the background Green tensor as the incoming field. Including specifically the possibility of $ {\bf r} $ inside the scattering geometry, we can write the Dyson equation for the electric field Green tensor as

(A29)$$\begin{split}{{\bf G}^{\text{EE}}}({\bf r},{{\bf r}^\prime},\omega ) &= {\bf G}_{\text{B}}^{\text{EE}}({\bf r},{{\bf r}^\prime},\omega ) + \frac{{{\omega ^2}}}{{{\text{c}^2}}}\int_{V - \delta V} {\bf G}_{\text{B}}^{\text{EE}}({\bf r},{\bf s},\omega )\Delta \epsilon ({\bf s},\omega ){{\bf G}^{\text{EE}}}({\bf s},{{\bf r}^\prime},\omega )\text{d}V \!\!\! \\ &\quad - {\bf L}\frac{{\Delta \epsilon ({\bf r},\omega )}}{{{ \epsilon _{\text{B}}}}}{{\bf G}^{\text{EE}}}({\bf r},{{\bf r}^\prime},\omega ).\!\!\!\!\!\end{split}$$

Appendix B: Practical Convergence Studies

Contemporary modeling of electromagnetic scattering relies heavily on numerical solutions of partial differential equations. These numerical solutions naturally come with associated numerical errors, and it is therefore of considerable interest to have a systematic way of estimating the accuracy. In this appendix, we discuss how one can use the mathematical definitions of convergence and consistency to assess the convergence and ultimately assign an estimated error to a numerical calculation in a systematic way.

B.1. Convergence and Consistency

Any discretization will introduce some kind of parameter $ h $ controlling how fine the discretization is. The implicit assumption is that, as the size of $ h $ is varied to make the discretization finer and finer, the resulting calculated values become closer and closer, so that the sequence of results is convergent.

Convergence (Cauchy criterion): A sequence of numbers $ {s_1},{s_2},{s_3},\ldots,{s_n} $ is convergent if for any $ \epsilon \gt 0 $, there exists a number $ N $, so that $ |{s_m} - {s_n}| \lt \epsilon $ for all $ n,m \gt N $.

In practice, however, we require not only that a sequence is convergent, but also that it tends to the correct result—a requirement known as consistency.

Consistency: A sequence of numbers $ {s_1},{s_2},{s_3},\ldots,{s_n} $ in an approximation scheme is convergent with limit $ S $ if, for any $ \epsilon \gt 0 $, there exists a number $ N $, so that $ |S - {s_n}| \lt \epsilon $ for all $ n \gt N $. If $ S $ is the correct function value of the underlying continuous problem, the approximation scheme is said to be consistent.

The ideas behind the two definitions are reflected in different approaches for assessing the convergence of a given sequence, which we shall refer to as convergence and consistency studies. Convergence studies rely on the Cauchy criterion and therefore can be performed on any set of data. Consistency studies, on the other hand, can only be applied to cases in which the correct result is known. In such (rare) cases, one can directly plot the absolute or relative error as a function of the parameter(s) limiting the accuracy. In both cases, if the norm $ |S - {s_n}| $ or $ |{s_m} - {s_n}| $ decreases in a systematic way as we increase the parameter(s) limiting the accuracy, we shall assume that the sequence is convergent. As noted in the introduction, however, we should also be interested in somehow estimating the error in the calculations. While this is straightforward in cases where the result is known, it requires a bit of work and additional assumptions for convergence studies based on the Cauchy criterion. Ideally, one should use both approaches to assess the accuracy. Indeed, there can easily be systematic errors that only appear in direct comparison to known solutions—this is the case, for example, with reflections from PML boundaries. One would then first check for consistency by use of an auxiliary problem with comparable physical dimensions and materials and comparing it to a high-accuracy reference calculation to identify the parameter(s) limiting the accuracy. Only thereafter does it make sense to worry about the convergence properties and the accuracy of the actual problem at hand.

In practical calculations, one can vary a (sometimes rather large) number of parameters, and the error will depend stronger on some of these parameters than on others. A prominent example of a parameter limiting the accuracy is the size of the discrete triangles making up the calculation mesh in a three-dimensional BEM calculation, in which case the error is expected to vanish only in the limit of vanishing discretization size. For a generic parameter $ h $ and a convergent series of function values $ f(h) $ with limit $ {f_0} $, we can write the calculated value at any $ h $ as the correct value plus an error term as

(B1)$$f(h) = {f_0} + {\cal E}(h),$$

where $ {\cal E}(h) \to 0 $ for $ h \to 0 $. Assuming that one can identify the parameters limiting the accuracy, one can typically also identify the functional form of the error term—even if one does not know the correct limit $ {f_0} $. To this end, the observations in Subsections B.2 and B.3 below may be of help.

In cases where we are able to obtain a model for the error term, we can estimate the limiting value $ {f_0} $ by rewriting Eq. (B1) as

(B2)$${f_0} \approx f({h_{\text{min}}}) - {\cal E}({h_{\text{min}}}),$$

where $ {h_{\text{min}}} $ is the smallest value of $ h $ used in the calculations. In addition to the estimated value of $ {f_0} $, we shall generally use the absolute value of $ {\cal E}({h_{\text{min}}}) $ as a conservative estimate of the numerical error on the estimated value. We write the estimated numerical error in parenthesis immediately behind the digit(s) to which it pertains. A value of $ \pi \approx 3.15(2) $ thus signifies that there is an estimated error of $ \Delta \pi = \pm 0.02 $, so that we expect the true value of $ \pi $ to lie in the interval $ 3.13 \lt \pi \lt 3.17 $.

B.2. Polynomial Convergence

In many numerical calculations, the error tends to zero in a polynomial fashion for which we can write the dominant term as

(B3)$${\cal E}(h) = {{\cal E}_0}{h^\alpha },$$

where $ {{\cal E}_0} $ and $ \alpha $ are the initially unknown parameters characterizing the dominant polynomial behavior. Assuming Eq. (B1) to hold, we can eliminate the unknown correct limit $ {f_0} $ by considering the difference

(B4)$$f(h) - f(xh) = {\cal E}(h) - {\cal E}(xh)$$

(B5)$$ = {{\cal E}_0}{h^\alpha }[1 - {x^\alpha }],$$

for $ x \lt 1 $. Taking the logarithm on both sides of the equation, we then find that

(B6)$$\log \left( {f(h) - f(xh)} \right) = \alpha \log (h) + \log \left( {{{\cal E}_0}[1 - {x^\alpha }]} \right),$$

so that one can conveniently extract the exponent in the (dominant) polynomial error term by a simple fit of the logarithm of the differences as a function of the logarithm of $ h $. One can conveniently limit the numerical work by only calculating the function at exponentially spaced values as $ f(h),f(xh),f({x^2}h),f({x^3}h),\ldots $. Traditionally, $ x = 1/2 $ is often used, corresponding to successively halving the parameter $ h $.

B.3. Exponential Convergence

In some cases, the numerical error tends to zero for increasing $ p $ in an exponential fashion as

(B7)$${\cal E}(p) = {{\cal E}_0}{k^p},$$

for $ k \lt 1 $. To make this model fit with Eq. (B1), we can set $ h = {k^p} $. As in the case of a polynomial convergence, we can eliminate the unknown limit $ {f_0} $ by forming the difference

(B8)$$f(p) - f(p + q) = {{\cal E}_0}[{k^p} - {k^{p + q}}]$$

(B9)$$ = {{\cal E}_0}{k^p}[1 - {k^q}],$$

for $ q \gt 0 $. Taking the logarithm on both sides, we find

(B10)$$\log (f(p) - f(p + q)) = p\log (k) + \log ({{\cal E}_0}[1 - {k^q}]),$$

so that one can conveniently extract the base $ k $ and subsequently the parameter $ {{\cal E}_0} $ from a simple fit of the logarithm of the differences as function of $ p $.

Appendix C: Calculation Details

In this Appendix, we provide explicit calculation details that we have deemed too lengthy to include in the main part of the Tutorial.

C.1. Calculation Details for the Dielectric Barrier

The one-dimensional example of a dielectric barrier represents the special case of a dielectric slab in three dimensions, when considering only fields of perpendicular incidence. For definiteness, we consider electromagnetic waves moving along the $ x $ axis and choose the polarizations of the electric and magnetic fields to be along the $ y $ and $ z $ axes, respectively. As stated in Subsection 2.1, we consider a dielectric barrier of width $ L $ and with frequency-independent refractive index $ {n_{\text{R}}} $, which is centered around the position $ x = 0 $ in a background of constant refractive index $ {n_{\text{B}}} = \sqrt {{ \epsilon _{\text{B}}}} $. Writing $ {\bf E}({\bf r},\omega ) = E(x,\omega )\hat{{\bf y}} $ and $ {\bf H}({\bf r},\omega ) = H(x,\omega )\hat{{\bf z}} $, we find that the Silver–Müller radiation condition takes the form

(C1)$$E(x,\omega ) \to \sqrt {\frac{{{\mu _0}}}{{{ \epsilon _{\text{B}}}{ \epsilon _0}}}} H(x,\omega ),\quad x \to \infty .$$

The general solutions to the wave equation are known to be a sum of forward and backward propagating plane waves of the form $ E(x,\omega ) \propto \exp \{ \pm \text{i}\omega x/\text{c}\} $. Therefore, if the radiation condition is satisfied in the limit $ x \to \infty $, it must also be satisfied at all finite values of $ x $ outside the scattering region. Considering also positions left of the barrier, we can therefore rewrite the radiation condition as a boundary condition in the exact form in Eq. (6), which is satisfied by both the QNMs and the Green function for the dielectric barrier.

C.1a. QNMs of the Dielectric Barrier

To calculate the QNMs of the dielectric barrier, we start with the ansatz for the electric field QNM,

(C2)$${\tilde{\text{f}}_m}(x) = \left\{ {\begin{array}{*{20}{l}}{{A_m}{\text{e}^{ - \text{i}{n_{\text{B}}}{{\tilde k}_m}x}}}&\quad{x \lt - L/2}\\[3pt]{{B_m}{\text{e}^{\text{i}{n_{\text{R}}}{{\tilde k}_m}x}} + {C_m}{\text{e}^{ - \text{i}{n_{\text{R}}}{{\tilde k}_m}x}}}&\quad{ - L/2 \lt x \lt L/2}\\[3pt]{{D_m}{\text{e}^{\text{i}{n_{\text{B}}}{{\tilde k}_m}x}}}&\quad{L/2 \lt x}\end{array}} \right.,$$

where $ {\tilde k_m} = {\tilde \omega _m}/\text{c} $. Because of the symmetry of the problem, the solutions will be either even or odd with respect to the point $ x = 0 $. Therefore, we can set $ {C_m} = ( - {1)^m}{B_m} $ and $ {D_m} = ( - {1)^m}{A_m} $. Moreover, since the QNMs are eigenfunctions, we can introduce a scaling of choice by setting $ {B_m} = 1 $. Finally, for a given $ m $, we can then calculate $ {A_m} $ using the continuity at $ x = - L/2 $ to find the general form for the electric field QNM in Eq. (2). From the Maxwell curl equations, it now follows that the corresponding magnetic field QNM must be of the form

(C3)$${\tilde{\text{g}}_m}(x) = \sqrt {\frac{{{ \epsilon _0}}}{{{\mu _0}}}} \left\{{\begin{array}{*{20}{l}}{ - {n_{\text{B}}}{A_m}{\text{e}^{ - \text{i}{n_{\text{B}}}{{\tilde k}_m}x}}}&\quad{x \lt - L/2}\\[3pt]{{n_{\text{R}}}\left( {{\text{e}^{\text{i}{n_{\text{R}}}{{\tilde k}_m}x}} - {{( - 1)}^m}{\text{e}^{ - \text{i}{n_{\text{R}}}{{\tilde k}_m}x}}} \right)}&\quad{ - L/2 \lt x \lt L/2}\\[3pt]{{{( - 1)}^m}{n_{\text{B}}}{A_m}{\text{e}^{\text{i}{n_{\text{B}}}{{\tilde k}_m}x}}}&\quad {L/2 \lt x}\end{array}} \right..$$

To calculate the complex QNM frequencies, we use the requirement of continuity of both $ \tilde{\text{f}}_m (x) $ and $ \tilde{\text{g}}_m (x) $. At $ x = - L/2 $, we find the conditions

(C4)$${A_m}{\text{e}^{\text{i}{n_{\text{B}}}{{\tilde k}_m}L/2}} = {\text{e}^{ - \text{i}{n_{\text{R}}}{{\tilde k}_m}L/2}} + {( - 1)^m}{\text{e}^{\text{i}{n_{\text{R}}}{{\tilde k}_m}L/2}}$$

and

(C5)$$ - {n_{\text{B}}}{A_m}{\text{e}^{\text{i}{n_{\text{B}}}{{\tilde k}_m}L/2}} = {n_{\text{R}}}\left( {{\text{e}^{ - \text{i}{n_{\text{R}}}{{\tilde k}_m}L/2}} - {{( - 1)}^m}{\text{e}^{\text{i}{n_{\text{R}}}{{\tilde k}_m}L/2}}} \right),$$

and combining the two, we can eliminate $ {A_m} $ to find the condition

(C6)$${({n_{\text{R}}} + {n_{\text{B}}})^2}{\text{e}^{ - \text{i}{n_{\text{R}}}{{\tilde k}_m}L}} = ({n_{\text{R}}} - {n_{\text{B}}}{)^2}{\text{e}^{\text{i}{n_{\text{R}}}{{\tilde k}_m}L}}.$$

Rearranging the terms by expanding the squares, this relation can be rewritten in the form of Eq. (4) for which the solutions are given in Eq. (5). Up to a normalization factor, this completely specifies the QNMs of the dielectric barrier.

C.1b. Electric Field Green Function for the Dielectric Barrier

To calculate $ G(x,x^{\prime},\omega ) $, we follow an approach similar to that adopted in Appendix A.1a for the background Green function. To this end, we choose a fixed $ x^{\prime} $ within the barrier and expand the function in forward and backward traveling waves in the various regions defined by $ x^{\prime} $ and the boundaries at $ x = \pm L/2 $. For the case of $ x^{\prime} $ inside the resonator, this gives the ansatz

(C7)$$G(x,x^{\prime},\omega ) = \left\{ {\begin{array}{*{20}{l}}{A(x^{\prime}){\text{e}^{ - \text{i}{n_{\text{B}}}kx}}}&\quad{x \lt - L/2}\\{B(x^{\prime}){\text{e}^{\text{i}{n_{\text{R}}}kx}} + C(x^{\prime}){\text{e}^{ - \text{i}{n_{\text{R}}}kx}}}&\quad{ - L/2 \lt x \lt x^{\prime}}\\{D(x^{\prime}){\text{e}^{\text{i}{n_{\text{R}}}kx}} + E(x^{\prime}){\text{e}^{ - \text{i}{n_{\text{R}}}kx}}}&\quad{x^{\prime} \lt x \lt L/2}\\{F(x^{\prime}){\text{e}^{\text{i}{n_{\text{B}}}kx}}}&\quad{L/2 \lt x}\end{array}} \right.,$$

where $ k = \omega /\text{c} $, and for which we must determine the six expansion coefficients. Demanding continuity of the Green function at $ x = \pm L/2 $, we can express $ A(x^{\prime}) $ and $ F(x^{\prime}) $ as

(C8)$$A(x^{\prime}) = {\text{e}^{\text{i}({n_{\text{R}}} - {n_{\text{B}}})kL/2}}C(x^{\prime}) + {\text{e}^{ - \text{i}({n_{\text{R}}} + {n_{\text{B}}})kL/2}}B(x^{\prime}),$$

(C9)$$F(x^{\prime}) = {\text{e}^{\text{i}({n_{\text{R}}} - {n_{\text{B}}})kL/2}}D(x^{\prime}) + {\text{e}^{ - \text{i}({n_{\text{R}}} + {n_{\text{B}}})kL/2}}E(x^{\prime}),$$

and by combining with the requirement of differentiability, we can eliminate $ A(x^{\prime}) $ and $ F(x^{\prime}) $ to find the relations

(C10)$$C(x^{\prime}) = \alpha {\text{e}^{ - \text{i}{n_{\text{R}}}L}}B(x^{\prime}),$$

(C11)$$D(x^{\prime}) = \alpha {\text{e}^{ - \text{i}{n_{\text{R}}}L}}E(x^{\prime}),$$

in which $ \alpha = ({n_{\text{R}}} + {n_{\text{B}}})/({n_{\text{R}}} - {n_{\text{B}}}) $. To find expressions for $ B(x^{\prime}) $ and $ E(x^{\prime}) $, we combine the requirement of continuity at $ x = x^{\prime} $ and the additional requirement obtained by integrating the defining equation around the point $ x = x^{\prime} $ as in Eq. (A8). In this way, we find that $ B(x^{\prime}) $ and $ E(x^{\prime}) $ solve the equation

(C12)$$\left[ {\begin{array}{*{20}{c}}{{\text{e}^{\text{i}{n_{\text{R}}}kx^{\prime}}}}&{ - \alpha {\text{e}^{ - \text{i}{n_{\text{R}}}kL}}{\text{e}^{\text{i}{n_{\text{R}}}kx^{\prime}}}}\\{\alpha {\text{e}^{ - \text{i}{n_{\text{R}}}kL}}{\text{e}^{ - \text{i}{n_{\text{R}}}kx^{\prime}}}}&{ - {\text{e}^{ - \text{i}{n_{\text{R}}}kx^{\prime}}}}\end{array}} \right]\left[ {\begin{array}{*{20}{c}}{B(x^{\prime})}\\{E(x^{\prime})}\end{array}} \right] = - \frac{\text{i}}{{2{n_{\text{R}}}k}}\left[ {\begin{array}{*{20}{c}}1\\{ - 1}\end{array}} \right],$$

for which we can write the solutions as

(C13)$$B(x^{\prime}) = - \frac{\text{i}}{{2{n_{\text{R}}}k}}\frac{{{\text{e}^{\text{i}{n_{\text{R}}}k(L - x^{\prime})}} + \alpha {\text{e}^{\text{i}{n_{\text{R}}}kx^{\prime}}}}}{{{\text{e}^{\text{i}{n_{\text{R}}}kL}} - {\alpha ^2} {\text{e}^{ - \text{i}{n_{\text{R}}}kL}}}},$$

(C14)$$E(x^{\prime}) = - \frac{\text{i}}{{2{n_{\text{R}}}k}}\frac{{{\text{e}^{\text{i}{n_{\text{R}}}k(L + x^{\prime})}} + \alpha {\text{e}^{ - \text{i}{n_{\text{R}}}kx^{\prime}}}}}{{{\text{e}^{\text{i}{n_{\text{R}}}kL}} - {\alpha ^2} {\text{e}^{ - \text{i}{n_{\text{R}}}kL}}}}.$$

Inserting in Eq. (C7), this completely specifies the electric field Green function for the dielectric barrier in the case where at least one of the two points $ x $ or $ x^{\prime} $ are inside the material.

From the expansion coefficients in Eqs. (C13) and (C14), we can find the poles of the Green function as the solutions to the equation

(C15)$${\text{e}^{\text{i}{n_{\text{R}}}kL}} = {\alpha ^2} {\text{e}^{ - \text{i}{n_{\text{R}}}kL}},$$

which is identical to the condition for the QNM resonance frequencies in Eq. (C6). In addition, the Green function has a pole at $ k = 0 $, which has consequences for the QNM expansions, as pointed out in Subsection 4.1a.

C.1c. Region of Consistency for the Dielectric Barrier

From the analytical form of the Green function, we can calculate the boundaries of the region of consistency for the dielectric barrier by analyzing the analytical continuation of the Green tensor as discussed in Subsection 4.3. For $ x^{\prime} $ inside the resonator, and $ x \gt L/2 $, we focus initially on the functional form of $ E(x^{\prime}) $ as given in Eq. (C14) and investigate the behavior for complex values of $ k $ in the lower half of the complex plane.

Simplifying the expression for $ E(x^{\prime}) $ by the factor $ \exp \{ \text{i}{n_{\text{R}}}kL\} $, we can rewrite it as

(C16)$$E(x^{\prime}) = R(k)\frac{{{\text{e}^{\text{i}{n_{\text{R}}}kx^{\prime}}} + \alpha {\text{e}^{\text{i}{n_{\text{R}}}k(x^{\prime} - L)}}}}{{1 - {\alpha ^2} {\text{e}^{ - 2\text{i}{n_{\text{R}}}kL}}}},$$

where $ R(k) = - \text{i}/2{n_{\text{R}}}k $ is unimportant in determining the region of consistency. The second term in the denominator tends to zero in an exponential fashion as $ k $ is varied along any straight line downwards in the lower half of the complex plane, and the first term in the numerator will be decisive in determining the region of consistency. Combining Eqs. (C11) and (C9), we can write

(C17)$$F(x^{\prime}){\text{e}^{i{n_{\text{B}}}kx}} = (\alpha + 1){\text{e}^{ - \text{i}({n_{\text{R}}} + {n_{\text{B}}})kL/2}}E(x^{\prime}),$$

and inserting Eq. (C16), we find that the numerator tends to zero when $ k $ is varied along any straight line downwards in the complex plane, provided

(C18)$${n_{\text{R}}}(x^{\prime} - L/2) + {n_{\text{B}}}(x - L/2) \lt 0.$$

As the position of $ x^{\prime} $ varies towards the boundary from the inside, the region of consistency closes in on the boundary from the outside. For $ x^{\prime} = L/2 $, a detailed analysis shows that $ x = L/2 $ is not in the region of consistency. Due to the symmetry of the problem, we can immediately infer a similar behavior for $ x \lt L/2 $, and rewriting slightly, we find the general condition

(C19)$$ - \frac{L}{2} - \frac{{{n_{\text{R}}}}}{{{n_{\text{B}}}}}\left( {\frac{L}{2} + x^{\prime}} \right) \lt x \lt \frac{L}{2} + \frac{{{n_{\text{R}}}}}{{{n_{\text{B}}}}}\left( {\frac{L}{2} - x^{\prime}} \right).$$

C.2. Calculation Details for the Plasmonic Dimer

In the VIE calculations for the plasmonic dimer, the spherical wave functions constitute a complete basis within the individual spheres. Therefore, we can assume the expansion to converge to the correct solution in the limit $ {l_{\text{max}}} \to \infty $, provided the numerical solution of the matrix problem in Eq. (8) does not itself introduce errors. In practice, there will be truncation errors associated with the representation of the data in finite precision, but these errors are expected to be at least 3 orders of magnitude smaller than the estimated errors below, so we will ignore them in this analysis. The calculations were performed by assuming that the material parameters of the Drude permittivity are given exactly by $ \hbar {\omega _{\text{p}}} = 7.9\,\text{eV} $ and $ \hbar \gamma = 0.06\,\text{eV} $, which are similar to values that have been found to provide a reasonable description for gold around the visible frequency range, see for example Refs. [173,174]. We note, however, that the Drude model is mostly a convenient abstraction that allows us to clearly and unambiguously define the mathematical problem. For high-accuracy calculations of real material systems, it may not be adequate, but one can often obtain good approximations in certain frequency ranges by working with material models made up from several resonance contributions as in the Drude–Lorentz model,

(C20)$${ \epsilon _{\text{R}}}(\omega ) = { \epsilon _{\text{const}}} + \sum\limits_j \frac{{A_j^2}}{{\omega _j^2 - {\omega ^2} - \text{i}{\gamma _j}\omega }},$$

where $ { \epsilon _{\text{const}}} $ is a constant background, which is typically larger than one, and $ {A_j} $ are amplitudes of the material resonances at the complex frequencies $ {\omega _j} - \text{i}{\gamma _j} $. Compared to Eq. (1), the Drude model is a special case of the Drude–Lorentz model with $ { \epsilon _{\text{const}}} = 1 $ and a single resonance with $ {\omega _j} = 0 $. Variants and generalizations of the Drude–Lorentz response have been used to model a wide variety of materials by fitting the unknown parameters in the model to tabulated material data such as those provided by Johnson and Christy [175].

For the calculation of $ {\tilde \omega _1} $, Fig. 4 shows, as a function of the cutoff parameter $ {l_{\text{max}}} $, the logarithm of $ {D_{ + 1}}({l_{\text{max}}}) = {\tilde \omega _{\text{num}}}({l_{\text{max}}}) - {\tilde \omega _{\text{num}}}({l_{\text{max}}} + 1) $. To a good approximation, the data points for both the real and the imaginary parts in Fig. 4 fall on a straight line, indicating an exponential convergence in both cases, cf. the discussion in Appendix B.3. For the other QNM resonance frequencies of interest, we find a similar behavior. Assuming that the finite value of $ {l_{\text{max}}} $ is the dominating source of error in the calculation, we can use fitted values of $ k $ and $ {{\cal E}_0} $ in Eq. (B7) to obtain a model for $ {\cal E}({l_{\text{max}}}) $. Rewriting Eq. (B1) by setting $ {f_0} = {\tilde \omega _n} $ and $ f(p) = {\tilde \omega _{\text{num}}}(p) $, we can then use the largest value of $ {l_{\text{max}}} = 8 $ to estimate the true value as

(C21)$${\tilde \omega _n} \approx {\tilde \omega _{\text{num}}}(8) - {{\cal E}_0}{k^8}.$$

The calculated values are listed in Table 1 along with the estimated errors $ \Delta{\tilde{\omega}_n} = |\Delta {\omega _n} - \text{i}\Delta {\gamma _n}| $, which grow as a function of $ n $, because all calculations were performed with a fixed number of spherical wave functions set by $ {l_{\text{max}}} = 8 $. The QNM wave functions of interest grow in complexity with increasing values of $ n $, and this translates into a slower rate of convergence and hence an increase in numerical error for a given value of $ {l_{\text{max}}} $. The associated generalized effective mode volumes of the five QNMs are listed in Table 2. They were all calculated using integration around the complex resonance frequency as described in Subsection 3.4c. with postprocessing of the data to estimate the numerical error as described above.

C.2a. Normalization Using the Boundary Element Method

The BEM enables calculations of the electric field Green function at complex frequencies. Therefore, we can use the normalization procedure in Subsection 3.4a to calculate the inverse norm, and by comparing it to the reference value in Eq. (81), we can investigate the convergence of the numerical scheme for the exact same meshes used in the investigation of the resonance frequency in Fig. 5. The resulting relative error in the calculated values is shown as a function of average side length in a double logarithmic plot in Fig. 32. Even though the overall convergence does not appear as convincing as for the resonance frequency in Fig. 5, we note that the data points corresponding to the finest meshes do appear to fall on a relatively straight line; by fitting to the last four data points, we find that the slope of this line is also close to one, albeit slightly below that of the corresponding fit to the resonance frequency data. Nevertheless, in view of the similar convergence rate observed for the resonance frequency, we can assume that the overall convergence will eventually tend to the same first-order polynomial behavior, as one would expect.

 

Figure 32. Convergence analysis for the BEM calculation of the inverse norm $ \langle \langle {\tilde {\underline {\bf F} }_1}({\bf r})|{\tilde {\underline {\bf F} }_1}({\bf r})\rangle \rangle _1^{ - 1} $ showing the relative error $ {{\cal E}_{\text{BEM}}} = |\langle \langle {\tilde {\underline {\bf F} }_1}({\bf r})|{\tilde {\underline {\bf F} }_1}({\bf r})\rangle \rangle _{1,\text{BEM}}^{ - 1} - \langle \langle {\tilde {\underline {\bf F} }_1}({\bf r})|{\tilde {\underline {\bf F} }_1}({\bf r})\rangle \rangle _{1,\text{ref}}^{ - 1}|/|\langle \langle {\tilde {\underline {\bf F} }_1}({\bf r})|{\tilde {\underline {\bf F} }_1}({\bf r})\rangle \rangle _{1,\text{ref}}^{ - 1}| $ as a function of average side length $ h $ of the triangular mesh. The reference value $ \langle \langle {\tilde {\underline {\bf F} }_1}({\bf r})|{\tilde {\underline {\bf F} }_1}({\bf r})\rangle \rangle _{1,\text{ref}}^{ - 1} $ is given in Eq. (81), while $ {{\cal E}_{\text{BEM}}} = |\langle \langle {\tilde {\underline {\bf F} }_1}({\bf r})|{\tilde {\underline {\bf F} }_1}({\bf r})\rangle \rangle _{1,\text{BEM}}^{ - 1} $ is the value obtained with the BEM approach. Circles indicate the calculated values, and the black line shows a fit to the last four data points.

Download Full Size | PDF

C.3. Dispersive Materials

Upon combining the Maxwell curl equations with Eq. (45) for the dispersive Drude material model, the governing matrix equation for the QNMs takes the form

(C22)$$\left[ {\begin{array}{*{20}{c}}0&{{{[{ \epsilon _0}]}^{ - 1}}\nabla \times }&{ - {{[{ \epsilon _0}]}^{ - 1}}}\\{ - {{[{\mu _0}]}^{ - 1}}\nabla \times }&0&0\\{{ \epsilon _0}\omega _{\text{p}}^2({\bf r})}&0&{ - \gamma ({\bf r})}\end{array}} \right]\left[ {\begin{array}{*{20}{c}}{{\bf E}({\bf r},\omega )}\\{{\bf H}({\bf r},\omega )}\\{{\bf J}({\bf r},\omega )}\end{array}} \right] = - \text{i}\omega \left[ {\begin{array}{*{20}{c}}{{\bf E}({\bf r},\omega )}\\{{\bf H}({\bf r},\omega )}\\{{\bf J}({\bf r},\omega )}\end{array}} \right],$$

with $ {\bf J}({\bf r},\omega ) = 0 $ at positions outside the material. In the Drude model, the current density is directly proportional to the electric field, so we do not introduce additional boundary conditions for $ {\bf J}({\bf r},\omega ) $ at the boundary of the material. The weight function in the inner product is extended to $ {\bf W} = \text{diag}\{ { \epsilon _0},{\mu _0},[{ \epsilon _0}\omega _{\text{p}}^2{]^{ - 1}}\} $, and the arguments of Subsection 3.2 can be repeated to find that if $ {\tilde {\underline {\bf F} }_m}({\bf r}) = [\tilde{{\bf f}}_m ({\bf r}),\tilde{{\bf g}}_m ({\bf r}),\tilde{{\bf j}}_m ({\bf r}{)]^{\text{T}}} $ is a solution to Eq. (25), where $ {\bf D} $ is the matrix in Eq. (C22) and the electromagnetic fields obey the Silver–Müller radiation condition, then $ \tilde{\underline {\bf F} }_m^{\ddagger }({\bf r}) = [\tilde{{\bf f}}_m ({\bf r}), - \tilde{{\bf g}}_m ({\bf r}), - \tilde{{\bf j}}_m ({\bf r}{)]^{\text{T}}} $ is a solution to Eq. (35), in which

(C23)$${\underline {\underline {\bf D} } ^{\ddagger }} = \left[ {\begin{array}{*{20}{c}}0&{ - {{[{ \epsilon _0}]}^{ - 1}}\nabla \times }&{{{[{ \epsilon _0}]}^{ - 1}}} \\ {{[\mu_0]}^{ - 1}\nabla \times }&0&0\\{ - { \epsilon _0}\omega _{\text{p}}^2({\bf r})}&0&{ - \gamma ({\bf r})}\end{array}} \right],$$

and the electromagnetic fields obey the adjoint radiation condition. The projection operator in Eq. (39) can now be generalized as

(C24)$$\begin{split}\langle \langle \tilde {\underline {\bf F} }_m^{\ddagger}({\bf r})|\underline {\bf F} ({\bf r},\omega )\rangle \rangle &= {\frac{1}{{2{ \epsilon _0}}}\int_V } {\big[{ \epsilon _0}{\tilde{{\bf f}}_m}({\bf r}) \cdot {\bf E}({\bf r},\omega ) - {\mu _0}{\tilde{{\bf g}}_m}({\bf r}) \cdot {\bf H}({\bf r},\omega )}\\&\quad- {\tilde{{\bf j}}_m}({\bf r}){\bf J}({\bf r},\omega )/{ \epsilon _0}\omega _{\text{p}}^2\big] \text{d}V\\ &\; + \frac{\text{i}}{{2{ \epsilon _0}({{\tilde \omega }_m} - \omega )}}\int_{\partial V} \big[{\bf E}({\bf r},\omega ) \times \tilde{\bf g}_m ({\bf r}) - \tilde{\bf f}_m ({\bf r}) \times {\bf H}({\bf r},\omega )\big] \cdot {\bf n} \text{d}A,\end{split}$$

and the rest of the arguments in Subsection 3.2 remain valid, as do the arguments surrounding the projection properties in Subsection 3.3. Note that the surface integral is unaffected, because the current density vanishes outside the material.

It is illustrative to rewrite the projection operator in terms of the electromagnetic fields only. From Eq. (C22), it is clear that the current density is directly proportional to the electric field,

(C25)$${\bf J}({\bf r},\omega ) = { \epsilon _0}\frac{{\omega _{\text{p}}^2({\bf r})}}{{\gamma ({\bf r}) - \text{i}\omega }}{\bf E}({\bf r},\omega ).$$

Inserting this expression in Eq. (C24), we can write the projection operator in the compact form of Eq. (46).

 

Figure 33. Schematic showing an electromagnetic resonator in the form of two spheres within an integration volume $ V $. For calculations of the projection of a field onto the QNMs, the additional integration through the volume $ V^{\prime}_{\text{B}} $ does not contribute.

Download Full Size | PDF

C.4. Independence of Integration Volume in Eq. (39)

Starting from Eq. (39) and considering any part of the volume with $ { \epsilon _{\text{R}}}({\bf r}) = { \epsilon _{\text{B}}} $, denoted by $ {V^{\prime}_{\text{B}}} $, we show below that the contribution to the entire integral from this part vanishes.

Considering two volumes $ V $ and $ V + {V^{\prime}_{\text{B}}} $, as illustrated in Fig. 33, the contribution to the projection operation from integration throughout $ {V^{\prime}_{\text{B}}} $ can be written as

(C26)$$\begin{split}{I_{{{V^{\prime}}_{\text{B}}}}}& = \frac{1}{{2{ \epsilon _0}}}\int_{{{V^{\prime}}_{\text{B}}}} [{ \epsilon _0}{ \epsilon _{\text{B}}}\tilde{{\bf f}}_m({\bf r}) \cdot {\bf E}({\bf r},\omega ) - {\mu _0}\tilde{{\bf g}}_m({\bf r}) \cdot {\bf H}({\bf r},\omega )]\text{d}{V^{\prime}_{\text{B}}} \\ &\quad + \frac{\text{i}}{{2{ \epsilon _0}({{\tilde \omega }_m} - \omega )}}\int_{\partial {{V^{\prime}}_{\text{B}}}} \left[ {{\bf H}({\bf r},\omega ) \times {\tilde{\bf f}}_m ({\bf r}) + {\bf E}({\bf r},\omega ) \times {\tilde{\bf g}}_m ({\bf r})} \right] \cdot \hat{{\bf n}} \text{d}A.\end{split}$$

Considering the first term in the surface integral, we can express the magnetic field in terms of the electric field using the curl equation and apply Green’s vector theorem of the first kind,

(C27)$$\int_V \Bigg[\nabla \times {\bf P}] \cdot [\nabla \times {\bf Q}] - {\bf P} \cdot \nabla \times \nabla \times {\bf Q}\Bigg] \text{d}V = \int_{\partial V} \hat{{\bf n}} \cdot \left[ {{\bf P} \times \nabla \times {\bf Q}} \right] \text{d}A,$$

to find

(C28)$$\begin{split}\int_{\partial {{V^{\prime}}_{\text{B}}}} \left[ {{\bf H}({\bf r},\omega ) \times {\tilde{{\bf f}}}_m ({\bf r})} \right] \cdot \hat{{\bf n}} \text{d}A &= \frac{\text{i}}{{{\mu _0}\omega }}\int_{{{V^{\prime}}_{\text{B}}}} \big[[\nabla \times \tilde{{\bf f}}_m ({\bf r})] \cdot [\nabla \times {\bf E}({\bf r},\omega )] \\ &\quad - \tilde{{\bf f}}_m ({\bf r}) \cdot \nabla \times \nabla \times {\bf E}({\bf r},\omega)\big]\text{d}V.\end{split}$$

Rewriting the first term using the curl equation, and the second term using the electric field wave equation, we can further simplify the expression as

(C29)$$\begin{split}\int_{\partial {{V^{\prime}}_{\text{B}}}} \left[ {{\bf H}({\bf r},\omega ) \times {\tilde{{\bf f}}}_m ({\bf r})} \right] \cdot \hat{{\bf n}} \text{d}A &= - \text{i}\int_{{{V^{\prime}}_{\text{B}}}} [{\tilde \omega _m}{\mu _0}\tilde{{\bf g}}_m ({\bf r}) \cdot {\bf H}({\bf r},\omega ) \\&\quad + \omega { \epsilon _0}{ \epsilon _{\text{B}}}\tilde{{\bf f}}_m ({\bf r}) \cdot {\bf E}({\bf r},\omega )] \text{d}V.\end{split}$$

In a similar fashion, we can rewrite the second term in the surface integral as

(C30)$$\begin{split}\int_{\partial {{V^{\prime}}_{\text{B}}}} \left[ {{\bf E}({\bf r},\omega ) \times {\tilde{{\bf g}}_m}({\bf r})} \right] \cdot \hat{{\bf n}} \text{d}A &= \text{i}\int_{{{V^{\prime}}_{\text{B}}}} [\omega {\mu _0}{\bf H}({\bf r},\omega ) \cdot \tilde{{\bf g}}_m ({\bf r}) \\ &\quad + {\tilde \omega _m}{ \epsilon _0}{ \epsilon _{\text{B}}}{\bf E}({\bf r},\omega ) \cdot \tilde{{\bf f}}_m ({\bf r})] \text{d}V.\end{split}$$

Finally, inserting the previous expressions in Eq. (C26), we find that the contributions from the surface integrals exactly cancel the volume integral, so that

(C31)$${I_{\partial {{V^{\prime}}_{\text{B}}}}} = 0.$$

We conclude that the regions outside the scatterers do not contribute to the value of the integral in Eq. (39), so one is free to choose any integration volume of convenience as long as it contains the scatterers for which $ { \epsilon _{\text{R}}}({\bf r}) \ne { \epsilon _{\text{B}}} $.

C.5. Region of Consistency for a Spherical Resonator

To investigate the region of consistency for the sphere, we will focus initially on the scalar electric field Green function in three dimensions, which solves the equation

(C32)$${\nabla ^2}g({\bf r},{{\bf r}^\prime},\omega ) + \frac{{{\omega ^2}}}{{{c^2}}}{ \epsilon _{\text{R}}}({\bf r},\omega )g({\bf r},{{\bf r}^\prime},\omega ) = - \delta ({\bf r} - {{\bf r}^\prime}).$$

As we shall see, the arguments leading to the region of consistency in the scalar case are identical to those in the vectorial case. The scalar Green function fulfills the Dyson equation

(C33)$$g({{\bf r}^\prime},{\bf r},\omega ) = {g_{\text{B}}}({{\bf r}^\prime},{\bf r},\omega ) + \frac{{{\omega ^2}}}{{{c^2}}}\int {g_{\text{B}}}({{\bf r}^\prime},{\bf s},\omega )\Delta \epsilon ({\bf s},\omega )g({\bf s},{\bf r},\omega )\text{d}V,$$

in which $ {g_{\text{B}}}({{\bf r}^\prime},{\bf r},\omega ) $ is the Green function of the homogeneous background material with permittivity $ { \epsilon _{\text{B}}} = n_{\text{B}}^2 $, and $ \Delta \epsilon ({\bf r},\omega ) = \epsilon ({\bf r},\omega ) - { \epsilon _{\text{B}}} $. For the present calculations, we shall assume that the material of the sphere is homogeneous, so that $ \Delta \epsilon ({\bf r},\omega ) = \Delta \epsilon (\omega ) $ inside the sphere, and zero outside.

As noted in Eq. (A24), The Green function in the homogeneous background material can be written in terms of the spherical Hankel function of the first kind. Therefore, given that we are ultimately interested in the limiting behavior of the Green function $ g({{\bf r}^\prime},{\bf r}, - \text{i}\xi ) $ in the limit $ \xi \to \infty $, we will make use of the fact that the spherical Hankel functions of the first kind, in general, can be written in the form [134]

(C34)$${h_n}(z) = {R_n}(z){\text{e}^{\text{i}z}},$$

where [134]

(C35)$${R_n}(z) = \frac{{{\text{i}^{ - (n + 1)}}}}{z}\sum\limits_{k = 0}^n \frac{{(n + k)!}}{{k!\Gamma (n - k + 1)( - 2\text{i}z{)^k}}}$$

is a rational function of polynomials; $ {R_1}(z) = - (z + \text{i})/{z^2} $, for example. Similarly, since the spherical bessel function $ {j_n}(z) $ can be written in terms of the spherical Hankel functions of the first and second kind, we can write

(C36)$${j_n}(z) = {P_n}(z){\text{e}^{\text{i}z}} + {Q_n}(z){\text{e}^{ - \text{i}z}},$$

where $ {P_n}(z) $ and $ {Q_n}(z) $ are rational functions of polynomials. Clearly, when $ z $ is varied in the direction of the negative imaginary axis, terms with $ \exp \{ \text{i}z\} $ will diverge (exponentially). As we will see, the exponentially divergent terms will be decisive in determining the region of consistency, whereas the terms with $ \exp \{ - \text{i}z\} $ will be subleading in the calculations below.

 

Figure 34. Schematic of a single homogeneous sphere with coordinates used for the derivation of the region of consistency. In the calculations using the Dyson equation, the background Green functions, which depend on the distances $ |{\bf r} - {{\bf r}^\prime}| $ and $ |{{\bf r}^\prime} - {\bf s}| $ (indicated by dashed lines), are reformulated in terms of vectors defined with respect to the center of the sphere $ {\cal O} $.

Download Full Size | PDF

We consider the case of a single spherical and homogeneous resonator with radius $ R $. For our analysis, $ {{\bf r}^\prime} $ is inside the sphere, and $ {\bf r} $ is outside, as illustrated in Fig. 34. The integral is over the volume of the sphere, and the general approach will therefore be to express the Green functions in terms of spherical Bessel and Hankel functions defined with respect to the center of the sphere. For $ s \gt r^{\prime} $, for example, we can write the spherical Hankel function as [132]

(C37)$${h_0}({n_{\text{B}}}k|{{\bf r}^\prime} - {\bf s}|) = 4\pi \sum\limits_{l,m} {h_l}({n_{\text{B}}}ks){\left\{ {Y_l^m(\hat{{\bf s}})} \right\}^*}{j_l}({n_{\text{B}}}kr^{\prime})Y_l^m(\hat{{\bf r}}^{\prime}),$$

where $ k = \omega/\text{c} $, and $ Y_l^m(\hat{{\bf v}}) $ denotes the spherical harmonic of order $ l,m $ evaluated at the direction of the vector $ {\bf v} $. Even if we do not know the explicit expression for $ g({{\bf r}^\prime},{\bf r},\omega ) $ at this point, we know that we can expand it in terms of spherical Bessel functions defined with respect to the center of the sphere as [46]

(C38)$$g({{\bf r}^\prime},{\bf r},\omega ) = \sum\limits_{\nu ,\mu } {j_n}({n_{\text{R}}}kr^{\prime})Y_\nu ^\mu (\hat {\bf r}^{\prime})\,{\alpha _{\nu ,\mu }}({\bf r},\omega ),$$

where the unknown expansion coefficients $ {\alpha _{\nu ,\mu }}({\bf r},\omega ) $ depend on both $ {\bf r} $ and $ \omega $ as indicated.

Inserting now Eqs. (C37) and (C38), we can express the integral in Eq. (C33) in terms of spherical wave functions defined with respect to the center of the sphere. For the angular part of the integration, the orthogonality relation [132]

(C39)$$\int_\Omega Y_n^m(\hat{{\bf s}}){\left\{ {Y_\nu ^\mu (\hat {{\bf s}})} \right\}^*}\text{d}\Omega = {\delta _{n\nu }}{\delta _{m\mu }},$$

where $ \Omega $ denotes the surface of the sphere, leads to considerable simplifications. For the integration over $ s $, we rewrite the spherical Bessel function as in Eq. (C36). Looking at the integrand of this expression, the limiting behavior is given by the terms in the integrand with $ s \gt r^{\prime} $ and divergent exponential factors of the form $ \exp \{ \text{i}{k_{\text{B}}}s\} \exp \{ \text{i}{n_{\text{R}}}{k_{\text{B}}}s\} $. In particular, by treating the integral as a Riemann sum, we can appreciate that there will be contributions to the integral with exponential factors of the form $ \exp \{ \text{i}{k_{\text{B}}}R\} \exp \{ \text{i}{n_{\text{R}}}{k_{\text{B}}}R\} $, and since $ s \le R $, terms with these factors will be the fastest growing terms. With this argumentation, we write the resulting expression in the form

(C40)$$\begin{split}g({{\bf r}^\prime},{\bf r},\omega ) = {g_{\text{B}}}({{\bf r}^\prime},{\bf r},\omega ) + \sum\limits_{l,m} R_{l,m}^{\text{max}}{\text{e}^{\text{i}{n_{\text{B}}}kR}}{\text{e}^{\text{i}{n_{\text{R}}}kR}}{\text{e}^{\text{i}{n_{\text{B}}}kr^{\prime}}}{\alpha _{lm}}({\bf r},\omega ) + \text{rest}\,\text{terms},\end{split}$$

where $ R_{l,m}^{\text{max}} $ contains factors from $ Y_l^m(\hat {\bf r}^{\prime}) $ and the rational functions multiplying onto the spherical Hankel and Bessel functions; they depend at most in a nonexponential manner on $ {k_{\text{B}}} $. As suggested by the superscript, the sum contains the fastest growing exponential terms that will be decisive in determining the limiting behavior of the Green function; the rest terms contain terms that grow slower or decay (exponentially) to zero when $ {k_{\text{B}}} $ is varied in the direction of the negative imaginary axis.

Next, we rewrite the expression as

(C41)$$g({{\bf r}^\prime},{\bf r},\omega ) \approx \frac{{{g_{\text{B}}}({{\bf r}^\prime},{\bf r},\omega )}}{{1 - \Xi }},$$

where

(C42)$$\Xi = \sum\limits_{l,m} R_{l,m}^{\text{max}}{\alpha _{lm}}({\bf r},\omega)\frac{{{\text{e}^{\text{i}{n_{\text{B}}}kR}}{\text{e}^{\text{i}{n_{\text{R}}}kR}}{\text{e}^{\text{i}{n_{\text{B}}}kr^{\prime}}}}}{{g({{\bf r}^\prime},{\bf r},\omega )}}.$$

The use of “$ \approx $” in Eq. (C41) signifies that we dropped the rest terms in Eq. (C40), which are unimportant in the limit of interest. From Eqs. (C36) and (C38), we can appreciate that if $ g({{\bf r}^\prime},{\bf r},\omega ) $ tends to zero as $ \omega $ is varied in the direction of the negative imaginary axis, then $ {\alpha _{lm}}({\bf r},\omega ) $ must also tend to zero in this limit—and it must do this faster than $ \exp \{ - \text{i}{n_{\text{R}}}kr^{\prime}\} $. One term in the expansion for $ g({{\bf r}^\prime},{\bf r},\omega ) $ will tend to zero slower than the other terms. Denoting this term $ l = L $ and $ m = M $, we can simplify the fraction in Eq. (C42) by $ \exp \{ \text{i}{n_{\text{R}}}kr^{\prime}\} {\alpha _{LM}} $ to find

(C43)$$\Xi = \frac{{R_{L,M}^{\text{max}}{\text{e}^{\text{i}k({n_{\text{B}}}R + {n_{\text{R}}}R + {n_{\text{B}}}r^{\prime} - {n_{\text{R}}}r^{\prime})}}}}{{1 + \text{rest}\,\text{terms}}} \\ \; + \sum_{{l \ne L \atop m\ne M}} R_{l,m}^{\text{max}}\frac{{{\alpha _{lm}}({\bf r},\omega )}}{{{\alpha _{LM}}({\bf r},\omega )}}\frac{{{\text{e}^{\text{i}k({n_{\text{B}}}R + {n_{\text{R}}}R + {n_{\text{B}}}r^{\prime} - {n_{\text{R}}}r^{\prime})}}}}{{1 + \text{rest}\,\text{terms}}},$$

where the rest terms in the numerators all tend to zero as $ {k_{\text{B}}} $ is varied in the direction of the negative imaginary axis. By construction, the limiting behavior of $ \Xi $ is now governed by the first term in Eq. (C43).

In the last step, we expand the background Green function $ {g_{\text{B}}}({{\bf r}^\prime},{\bf r},\omega ) $ in spherical Hankel and Bessel functions defined with respect to the center of the sphere using Eq. (C37) and simplify the fraction in Eq. (C41) by $ \exp \{ \text{i}k({n_{\text{B}}}R + {n_{\text{R}}}R + {n_{\text{B}}}r^{\prime} - {n_{\text{R}}}r^{\prime})\} $. In this way, we can see that the Green function $ g({{\bf r}^\prime},{\bf r}, - \text{i}\xi ) $ tends to zero in the limit $ \xi \to \infty $ if

(C44)$$r \lt R + \frac{{{n_{\text{R}}}}}{{{n_{\text{B}}}}}(R - r^{\prime}).$$

Thus, the region of consistency for the spherical resonator is itself a sphere with the same center and with a radius that depends on $ r^{\prime} $. As $ r \to R $, the radius of the region of consistency tends to the radius of the spherical resonator.

Generalization to the tensor case. For $ {\bf r} $ outside the sphere, the electric field Green tensor satisfies the Dyson equation in Eq. (A29) where the Green tensor in the homogeneous background material can be written in terms of the scalar background Green function as in Eq. (A22). As in the scalar case, even if we do not know the explicit expression, we know that we can expand all components of $ {{\bf G}^{\text{EE}}}({\bf r},{{\bf r}^\prime},\omega ) $ in terms of spherical wave functions as in Eq. (C38). In this way, we can expand all terms in the integrand to end up with nine separate equations of the form in Eq. (C40). The rest of the argumentation remains valid, so we conclude that the region of consistency is the same as in the scalar case in Eq. (C44).

Funding

Deutsche Forschungsgemeinschaft (SFB 951 HIOS B10, Project 182087777, DIP FO 703/2-1, DIP SCHM 1049/7-1); Bundesministerium für Bildung und Forschung (13N14149).

Acknowledgment

We would like to express our sincere gratitude to Dirk Hundertmark (Karlsruhe Institute of Technology), Herbert Koch (Bonn University), Stephen Hughes (Queen’s University), Felix Binkowski and Lin Zschiedrich (Zuse Institute Berlin), and Christian Wolff (University of Southern Denmark) for numerous inspiring discussions. Also, we thank Julia Werra and Benjamin Zaslansky for help with the initial BEM calculations and Thomas Kiel for Blender support . A special thanks to Jakob Rosenkrantz de Lasson, who developed and implemented the code for the VIE formulation calculations.

References

1. G. García-Calderón and R. Peierls, “Resonant states and their uses,” Nucl. Phys. A265, 443–460 (1976). [CrossRef]  

2. G. García-Calderón, “Theory of resonant states: an exact analytical approach for open quantum systems,” in Unstable States in the Continuous Spectra, Part I: Analysis, Concepts, Methods, and Results, Advances in Quantum Chemistry, C. A. Nicolaides and E. Brändas, eds. (Academic, 2010), Vol. 60, pp. 407–455.

3. S. C. Hill and R. E. Benner, “Morphology-dependent resonances,” in Optical Effects Associated with Small Particles, R. K. Chang and P. W. Barber, eds. (World Scientific, 1988), Chap. 1, pp. 1–61.

4. B. R. Johnson, “Theory of morphology-dependent resonances: shape resonances and width formulas,” J. Opt. Soc. Am. A 10, 343–352 (1993). [CrossRef]  

5. E. S.-C. Ching, P.-T. Leung, and K. Young, “The role of quasi-normal modes, in Optical Processes in Microcavities, R. K. Chang and A. J. Campillo, eds. (World Scientific, 1996), Chap. 1, pp. 1–76.

6. E. S. C. Ching, P. T. Leung, A. M. van den Brink, W. M. Suen, S. S. Tong, and K. Young, “Quasinormal-mode expansion for waves in open systems,” Rev. Mod. Phys. 70, 1545–1554 (1998). [CrossRef]  

7. G. Gamow, “Zur quantentheorie des atomkernes,” Zeitschrift für Physik 51, 204–212 (1928). [CrossRef]  

8. Y. B. Zel’dovich, “On the theory of unstable states,” Sov. Phys. JETP 12, 542 (1961).

9. L. A. Vaynshteyn, Open Resonators and Open Waveguides (The Golem, 1969).

10. H. C. van de Hulst, Light Scattering by Small Particles (Dover, 1981).

11. R. Fuchs and K. L. Kliewer, “Optical modes of vibration in an ionic crystal sphere,” J. Opt. Soc. Am. 58, 319–330 (1968). [CrossRef]  

12. C. V. Vishveshwara, “Scattering of gravitational radiation by a Schwarzschild black-hole,” Nature 227, 936 (1970). [CrossRef]  

13. V. P. Frolov and I. D. Novinov, “Black hole perturbations,” in Black Hole Physics, A. van der Merwe, ed. (Springer, 1998), Chap. 4, pp. 87–149.

14. W. B. Campbell and T. A. Morgan, “Maxwell form of the linear theory of gravitation,” Am. J. Phys 44, 356–365 (1976). [CrossRef]  

15. E. S. C. Ching, P. T. Leung, W. M. Suen, and K. Young, “Quasinormal mode expansion for linearized waves in gravitational systems,” Phys. Rev. Lett. 74, 4588–4591 (1995). [CrossRef]  

16. E. M. Purcell, “Spontaneous emission probabilities at radio frequencies,” Phys. Rev. 69, 681 (1946). [CrossRef]  

17. P. T. Kristensen, C. V. Vlack, and S. Hughes, “Generalized effective mode volume for leaky optical cavities,” Opt. Lett. 37, 1649–1651 (2012). [CrossRef]  

18. C. Sauvan, J. P. Hugonin, I. S. Maksymov, and P. Lalanne, “Theory of the spontaneous optical emission of nanosize photonic and plasmon resonators,” Phys. Rev. Lett. 110, 237401 (2013). [CrossRef]  

19. H. M. Lai, P. T. Leung, K. Young, P. W. Barber, and S. C. Hill, “Time-independent perturbation for leaking electromagnetic modes in open systems with application to resonances in microdroplets,” Phys. Rev. A 41, 5187–5198 (1990). [CrossRef]  

20. P. T. Leung, S. Y. Liu, and K. Young, “Completeness and orthogonality of quasinormal modes in leaky optical cavities,” Phys. Rev. A 49, 3057–3067 (1994). [CrossRef]  

21. P. T. Leung and K. M. Pang, “Completeness and time-independent perturbation of morphology-dependent resonances in dielectric spheres,” J. Opt. Soc. Am. B 13, 805–817 (1996). [CrossRef]  

22. P. T. Leung, S. S. Tong, and K. Young, “Two-component eigenfunction expansion for open systems described by the wave equation I: completeness of expansion,” J. Phys. A 30, 2139–2151 (1997). [CrossRef]  

23. P. T. Leung, S. Y. Liu, S. S. Tong, and K. Young, “Time-independent perturbation theory for quasinormal modes in leaky optical cavities,” Phys. Rev. A 49, 3068–3073 (1994). [CrossRef]  

24. P. T. Leung, S. S. Tong, and K. Young, “Two-component eigenfunction expansion for open systems described by the wave equation II: linear space structure,” J. Phys. A 30, 2153–2162 (1997). [CrossRef]  

25. P. T. Leung, S. Y. Liu, and K. Young, “Completeness and time-independent perturbation of the quasinormal modes of an absorptive and leaky cavity,” Phys. Rev. A 49, 3982–3989 (1994). [CrossRef]  

26. K. M. Lee, P. T. Leung, and K. M. Pang, “Dyadic formulation of morphology-dependent resonances. I. Completeness relation,” J. Opt. Soc. Am. B 16, 1409–1417 (1999). [CrossRef]  

27. K. M. Lee, P. T. Leung, and K. M. Pang, “Dyadic formulation of morphology-dependent resonances. II. Perturbation theory,” J. Opt. Soc. Am. B 16, 1418–1430 (1999). [CrossRef]  

28. S. W. Ng, P. T. Leung, and K. M. Lee, “Dyadic formulation of morphology-dependent resonances. III. Degenerate perturbation theory,” J. Opt. Soc. Am. B 19, 154–164 (2002). [CrossRef]  

29. S. Both and T. Weiss, “First-order perturbation theory for changes in the surrounding of open optical resonators,” Opt. Lett. 44, 5917–5920 (2019). [CrossRef]  

30. W. Yan, P. Lalanne, and M. Qiu, “Perturbation theory of quasinormal modes for geometrically deformed nanoresonators,” arXiv:1909.03386v1 (2019).

31. E. A. Muljarov, W. Langbein, and R. Zimmermann, “Brillouin-Wigner perturbation theory in open electromagnetic systems,” Eur. Phys. Lett. 92, 50010 (2010). [CrossRef]  

32. P. Lind, “Completeness relations and resonant state expansions,” Phys. Rev. C 47, 1903–1920 (1993). [CrossRef]  

33. E. A. Muljarov and W. Langbein, “Resonant-state expansion of dispersive open optical systems: creating gold from sand,” Phys. Rev. B 93, 075417 (2016). [CrossRef]  

34. S. Neale and E. Muljarov, “Resonant-state expansion for planar photonic-crystal structures,” Phys. Rev. B 101, 155128 (2020). [CrossRef]  

35. S. V. Lobanov, W. Langbein, and E. A. Muljarov, “Resonant-state expansion applied to three-dimensional open optical systems: complete set of static modes,” Phys. Rev. A 100, 063811 (2019). [CrossRef]  

36. A. Settimi, S. Severini, N. Mattiucci, C. Sibilia, M. Centini, G. D’Aguanno, M. Bertolotti, M. Scalora, M. Bloemer, and C. M. Bowden, “Quasinormal-mode description of waves in one-dimensional photonic crystals,” Phys. Rev. E 68, 026614 (2003). [CrossRef]  

37. A. Settimi, S. Severini, and B. J. Hoenders, “Quasi-normal-modes description of transmission properties for photonic bandgap structures,” J. Opt. Soc. Am. B 26, 876–891 (2009). [CrossRef]  

38. M. B. Doost, W. Langbein, and E. A. Muljarov, “Resonant-state expansion applied to planar open optical systems,” Phys. Rev. A 85, 023835 (2012). [CrossRef]  

39. M. B. Doost, W. Langbein, and E. A. Muljarov, “Resonant state expansion applied to two-dimensional open optical systems,” Phys. Rev. A 87, 043827 (2013). [CrossRef]  

40. M. B. Doost, W. Langbein, and E. A. Muljarov, “Resonant-state expansion applied to three-dimensional open optical systems,” Phys. Rev. A 90, 013834 (2014). [CrossRef]  

41. P. T. Leung, W. M. Suen, C. P. Sun, and K. Young, “Waves in open systems via a biorthogonal basis,” Phys. Rev. E 57, 6101–6104 (1998). [CrossRef]  

42. Z. Li, G. Yi-Bo, and W. Cheng, “Green function and perturbation method for dissipative systems based on biorthogonal basis,” Commun. Theor. Phys. 51, 1017–1022 (2009). [CrossRef]  

43. M. Mansuripur, M. Kolesik, and P. Jakobsen, “Leaky modes of solid dielectric spheres,” Phys. Rev. A 96, 013846 (2017). [CrossRef]  

44. L. J. Armitage, M. B. Doost, W. Langbein, and E. A. Muljarov, “Resonant-state expansion applied to planar waveguides,” Phys. Rev. A 89, 053832 (2014). [CrossRef]  

45. S. V. Lobanov, W. Langbein, and E. A. Muljarov, “Resonant-state expansion of three-dimensional open optical systems: light scattering,” Phys. Rev. A 98, 033820 (2018). [CrossRef]  

46. J. R. de Lasson, J. Mørk, and P. T. Kristensen, “Three-dimensional integral equation approach to light scattering, extinction cross sections, local density of states, and quasi-normal modes,” J. Opt. Soc. Am. B 30, 1996–2007 (2013). [CrossRef]  

47. J. Wiersig, “Boundary element method for resonances in dielectric microcavities,” J. Opt. A 5, 53–60 (2002). [CrossRef]  

48. J. Mäkitalo, M. Kauranen, and S. Suuriniemi, “Modes and resonances of plasmonic scatterers,” Phys. Rev. B 89, 165429 (2014). [CrossRef]  

49. F. Alpeggiani, S. D. Agostino, D. Sanvitto, and D. Gerace, “Visible quantum plasmonics from metallic nanodimers,” Sci. Rep. 6, 34772 (2016). [CrossRef]  

50. T. Weiss, N. A. Gippius, S. G. Tikhodeev, G. Granet, and H. Giessen, “Derivation of plasmonic resonances in the Fourier modal method with adaptive spatial resolution and matched coordinates,” J. Opt. Soc. Am. A 28, 238–244 (2011). [CrossRef]  

51. T. Weiss, M. Mesch, M. Schäferling, H. Giessen, W. Langbein, and E. A. Muljarov, “From dark to bright: first-order perturbation theory with analytical mode normalization for plasmonic nanoantenna arrays applied to refractive index sensing,” Phys. Rev. Lett. 116, 237401 (2016). [CrossRef]  

52. T. Weiss, M. Schäferling, H. Giessen, N. A. Gippius, S. G. Tikhodeev, W. Langbein, and E. A. Muljarov, “Analytical normalization of resonant states in photonic crystal slabs and periodic arrays of nanoantennas at oblique incidence,” Phys. Rev. B 96, 045129 (2017). [CrossRef]  

53. C. Sauvan, J. P. Hugonin, R. Carminati, and P. Lalanne, “Modal representation of spatial coherence in dissipative and resonant photonic systems,” Phys. Rev. A 89, 043825 (2014). [CrossRef]  

54. Z. Hu and Y. Y. Lu, “Efficient analysis of photonic crystal devices by Dirichlet-to-Neumann maps,” Opt. Express 16, 17383–17399 (2008). [CrossRef]  

55. J. R. de Lasson, P. T. Kristensen, J. Mørk, and N. Gregersen, “Roundtrip matrix method for calculating the leaky resonant modes of open nanophotonic structures,” J. Opt. Soc. Am. A 31, 2142–2151 (2014). [CrossRef]  

56. P. T. Kristensen, J. R. de Lasson, and N. Gregersen, “Calculation, normalization, and perturbation of quasinormal modes in coupled cavity-waveguide systems,” Opt. Lett. 39, 6359–6362 (2014). [CrossRef]  

57. F. Römer and B. Witzigmann, “Spectral and spatial properties of the spontaneous emission enhancement in photonic crystal cavities,” J. Opt. Soc. Am. B 25, 31–39 (2008). [CrossRef]  

58. B. Vial, F. Zolla, A. Nicolet, and M. Commandré, “Quasimodal expansion of electromagnetic fields in open two-dimensional structures,” Phys. Rev. A 89, 023829 (2014). [CrossRef]  

59. W. Yan, R. Faggiani, and P. Lalanne, “Rigorous modal analysis of plasmonic nanoresonators,” Phys. Rev. B 97, 205422 (2018). [CrossRef]  

60. E. A. Muljarov and W. Langbein, “Exact mode volume and Purcell factor of open optical systems,” Phys. Rev. B 94, 235438 (2016). [CrossRef]  

61. P. T. Kristensen, R.-C. Ge, and S. Hughes, “Normalization of quasinormal modes in leaky optical cavities and plasmonic resonators,” Phys. Rev. A 92, 053810 (2015). [CrossRef]  

62. E. A. Muljarov and W. Langbein, “Comment on ‘normalization of quasinormal modes in leaky optical cavities and plasmonic resonators’,” Phys. Rev. A 96, 017801 (2017). [CrossRef]  

63. P. T. Kristensen, R.-C. Ge, and S. Hughes, “Reply to ‘comment on “normalization of quasinormal modes in leaky optical cavities and plasmonic resonators’”,” Phys. Rev. A 96, 017802 (2017). [CrossRef]  

64. B. Stout, R. Colom, N. Bonod, and R. McPhedran, “Eigenstate normalization for open and dispersive systems,” arXiv:1903.07183v1 (2019).

65. Q. Bai, M. Perrin, C. Sauvan, J.-P. Hugonin, and P. Lalanne, “Efficient and intuitive method for the analysis of light scattering by a resonant nanostructure,” Opt. Express 21, 27371–27382 (2013). [CrossRef]  

66. M. Perrin, “Eigen-energy effects and non-orthogonality in the quasi-normal mode expansion of Maxwell equations,” Opt. Express 24, 27137–27151 (2016). [CrossRef]  

67. L. Zschiedrich, F. Binkowski, N. Nikolay, O. Benson, G. Kewes, and S. Burger, “Riesz-projection-based theory of light-matter interaction in dispersive nanoresonators,” Phys. Rev. A 98, 043806 (2018). [CrossRef]  

68. F. Binkowski, L. Zschiedrich, and S. Burger, “A Riesz-projection-based method for nonlinear eigenvalue problems,” arXiv:1811.11624v2 (2018).

69. F. Binkowski, L. Zschiedrich, M. Hammerschmidt, and S. Burger, “Modal analysis for nanoplasmonics with nonlocal material properties,” Phys. Rev. B 100, 155406 (2019). [CrossRef]  

70. R.-C. Ge, P. T. Kristensen, J. F. Young, and S. Hughes, “Quasinormal mode approach to modelling light-emission and propagation in nanoplasmonics,” New J. Phys. 16, 113048 (2014). [CrossRef]  

71. R.-C. Ge and S. Hughes, “Design of an efficient single photon source from a metallic nanorod dimer: a quasi-normal mode finite-difference time-domain approach,” Opt. Lett. 39, 4235–4238 (2014). [CrossRef]  

72. M. K. Dezfouli and S. Hughes, “Regularized quasinormal modes for plasmonic resonators and open cavities,” Phys. Rev. B 97, 115302 (2018). [CrossRef]  

73. M. I. Abdelrahman and B. Gralak, “Completeness and divergence-free behavior of the quasi-normal modes using causality principle,” OSA Continuum 1, 340–348 (2018). [CrossRef]  

74. R. Colom, R. McPhedran, B. Stout, and N. Bonod, “Modal expansion of the scattered field: causality, nondivergence, and nonresonant contribution,” Phys. Rev. B 98, 085418 (2018). [CrossRef]  

75. P. T. Kristensen, J. R. de Lasson, M. Heuck, N. Gregersen, and J. Mørk, “On the theory of coupled modes in optical cavity-waveguide structures,” J. Lightwave Technol. 35, 4247–4259 (2017). [CrossRef]  

76. A. B. Akimov, N. A. Gippius, and S. G. Tikhodeev, “Optical fano resonances in photonic crystal slabs near diffraction threshold anomalies,” JETP Lett. 93, 427–430 (2011). [CrossRef]  

77. V. Grigoriev, A. Tahri, S. Varault, B. Rolly, B. Stout, J. Wenger, and N. Bonod, “Optimization of resonant effects in nanostructures via Weierstrass factorization,” Phys. Rev. A 88, 011803 (2013). [CrossRef]  

78. F. Alpeggiani, N. Parappurath, E. Verhagen, and L. Kuipers, “Quasinormal-mode expansion of the scattering matrix,” Phys. Rev. X 7, 021035 (2017). [CrossRef]  

79. T. Weiss and E. A. Muljarov, “How to calculate the pole expansion of the optical scattering matrix from the resonant states,” Phys. Rev. B 98, 085433 (2018). [CrossRef]  

80. G. Unger, A. Trügler, and U. Hohenester, “Novel modal approximation scheme for plasmonic transmission problems,” Phys. Rev. Lett. 121, 246802 (2018). [CrossRef]  

81. P. T. Kristensen, M. Heuck, and J. Mørk, “Optimal switching using coherent control,” Appl. Phys. Lett. 102, 041107 (2013). [CrossRef]  

82. C. Gigli, T. Wu, G. Marino, A. Borne, G. Leo, and P. Lalanne, “Quasinormal-mode modeling and design in nonlinear nano-optics,” arXiv:1911.06373v1 (2019).

83. J. Yang, H. Giessen, and P. Lalanne, “Simple analytical expression for the peak-frequency shifts of plasmonic resonances for sensing,” Nano Lett. 15, 3439–3444 (2015). [CrossRef]  

84. M. K. Dezfouli, C. Tserkezis, N. A. Mortensen, and S. Hughes, “Nonlocal quasinormal modes for arbitrarily shaped three-dimensional plasmonic resonators,” Optica 4, 1503–1509 (2017). [CrossRef]  

85. J. R. de Lasson, P. T. Kristensen, J. Mørk, and N. Gregersen, “Semianalytical quasi-normal mode theory for the local density of states in coupled photonic crystal cavity–waveguide structures,” Opt. Lett. 40, 5790–5793 (2015). [CrossRef]  

86. T. Malhotra, R.-C. Ge, M. K. Dezfouli, A. Badolato, N. Vamivakas, and S. Hughes, “Quasinormal mode theory and design of on-chip single photon emitters in photonic crystal coupled-cavity waveguides,” Opt. Express 24, 13574–13583 (2016). [CrossRef]  

87. M. K. Dezfouli, R. Gordon, and S. Hughes, “Modal theory of modified spontaneous emission of a quantum emitter in a hybrid plasmonic photonic-crystal cavity system,” Phys. Rev. A 95, 013846 (2017). [CrossRef]  

88. R.-C. Ge and S. Hughes, “Quasinormal mode theory and modelling of electron energy loss spectroscopy for plasmonic nanostructures,” J. Opt. 18, 054002 (2016). [CrossRef]  

89. A. Hörl, A. Trügler, and U. Hohenester, “Full three-dimensonal reconstruction of the dyadic green tensor from electron energy loss spectroscopy of plasmonic nanoparticles,” ACS Photon. 2, 1429–1435 (2015). [CrossRef]  

90. J. Andreasen, A. A. Asatryan, L. C. Botten, M. A. Byrne, H. Cao, L. Ge, L. Labonté, P. Sebbah, A. D. Stone, H. E. Türeci, and C. Vanneste, “Modes of random lasers,” Adv. Opt. Photon. 3, 88–127 (2011). [CrossRef]  

91. H. E. Türeci, A. D. Stone, and B. Collier, “Self-consistent multimode lasing theory for complex or random lasing media,” Phys. Rev. A 74, 043822 (2006). [CrossRef]  

92. W. Cartar, J. Mørk, and S. Hughes, “Self-consistent Maxwell-Bloch model of quantum-dot photonic-crystal-cavity lasers,” Phys. Rev. A 96, 023859 (2017). [CrossRef]  

93. P. T. Kristensen, J. E. Mortensen, P. Lodahl, and S. Stobbe, “Shell theorem for spontaneous emission,” Phys. Rev. B 88, 205308 (2013). [CrossRef]  

94. X. Zambrana-Puyalto and N. Bonod, “Purcell factor of spherical Mie resonators,” Phys. Rev. B 91, 195422 (2015). [CrossRef]  

95. C. Carlson and S. Hughes, “Dissipative modes, Purcell factors and directional beta factors in gold bowtie nanoantenna structures,” arXiv:1910.10110v1 (2019).

96. R.-C. Ge and S. Hughes, “Quantum dynamics of two quantum dots coupled through localized plasmons: an intuitive and accurate quantum optics approach using quasinormal modes,” Phys. Rev. B 92, 205420 (2015). [CrossRef]  

97. J. Mørk, P. T. Kristensen, P. Kaer, M. Heuck, Y. Yu, and N. Gregersen, “Cavity photonics,” in Photonics: Scientific Foundations, Technology and Applications, D. L. Andrews, ed. (Wiley, 2015), Vol. II, Chap. 2, pp. 21–51.

98. M. K. Dezfouli, R. Gordon, and S. Hughes, “Molecular optomechanics in the strong coupling regime using hybrid metal-dielectric cavity modes,” ACS Photon. 6, 1400–1408 (2019). [CrossRef]  

99. K. C. Ho, P. T. Leung, A. Maassen van den Brink, and K. Young, “Second quantization of open systems using quasinormal modes,” Phys. Rev. E 58, 2965–2978 (1998). [CrossRef]  

100. S. M. Dutra and G. Nienhuis, “Quantized mode of a leaky cavity,” Phys. Rev. A 62, 063805 (2000). [CrossRef]  

101. S. Severini, A. Settimi, C. Sibilia, M. Bertolotti, A. Napoli, and A. Messina, “Second quantization and atomic spontaneous emission inside one-dimensional photonic crystals via a quasinormal-modes approach,” Phys. Rev. E 70, 056614 (2004). [CrossRef]  

102. S. Franke, S. Hughes, M. K. Dezfouli, P. T. Kristensen, K. Busch, A. Knorr, and M. Richter, “Quantization of quasinormal modes for open cavities and plasmonic cavity quantum electrodynamics,” Phys. Rev. Lett. 122, 213901 (2019). [CrossRef]  

103. R. Sammut and A. W. Snyder, “Leaky modes on a dielectric waveguide: orthogonality and excitation,” Appl. Opt. 15, 1040–1044 (1976). [CrossRef]  

104. S. V. Lobanov, G. Zoriniants, W. Langbein, and E. A. Muljarov, “Resonant-state expansion of light propagation in nonuniform waveguides,” Phys. Rev. A 95, 053848 (2017). [CrossRef]  

105. S. Upendar, I. Allayarov, M. A. Schmidt, and T. Weiss, “Analytical mode normalization and resonant state expansion for bound and leaky modes in optical fibers—an efficient tool to model transverse disorder,” Opt. Express 26, 22536–22546 (2018). [CrossRef]  

106. I. Allayarov, S. Upendar, M. A. Schmidt, and T. Weiss, “Analytic mode normalization for the Kerr nonlinearity parameter: prediction of nonlinear gain for leaky modes,” Phys. Rev. Lett. 121, 213905 (2018). [CrossRef]  

107. N. Moiseyev, Non-Hermitian Quantum Mehcanics (Cambridge Universtity, 2011).

108. R. A. Konoplya and A. Zhidenko, “Quasinormal modes of black holes: from astrophysics to string theory,” Rev. Mod. Phys. 83, 793–836 (2011). [CrossRef]  

109. J. Kergomard, V. Debut, and D. Matignon, “Resonance modes in a one-dimensional medium with two purely resistive boundaries: calculation methods, orthogonality, and completeness,” J. Acoust. Soc. Am. 119, 1356–1367 (2006). [CrossRef]  

110. P. T. Kristensen and S. Hughes, “Modes and mode volumes of leaky optical cavities and plasmonic nanoresonators,” ACS Photon. 1, 2–10 (2014). [CrossRef]  

111. P. Lalanne, W. Yan, K. Vynck, C. Sauvan, and J.-P. Hugonin, “Light interaction with photonic and plasmonic resonances,” Laser Photon. Rev. 12, 1700113 (2018). [CrossRef]  

112. R. H. Cole, Theory of Ordinary Differential Equations (Appleton-Century-Crofts, 1968).

113. L. Botten, M. Craig, R. McPhedran, J. Adams, and J. Andrewartha, “The finitely conducting lamellar diffraction grating,” Opt. Acta 28, 1087–1102 (1981). [CrossRef]  

114. C. Wolff, K. Busch, and N. A. Mortensen, “Modal expansions in periodic photonic systems with material loss and dispersion,” Phys. Rev. B 97, 104203 (2018). [CrossRef]  

115. S. G. Tikhodeev, A. L. Yablonskii, E. A. Muljarov, N. A. Gippius, and T. Ishihara, “Quasiguided modes and optical properties of photonic crystal slabs,” Phys. Rev. B 66, 045102 (2002). [CrossRef]  

116. G. Demésy, A. Nicolet, B. Gralak, C. Geuzaine, C. Campos, and J. E. Roman, “Non-linear eigenvalue problems with GetDP and SLEPc: eigenmode computations of frequency-dispersive photonic open structures,” arXiv:1802.02363v3 (2018).

117. P. T. Kristensen, MATLAB Code for Dielectric Barrier Calculations, figshare (2020), https://doi.org/10.6084/m9.figshare.9810383.

118. P. T. Kristensen, MATLAB Code for Plasmonic Dimer Calculations Using MNPBEM, figshare (2020), https://doi.org/10.6084/m9.figshare.9810386.

119. U. Hohenester and A. Trügler, “MNPBEM—a Matlab toolbox for the simulation of plasmonic nanoparticles,” Comput. Phys. Commun. 183, 370–381 (2012). [CrossRef]  

120. J. Waxenegger, A. Trugler, and U. Hohenester, “Plasmonics simulations with the MNPBEM toolbox: consideration of substrates and layer structures,” Comput. Phys. Commun. 193, 138–150 (2015). [CrossRef]  

121. U. Hohenester, “Making simulations with the MNPBEM toolbox big: hierarchical matrices and iterative solvers,” Comput. Phys. Commun. 222, 209–228 (2018). [CrossRef]  

122. B. Maes, J. Petráček, S. Burger, P. Kwiecien, J. Luksch, and I. Richter, “Simulations of high-Q optical nanocavities with a gradual 1D bandgap,” Opt. Express 21, 6794–6806 (2013). [CrossRef]  

123. J. R. de Lasson, L. H. Frandsen, P. Gutsche, S. Burger, O. S. Kim, O. Breinbjerg, A. Ivinskaya, F. Wang, O. Sigmund, T. Häyrynen, A. V. Lavrinenko, J. Mørk, and N. Gregersen, “Benchmarking five numerical simulation techniques for computing resonance wavelengths and quality factors in photonic crystal membrane line defect cavities,” Opt. Express 26, 11366–11392 (2018). [CrossRef]  

124. P. Lalanne, W. Yan, A. Gras, C. Sauvan, J.-P. Hugonin, M. Besbes, G. Demésy, M. D. Truong, B. Gralak, F. Zolla, A. Nicolet, F. Binkowski, L. Zschiedrich, S. Burger, J. Zimmerling, R. Remis, P. Urbach, H. T. Liu, and T. Weiss, “Quasinormal mode solvers for resonators with dispersive materials,” J. Opt. Soc. Am. A 36, 686–704 (2019). [CrossRef]  

125. P. Lalanne, C. Sauvan, and J. Hugonin, “Photon confinement in photonic crystal nanocavities,” Laser Photon. Rev. 2, 514–526 (2008). [CrossRef]  

126. T. H. Taminiau, F. D. Stefani, and N. F. van Hulst, “Optical nanorod antennas modeled as cavities for dipolar emitters: evolution of sub- and super-radiant modes,” Nano Lett. 11, 1020 (2011). [CrossRef]  

127. F. Yang, H. Liu, H. Jia, and Y. Zhong, “Analytical description of quasi-normal mode in resonant plasmonic nano cavities,” J. Opt. 18, 035003 (2016). [CrossRef]  

128. M. Oxborrow, “Ex-house 2D finite-element simulation of the whispering-gallery modes of arbitrarily shaped axisymmetric electromagnetic resonators,” arXiv:quant-ph/0607156v2 (2006).

129. F. J. García de Abajo and A. Howie, “Retarded field calculation of electron energy loss in inhomogeneous dielectrics,” Phys. Rev. B 65, 115418 (2002). [CrossRef]  

130. E. A. Muljarov and T. Weiss, “Resonant-state expansion for open optical systems: generalization to magnetic, chiral, and bi-anisotropic materials,” Opt. Lett. 43, 1978–1981 (2018). [CrossRef]  

131. S. Silver, “Radiation from current distributions,” in Microwave Antenna Theory and Design, S. Silver, ed. (McGraw-Hill, 1949), Chap. 3, pp. 61–106.

132. P. Martin, Multiple Scattering. Interaction of Time-Harmonic Waves with N Obstacles (Cambridge University, 2006).

133. R. G. Barrera, G. A. Estévez, and J. Giraldo, “Vector spherical harmonics and their application to magnetostatics,” Eur. J. Phys. 6, 287–294 (1985). [CrossRef]  

134. M. Abromowitz, E. Irene, and A. Stegun, Handbook of Mathematical Functions with Formulas, Graphs, and Mathematical Tables, 10th printing (U.S. Government Printing Office, 1972).

135. J. Wiersig, “Structure of whispering-gallery modes in optical microdisks perturbed by nanoparticles,” Phys. Rev. A 84, 063828 (2011). [CrossRef]  

136. A. Raman and S. Fan, “Photonic band structure of dispersive metamaterials formulated as a Hermitian eigenvalue problem,” Phys. Rev. Lett. 104, 087401 (2010). [CrossRef]  

137. L. Trefethen and J. Weideman, “The exponentially convergent trapezoidal rule,” SIAM Rev. 56, 385–458 (2014). [CrossRef]  

138. P. T. Kristensen, P. Lodahl, and J. Mørk, “Light propagation in finite-sized photonic crystals: multiple scattering using an electric field integral equation,” J. Opt. Soc. Am. B 27, 228–237 (2010). [CrossRef]  

139. J.-M. Jin, Theory and Computation of Electromagnetic Fields (Wiley, 2010).

140. N. A. Gippius, T. Weiss, S. G. Tikhodeev, and H. Giessen, “Resonant mode coupling of optical resonances in stacked nanostructures,” Opt. Express 18, 7569–7574 (2010). [CrossRef]  

141. T. Weiss, N. A. Gippius, G. Granet, S. G. Tikhodeev, R. Taubert, L. Fu, H. Schweizer, and H. Giessen, “Strong resonant mode coupling of Fabry–Perot and grating resonances in stacked two-layer systems,” Photon. Nanostruct. Fundam. Appl. 9, 390–397 (2011). [CrossRef]  

142. H. A. Haus, Waves and Fields in Optoelectronics (Prentice Hall, 1984).

143. W. Suh, Z. Wang, and S. Fan, “Temporal coupled-mode theory and the presence of non-orthogonal modes in lossless multimode cavities,” IEEE J. Quantum Electron. 40, 1511–1518 (2004). [CrossRef]  

144. J. D. Joannopoulos, S. G. Johnson, J. N. Winn, and R. D. Meade, Photonic Crystals: Molding the Flow of Light, 2nd ed. (Princeton University, 2008).

145. G. Mur, “Absorbing boundary conditions for the finite-difference approximation of the time-domain electromagnetic-field equations,” IEEE Trans. Electromagn. Compat. EMC-23, 377–382 (1981). [CrossRef]  

146. A. Tavlove and S. C. Hagness, Computational Electromagnetics: The Finite-Difference Time-Domain Method, 3rd ed. (Artech House, 2003).

147. K. Busch, M. König, and J. Niegemann, “Discontinuous Galerkin methods in nanophotonics,” Laser Photon. Rev. 5, 773–809 (2011). [CrossRef]  

148. S. Fiedler, S. Raza, R. Ai, J. Wang, K. Busch, N. Stenger, N. A. Mortensen, and C. Wolff, “Importance of substrates for the visibility of 'dark' plasmonic modes,” Opt. Express 28, 13938–13948 (2020). [CrossRef]  

149. M. Moeferdt, T. Kiel, T. Sproll, F. Intravaia, and K. Busch, “Plasmonic modes in nanowire dimers: a study based on the hydrodynamic Drude model including nonlocal and nonlinear effects,” Phys. Rev. B 97, 075431 (2018). [CrossRef]  

150. V. P. Gupta, Principles and Applications of Quantum Chemistry (Academic, 2016).

151. L. E. Ballentine, Quantum Mechanics: a Modern Development (World Scientific, 1998).

152. S. G. Johnson, M. Ibanescu, M. A. Skorobogatiy, O. Weisberg, J. D. Joannopoulos, and Y. Fink, “Perturbation theory for Maxwell’s equations with shifting material boundaries,” Phys. Rev. E 65, 066611 (2002). [CrossRef]  

153. H. T. Dung, L. Knöll, and D.-G. Welsch, “Spontaneous decay in the presence of dispersing and absorbing bodies: general theory and application to a spherical cavity,” Phys. Rev. A 62, 053804 (2000). [CrossRef]  

154. L. Novotny and B. Hecht, Principles of Nano-Optics (Cambridge University, 2007).

155. K. Drexhage, “Influence of a dielectric interface on fluorescence decay time,” J. Lumin. 1–2, 693–701 (1970). [CrossRef]  

156. P. Goy, J. M. Raimond, M. Gross, and S. Haroche, “Observation of cavity-enhanced single-atom spontaneous emission,” Phys. Rev. Lett. 50, 1903–1906 (1983). [CrossRef]  

157. J. M. Gérard, B. Sermage, B. Gayral, B. Legrand, E. Costard, and V. Thierry-Mieg, “Enhanced spontaneous emission by quantum boxes in a monolithic optical microcavity,” Phys. Rev. Lett. 81, 1110–1113 (1998). [CrossRef]  

158. P. de Vries and A. Lagendijk, “Resonant scattering and spontaneous emission in dielectrics: microscopic derivation of local-field effects,” Phys. Rev. Lett. 81, 1381–1384 (1998). [CrossRef]  

159. M. Fleischhauer, “Spontaneous emission and level shifts in absorbing disordered dielectrics and dense atomic gases: a Green’s-function approach,” Phys. Rev. A 60, 2534–2539 (1999). [CrossRef]  

160. M. S. Tomaš, “Local-field corrections to the decay rate of excited molecules in absorbing cavities: the Onsager model,” Phys. Rev. A 63, 053811 (2001). [CrossRef]  

161. H. T. Dung, S. Y. Buhmann, and D.-G. Welsch, “Local-field correction to the spontaneous decay rate of atoms embedded in bodies of finite size,” Phys. Rev. A 74, 023803 (2006). [CrossRef]  

162. F. Intravaia and K. Busch, “Fluorescence in nonlocal dissipative periodic structures,” Phys. Rev. A 91, 053836 (2015). [CrossRef]  

163. P. M. Morse and H. Feshbach, Methods of Theoretical Physics (McGraw-Hill, 1953).

164. G. B. Arfken and H. J. Weber, Mathematical Methods for Physicists, 5th ed. (Academic, 2001).

165. C.-T. Tai, Dyadic Green Functions in Electromagnetic Theory, 2nd ed. (IEEE, 1994).

166. B. A. Lippmann and J. Schwinger, “Variational principles for scattering processes. I,” Phys. Rev. 79, 469–480 (1950). [CrossRef]  

167. F. J. Dyson, “The S matrix in quantum electrodynamics,” Phys. Rev. 75, 1736–1755 (1949)

168. J. Van Bladel, “Some remarks on Green’s dyadic for infinite space,” IRE Trans. Antennas Propag. 9, 563–566 (1961). [CrossRef]  

169. A. D. Yaghjian, “Electric dyadic Green’s functions in the source region,” Proc. IEEE 68, 248–263 (1980). [CrossRef]  

170. P. de Vries, D. V. van Coevorden, and A. Lagendijk, “Point scatterers for classical waves,” Rev. Mod. Phys. 70, 447–466 (1998). [CrossRef]  

171. M. Wubs, L. G. Suttorp, and A. Lagendijk, “Multiple-scattering approach to interatomic interactions and superradiance in inhomogeneous dielectrics,” Phys. Rev. A 70, 053823 (2004). [CrossRef]  

172. O. J. F. Martin and N. B. Piller, “Electromagnetic scattering in polarizable backgrounds,” Phys. Rev. E 58, 3909–3915 (1998). [CrossRef]  

173. M. Kreiter, S. Mittler, W. Knoll, and J. R. Sambles, “Surface plasmon-related resonances on deep and asymmetric gold gratings,” Phys. Rev. B 65, 125415 (2002). [CrossRef]  

174. N. Grady, N. Halas, and P. Nordlander, “Influence of dielectric function properties on the optical response of plasmon resonant metallic nanoparticles,” Chem. Phys. Lett. 399, 167–171 (2004). [CrossRef]  

175. P. B. Johnson and R. W. Christy, “Optical constants of the noble metals,” Phys. Rev. B 6, 4370–4379 (1972).

Philip Trøst Kristensen received the degree of civilingeniør in 2006 and a Ph.D. in nanophotonics in 2010 from DTU Fotonik at the Technical University of Denmark. Subsequently, he continued as a postdoctoral researcher at DTU Fotonik with a number of research stays at Queen’s University in Kingston Ontario. In 2014, he moved to the Humboldt-Universität zu Berlin as a postdoctoral researcher enabled by a Carlsberg Foundation internationalization fellowship. In 2020 he returned to DTU Fotonik to take up a position as a senior researcher. His research has so far been concerned with theoretical and computational aspects of light scattering and light–matter interaction in micro- and nano-structured materials.

Kathrin Herrmann received her master’s degree in physics from the Humboldt-Universität zu Berlin in 2013, where she subsequently worked in the theoretical optics and photonics group of Prof. Busch until 2018. Her research interests include the modal description of leaky resonators.

Francesco Intravaia received his laurea in theoretical physics in 2002 from the Università degli Studi di Palermo, Italy. In 2005, he was awarded with a doctoral degree in physics from the Université Paris VI for his work at the Laboratoire Kastler-Brossel in Paris, France. Later, he was an Alexander von Humboldt Postdoctoral Fellow at the Universität Potsdam, Germany, and a Director’s Postdoctoral Fellow at the Los Alamos National Laboratory, New Mexico, USA. From 2014–2019, he worked at the Max Born Institute in Berlin and subsequently at the Humboldt-Universität zu Berlin, where he became a permanent staff member in 2019. His research interests focus on the theoretical investigation of light–matter interaction in micro- and mesoscopic systems, with special attention to quantum fluctuation-induced phenomena.

Kurt Busch received the Diplom and Dr. rer. nat degree in physics from the Universität Karlsruhe, Karlsruhe, Germany, in 1993 and 1996, respectively. After a postdoctoral stay at the University of Toronto, in 2000, he became a Junior Research Group Leader within the Emmy-Noether Program, Deutsche Forschungsgemeinschaft (DFG), Universität Karlsruhe. In 2004, he joined the University of Central Florida as an Associate Professor. From 2005–2011, he has been a Professor of Theoretical Physics at the Universität Karlsruhe (now Karlsruhe Institute of Technology). Since 2011, he has been Professor of Theoretical Optics & Photonics at the Humboldt-Universität zu Berlin, and he is a group leader at the Max Born Institute, Berlin, Germany. His research interests lie in the theory of light propagation and light–matter interaction in strongly scattering and complex systems. Dr. Busch is a fellow of The Optical Society and the recipient of the 2006 Carl Zeiss Research Award.