1. Introduction

1.1 Modes in conservative and non-conservative systems

Modes are central in physics, since the most general evolution of a system is a superposition of its modes, or eigenstates. From textbooks, everyone knows the special case of conservative systems studied in the mathematical framework of Hilbert spaces, e.g., particles in a potential with infinite barriers in quantum mechanics, cavities with perfectly-conducting metallic walls in electromagnetics, non-dissipative vibrating ropes with a finite length, or photonic crystals made with a perfectly periodic arrangement of non-absorbing materials.

Conservative systems have simple imposed boundary conditions (Dirichlet, Neuman, or Born-von Karman conditions) and the operators describing them are Hermitian (or self-adjoint). The eigenstates of conservative systems are usually called normal modes. They form a complete set, while the corresponding eigenvalues (frequencies or energies) are real. Moreover, it is easy to define an inner product for normal modes, namely the usual scalar product of the corresponding Hilbert space. Many beautiful consequences immediately follow, such as the mode orthogonality or the well-developed perturbation theory for linear problems. For a didactical introduction to Hermitian or conservative electromagnetic systems, we refer the interested reader to Chapters 7 and 8 in [1] and Chapter 2 in [2].

Conservative systems are obviously idealized systems. In reality, there exists flows (of energy, particle, information…) into external degrees of freedom that are outside the system subspace under consideration. In electromagnetism, there is always a given amount of loss, either by absorption or by radiative leakage in an open space. Such real open systems are referred to as non-Hermitian or non-conservative and are described by non-self-adjoint operators (non-Hermitian matrices).

The properties of non-conservative systems are much richer than those of conservative ones. The other side of the coin is that the mathematical tools used for their conceptualization are more difficult [3]. The nice properties of Hilbert spaces are lost and the eigenstates are no longer normal modes: eigenvalues are now complex and sometimes unstable with respect to small perturbations of the system [4], the modes form a complete set only under certain conditions, and the usual inner product of Hermitian systems is unfortunately invalid. The latter cannot be defined because the eigenstate field exponentially diverges at long distance. Therefore, non-Hermitian systems require more mathematical rigor and heuristics should be considered with care. However, it is still possible to prove that the eigenstates are orthonormal with a suitable regularization of the diverging fields and a redefinition of the inner product. Defining a correct inner product that allows to recover the property of orthonormality is extremely helpful for extending the well-known tools of conventional Hermitian electrodynamics to non-Hermitian systems.

Note that the notion of non-Hermitian systems spreads in physics well beyond the present electromagnetic context, e.g., in mechanics, gravitation, acoustics, or atomic, nuclear, and molecular physics. The field of non-Hermitian physics has garnered widespread interest in recent years in relation with prominent topics such as unidirectional invisibility, enhanced sensitivity, topological energy transfer, coherent perfect absorption, single-mode lasing, and robust biological transport, see [5] for a recent and extensive review on this broad topic.

1.2 Quasinormal modes, resonant states, leaky modes…

The interest in the eigenstates of non-conservative systems probably started in quantum mechanics with the calculation of the alpha decay [6] or the cross-section for nuclear reactions [7]. However, since they are the solutions of a source-free wave equation, these states can be encountered in any field of wave physics. Over the years, different communities have given them various names, each one putting the emphasis on a different characteristic.

Quantum mechanics chose to underline their finite lifetime (unstable states [8], decaying states [9]), their link with the system resonances (resonant states [10]), or the fact that they are the natural modes of the scatterer [11]. The term “quasinormal mode” was likely coined in the context of black holes and gravitational waves [4,12–14] before being used for optical cavities [15,16]. In general relativity, the scattering potential has an infinite spatial extent and not a compact support. As a result, the modes do not form a complete set; that was the main reason to call them “quasi”. In electromagnetism, the scattering potential (e.g., a permittivity change $\mathrm{\Delta }\varepsilon $) has usually a compact support and instability problems are rarely encountered. In this context, the prefix “quasi” was chosen to highlight the differences but also the similarities between Hermitian systems, which support normal modes, and non-Hermitian systems, which support modes that can be seen as almost normal under a suitable regularization.

The eigenmodes of open systems are also known as leaky modes in optical waveguides theory [17] and in acoustics [18]. Finally, let us mention the name “morphology-dependent resonance” that is sometimes given to the optical modes of microparticles [19,20].

Hereafter, we refer to the eigenmodes of non-Hermitian open systems as quasinormal modes, abbreviated in QNMs.

1.3 Scope of the review

This review focuses on the normalization, the orthogonality, and the completeness of electromagnetic QNMs. For a review of other theoretical issues (e.g., field reconstruction) or of the numerous and various applications of QNMs in photonics, we refer the interested reader to [3,21–23].

Let us recall that completeness and orthonormality are essential properties in the construction of many theories. First, completeness underpins everything: the eigenstates ${\boldsymbol{\Psi}_m}$ constitute a relevant set to decompose scattered waves only if they form a complete set. One usually writes that any scattered wave ${\boldsymbol{\Psi}_s}$ can be written as a superposition of modes, ${\boldsymbol{\Psi}_s} = {\mathrm{\Sigma }_m}{\alpha _m}{\boldsymbol{\Psi}_m}$, where ${\alpha _m}$ is the amplitude (or excitation coefficient) of the $m$th mode. In classical textbooks on conservative systems, normalization is directly related to orthogonality through the concept of inner product. With the standard bra-ket notation, $\langle{\boldsymbol{\Psi}_m}\textrm{|}{\boldsymbol{\Psi}_p}\rangle$ is the inner product between eigenstates ${\boldsymbol{\Psi}_m}$ and ${\boldsymbol{\Psi}_p}$. The latter are orthogonal if $\langle{\boldsymbol{\Psi}_m}\textrm{|}{\boldsymbol{\Psi}_p}\rangle = 0$ and $\langle{\boldsymbol{\Psi}_m}\textrm{|}{\boldsymbol{\Psi}_m}\rangle$ represents the norm of ${\boldsymbol{\Psi}_m}$. The inner product constitutes a fundamental building block of standard (Hermitian) electrodynamics [1,2] since it allows the projection of any scattered wave on the basis of modes, $\langle{\boldsymbol{\Psi}_s}\textrm{|}{\boldsymbol{\Psi}_p}\rangle = {\alpha _p}\langle{\boldsymbol{\Psi}_p}\textrm{|}{\boldsymbol{\Psi}_p}\rangle$. In particular, the essential orthonormality of normal modes is rooted in all variational and perturbational methods. Physically, orthogonal and normalized field distributions enable a direct evaluation of the strength ${\alpha _m}$ with which a mode interacts with an incident wave.

Conservative systems are described by Hermitian (self-adjoint) operators, whose eigenstates form a complete orthonormal basis. They are studied within the mathematical framework of Hilbert spaces. However, non-Hermitian systems are described by non-self-adjoint operators whose eigenstates (the QNMs) do not form a Hilbert space. Therefore, the issues of completeness and orthonormality with a suitable inner product are more complex. Note that, once a suitable inner product has been defined, the quantity $\langle{\boldsymbol{\Psi}_m}\textrm{|}{\boldsymbol{\Psi}_m}\rangle$ is not positive definite as a norm should be (it is a complex number). Nevertheless, we still refer to this expression as the norm of ${\boldsymbol{\Psi}_m}$.

The question of decomposing the response of a non-Hermitian system as a superposition of QNMs has a venerable history, which probably started in quantum mechanics [6–11], see [24,25] for a recent review. Important developments have also been provided in the context of gravitational physics and general relativity [13,14,26]. Note however that in this case the scattering potential has an infinite spatial extent, which causes specific difficulties not encountered in other fields of non-Hermitian physics. On the other hand, electromagnetism is a more favorable environment for QNM theories since the scattering potential (e.g., a permittivity change $\mathrm{\Delta }\varepsilon $) has usually a compact support. The first studies in electromagnetism can be traced back to the seventies [26,27]. Let us emphasize that the vectorial nature of electromagnetic fields and the presence of a substrate or other complex environments raise specific issues that will be considered with care later on.

From the beginning, two fundamental questions have received considerable attention. Under which conditions are QNM expansions complete? How can we evaluate the amplitude of each mode in a superposition of QNMs considering that it is challenging to properly define an inner product and an orthonormalization for modes whose field diverges as $|\textbf{r} |\to \infty $? Two different frameworks have been developed for dealing with these fundamental issues.

Since the nice properties of Hilbert space are lost, the derivation of orthonormality is a long and winding road. Therefore, a first class of methods for deriving QNM expansions chose a framework radically different from the routine of Hermitian systems. Instead of relying on the definition of an inner product, this approach is based on the theory of functions of a complex variable, more particularly the residue theorem and one of its consequences, the Mittag-Leffler theorem [9,10,28–31]. Under certain conditions, a representative function of the system response (often its Green tensor) is a meromorphic function that can be decomposed as a sum of its poles, which are the QNMs complex frequencies. It is interesting to note that this approach does not explicitly rely on the definition of an inner product. Consistent with the terminology adopted in [3], we refer to this approach as the residue-decomposition approach.

On the contrary, the second approach for deriving QNM expansions closely follows the routine of Hermitian systems (the simple geometrical idea of projecting the system response), even if the vector space formed by the QNMs is not a Hilbert space. The regularization of the diverging fields is of particular importance in this approach for defining a correct inner product that leads to orthogonality and normalization [8,16,24,25,32–35]. The guiding principle is to “map” the exponentially divergent QNMs into bounded square-integrable modes, which can be considered as new vectors of a sort of generalized Hilbert space. The important (and physically appealing) consequence of this mapping is that regularized QNMs share many similarities with normal modes of conservative systems. Consistent with the terminology adopted in [3], we refer to this approach as the orthogonality-decomposition approach.

Before going any further, let us draw the reader’s attention to the existence of different modal theories that do not use QNMs. The main motivation for not using the natural resonances of the system is to retain a Hilbert space approach that avoids the issue of the divergence of the QNM field. The downside is that the ‘modes’ used in these alternative theories have less physical meaning. For instance, ‘permittivity modes’ can be used to calculate the scattered field; they are the eigenstates of the source-free wave equation with permittivity rather than frequency as the eigenvalue, see [36] and references therein. In another formalism, the system under consideration is not the scatterer alone but the closed ensemble formed by the input (source) space + the scatterer + the output (receiving) space. Then, the modes are pairs of functions that are calculated by the singular value decomposition of the scattering operator [37]. We will not further discuss these topics in this review.

1.4 Motivation

QNM completeness, orthogonality, and normalization have been recently addressed, see the Sections 4 and 5 of the review article [3] and the Sections 3 and 4 of the tutorial [21]. However, the topics covered in these articles are much broader than in this focused review, which is thus expected to deliver a more comprehensive message on these issues. More specifically, the main points that we develop are: (i) practical difficulties and limitations of the residue-decomposition approach (see Section 2.3), (ii) necessity to regularize the QNM field to allow the use of orthogonality-decomposition approaches (see Section 3), and (iii) comparison of different methods for normalizing QNMs (see Section 4). The latter is accompanied by new mathematical demonstrations and numerical evidences. The global vision we intend to provide is, to our knowledge, unique in the field of electromagnetic QNMs. It may didactically lead newcomers towards sound, concrete, and practical methods to calculate and use QNMs with their own approaches or with available software [38,39].

The puzzling issue of normalizing diverging fields that are not square integrable has lastingly hindered the development of QNM formalisms for analyzing nanoscale light confinement. The issue is theoretically solved nowadays, see for instance [3], at least in broad terms, owing to the intense efforts provided during the last decade. However, this has not been achieved in a progressive and linear way. As often in science, controversies have arisen along the way. For instance, the authors of [40] discussed different normalization methods and concluded that they provide the same result. This conclusion was vividly debated in a series of articles, comments, and replies [41–44]. When preparing this review, we realized that the literature, including earlier milestone publications, is quite confusing, sometimes misleading, and not free of mistakes. As a result, more recent works take some erroneous conclusions for granted, while correct results seem to be already forgotten. We try to clarify these points, while keeping the discussion accessible for a broad audience. We believe that it is instructive to look at how ideas have evolved to help newcomers who want to apply the QNM machinery and to promote further developments.

1.5 Structure

Besides the introduction, this review is divided up into five sections and is enriched by five appendices. Sections 2 and 3 address the construction of QNM expansions with the residue-decomposition approach and the orthogonality-decomposition approach. In general, the residue-decomposition approach works with the divergent QNMs of an open space (Section 2) while the orthogonality-decomposition approach requires regularized (i.e., non-divergent and square-integrable) QNMs (Section 3). Section 4 is devoted to the normalization issue. We first provide a historical perspective on the development of different normalization methods and then we compare their domain of validity. Section 5 summarizes and concludes the discussion.

Section 2 is devoted to the divergent QNMs of an open space. Section 2.1 first recalls the divergence issue that hinders the use of the classical inner product of conservative systems. Then, Section 2.2 presents the complexity of the spectrum (eigenvalues distribution in the complex frequency plane) of the open-space Maxwell operator. Section 2.3 describes the main lines of the construction of QNM expansions with the residue-decomposition approach. The conditions under which completeness is achieved are discussed. Finally, Section 2.4 further reviews the completeness issue.

In Section 3, we review the theoretical tools to regularize the divergent field of QNMs. The regularization enables us to define an inner product that leads to orthogonality and normalization. Importantly, we bridge the latest developments with earlier related methods explored in other fields than electromagnetism, especially in quantum mechanics. Section 3.1 first provides an overview from the non-Hermitian quantum mechanics literature, which has promoted similarity transformations to regularize QNMs. In electromagnetism, perfectly matched layers (PMLs) are extensively used since their birth in the 1990’s as a mathematical and numerical tool to implement outgoing wave boundary conditions in scattering problems at real frequencies. Using PMLs for the eigenproblem at complex frequencies is a quite natural road to build regularized QNMs. Section 3.2 presents the main ingredients of PML regularization in the continuous infinite space, while Section 3.3 addresses the key issue of PML regularization in a finite discretized space. Section 3.4 discusses important properties of regularized QNMs: norm, (bi)orthogonality, and completeness. The orthogonality-decomposition approach can be extended to systems involving coupling between Maxwell equations and other physical or chemical equations. Such problematics will probably occupy an increasingly important place in nanophotonics studies. The use of regularized QNMs for these multiphysics situations is addressed in Section 3.5.

Section 4 goes through the history of QNM normalization in optics, which can be separated in two main periods. In the 1990’s, P. T. Leung, K. Young, and coworkers at the Chinese University of Hong Kong played an important role in the development of QNM theory for optical cavities. They studied scalar fields in 1D cavities as well as vector fields in 3D systems with a spherical symmetry. Then, the 2010’s witnessed a second active period where different groups addressed the issue of QNM normalization, for arbitrary 3D geometries for which no analytical solution exists.

In view of this historical perspective, one can underline that a few general normalization methods have existed for about ten years to normalize the QNMs of any general structure, potentially made of dispersive, anisotropic, magnetic materials and surrounded by non-uniform media, as in the case with resonators on substrates for instance. The domain of validity of these different normalization methods is studied in detail. In addition, we correct inaccuracies or inconsistencies of the recent literature. In particular, we recall that the normalization method inherited from the works of the Hong Kong group is only valid for 1D systems and provides erroneous results for any other case, even 3D systems with a spherical symmetry [42,43].

Finally, let us briefly describe the content of the Appendices. Appendix A recalls the Lorentz reciprocity theorem, sometimes called divergence theorem, because it plays a central role in the review. We provide in Appendix B a proof of the invariance of the QNM norm calculated with PMLs. The review mostly considers the case of reciprocal materials. However, we show in Appendix C that the PML regularization allows the treatment of non-reciprocal materials as well. Appendix D addresses the numerical implementation of the normalization.

2. Divergent QNMs of the open space

In optics and more generally in electromagnetism, the QNMs of any open cavity are formally defined as the time-harmonic solutions to the source-free Maxwell equations with outgoing waves boundary conditions [3]

We use throughout the article the $\textrm{exp}({ - i\omega t} )$ convention and denote with a tilde QNM complex-valued frequencies, $\tilde{\omega } = \mathrm{\Omega } - i\mathrm{\Gamma }/2$, where $\mathrm{\Gamma }$ stands for the decay rate, i.e., the inverse of the mode lifetime $\tau = 1/\mathrm{\Gamma }$. The factor of $1/2$ accounts for the difference between amplitude and energy decays.

In Eq. (1), $\boldsymbol{\varepsilon }({\textbf{r},\omega } )$ and $\boldsymbol{\mathrm{\mu}}({\textbf{r},\omega } )$ are the permittivity and permeability tensors of the system under study. They are in general frequency-dependent and complex-valued tensors characterizing dispersive, absorptive, and anisotropic materials. We assume them to be equal to their transpose (reciprocal materials) for simplicity reasons. The specific case of non-reciprocal material is studied in Appendix C. We do not consider situations involving coupling between Maxwell equations and other physical or chemical equations, e.g., optomechanics or thermoplasmonics, but the following developments remain valid if the joint operator combining all differential equations under consideration admits an analytical continuation at complex-valued frequencies. Incidentally, note that when Drude or Lorentz permittivities are considered, multiphysics (electromagnetism and charge motion) is already enrolled in Eq. (1).

2.1 Divergence of the QNM field

The non-Hermitian character of QNMs has a thorny consequence: the exponential growth of the QNM field away from the resonator. It is easy to understand the origin of this divergence by considering a resonator embedded in a uniform medium with a refractive index n. On the one hand, the QNMs satisfy the outgoing-wave condition and, for large $r$’s, their fields take the asymptotic simple form of a progressive spherical wave, $\textbf{A}({\theta ,\varphi } ){r^{ - 1}}\exp [{ - i\tilde{\omega }({t - nr/c} )} ]$. On the other hand, the QNM fields have to exponentially decay in time. Then, $\textrm{Im}({\tilde{\omega }} )< 0$ and the fields exponentially diverge in space as $r \to \infty $ [45].

Despite what is occasionally read in the literature, the exponential growth does not raise any physical inconsistency. It follows from causality: QNMs are exponentially damped in time. Thus, they diverge as $t \to - \infty $, and since the QNM field radiated at $t ={-} \infty $ is perceived at $r \to \infty $ at a finite t, the QNM field diverges in space. The spatial divergence for $r \to \infty $ is not more unphysical than the temporal divergence for $t ={-} \infty $. In practice, the signal does not start at $t ={-} \infty $, the field is actually observed over finite spatial intervals, and the divergence is not met.

On another note, the divergence of the electromagnetic field represents a difficulty for the development of a non-Hermitian electrodynamic formalism. The concept of inner product constitutes a fundamental building block of standard (Hermitian) electrodynamics [1,2]. In particular, it leads to the essential orthonormality of normal modes that is rooted in all variational and perturbational methods based on the field expansion over a set of eigenstates. Indeed, the inner “wavefunctions” product of Hermitian electrodynamics $\smallint {\tilde{\textbf{E}}_m} \cdot \tilde{\textbf{E}}_n^\ast {d^3}r$ cannot be used in the non-Hermitian context. Furthermore, note that the unconjugated form $\smallint {\tilde{\textbf{E}}_m} \cdot {\tilde{\textbf{E}}_n}{d^3}r$ often used for the orthonormality of guided (bounded) modes in absorbing waveguides [17] cannot be used as such either, since the QNMs diverge as $r \to \infty $.

2.2 Spectrum of Maxwell operators

Figure 1(a) shows the spectrum (eigenvalues distribution in the complex frequency plane) of the open-space operator of Eq. (1). The spectrum comprises a discrete set of QNMs represented with red dots. They are often termed natural frequencies since they are frequencies for which the object can have a finite response even for an infinitesimal excitation by an incident wave. They are thus poles of the scattering operator. Owing to the causality principle, and ${\boldsymbol{\mathrm{\mu}}^\ast }({\textbf{r},\omega } )= \boldsymbol{\mathrm{\mu}}({\textbf{r}, - {\omega^\ast }} )$. Thus, if $({{{\tilde{\textbf{E}}}_m},{{\tilde{\textbf{H}}}_m},{{\tilde{\omega }}_m}} )$ is a source-free solution of Maxwell equations, $({\tilde{\textbf{E}}_m^\ast ,\tilde{\textbf{H}}_m^\ast , - \tilde{\omega }_m^\ast } )$ is another solution for the same permittivity and permeability distributions, see the Section 2.3 in [3]. This explains why every eigenfrequency $\tilde{\omega }$ in the lower right quadrant of the complex plane is paired with another twin frequency $- {\tilde{\omega }^\ast }$ in the lower-left quadrant, except those with a null real part.

 

Fig. 1. Complex frequency plane for a resonant open system. The problem is defined by the Maxwellian operator of Eq. (1) in an unbounded space with real spatial coordinates. The spectrum features accumulation points (black plusses), QNMs (red dots), branch cuts (blue solid curve). ${C_0}$ is a closed and finite counterclockwise oriented curve. It is chosen to encircle the point $z = \omega $ (often on the real axis) and a finite set of poles, ${\tilde{\omega }_m},m = 1, \ldots M$, but to avoid branch cuts.

Download Full Size | PDF

The spectrum may feature other singularities, such as accumulation points arising from either poles of the $\boldsymbol{\varepsilon }(\omega )$ or $\boldsymbol{\mathrm{\mu}}(\omega )$ (material resonance) [30,34,35,46–48], corner modes [39], or plasmonic resonances [34,35,46]. The latter appear in the case of metallic particles surrounded by a dielectric material and correspond to complex frequencies for which $\varepsilon (\omega )={-} {\varepsilon _b}$, with ${\varepsilon _b}$ the permittivity of the dielectric medium. Such accumulation points are not a mathematical curiosity. They are responsible for physical effects such as quenching, which is a fundamental and detrimental process for a molecule that emits light in the vicinity of a metallic material, see Fig. 3 in [34].

The spectrum may also feature branch cuts when the resonator is not surrounded by a uniform medium, a typical example being a resonator laid on a substrate. Different outgoing wave conditions should be satisfied in two different half spaces, the substrate and superstrate. The difficulty further increases for resonators surrounded by complex environments that potentially support guided waves. Typical examples are ring resonators, photonic-crystal cavities, or plasmonic nanoantennas coupled to a waveguide. The last case is analyzed in Figs. 3 and 4.

Finally, let us also mention spectral singularities, known as exceptional points, for which two or more QNMs coalesce for specific values of an opto-geometrical parameter. Such degeneracies have recently raised considerable attention in photonics, given that optical gain and loss can be used as nonconservative ingredients to create the right conditions for the apparition of exceptional points, which in turn result in new optical properties [49].

2.3 QNM expansion and completeness

To interpret experimental measurements, physicists like to consider how the resonance explicitly impacts the system, reconstructing the resonator response (at least in a compact subspace of ${\mathrm{\mathbb{R}}^3}$) with a QNM expansion of the scattered field

There is a large literature on QNM expansions, and whether the QNMs form a complete set or not (i.e., whether Eq. (2) converges or not) is considered as a crucial question that has commanded considerable attention in foundational electromagnetic studies [9,10,26,27]. It is interesting to summarize the main results, and for that purpose, we follow the residue-decomposition approach based on the pole expansion of meromorphic functions (Mittag-Leffler theorem) [50]. Equation (2) links a function of a real variable $\omega $ to its complex poles, and in complex analysis, this is achieved by evaluating line integrals of analytic functions over closed curves. Let us consider the integral

To analyze the question of completeness, one enlarges the contour ${C_0}$ so that all the QNM poles are encircled. The line integral of the first term on the right-hand side of Eq. (5) is now computed over a circle denoted ${C_\infty }$ with a center at $\omega = 0$ and an infinitely large radius [9]. If the integral over the circle ${C_\infty }$ is zero, i.e., if the function $\frac{{{\textbf{E}_S}({\textbf{r},z} )}}{{z - \omega }}$ vanishes sufficiently fast as $|z |\to \infty $ and if there are no branch cuts, we then end up with the desired expansion for which the response is set as a sum of independent pole contributions,

In general, the integral over ${C_\infty }$ is zero as long as the parameter $\textbf{r}$ belongs to a certain compact region and the QNM basis is then complete only for positions $\textbf{r}$ inside this region. The region for which the QNM expansion is complete has been extensively studied for simple problems for which the scattered field or the Green tensor is known with a closed form expression [9,10,11,16,19,26,27]. 1D Fabry-Perot resonators (e.g., potential wells in quantum mechanics) or 3D spherical Mie-resonators have been primarily analyzed. For 1D Fabry-Perot resonators, typically a 1D slab of permittivity $\varepsilon (x )= {\varepsilon _0}$ for $|x |< L/2$ in a uniform background, $\varepsilon (x )= {\varepsilon _1}$ for $|x |> L/2$, the QNM basis is complete inside the slab in the interval $\left] { - L/2;L/2} \right[$ [16]. For spherical resonators defined by $\varepsilon (r )= {\varepsilon _0}$ for $|r |< R$ and $\varepsilon (r )= {\varepsilon _1}$ for $|r |> R$, completeness is achieved inside the sphere of radius R [19]. These results are quite intuitive; at least, it is obvious to see that the field outside the cavity behaves as $\exp [{i\tilde{\omega }nr/c)} ]$, and thus exponentially diverges as $z \to \infty - i\infty $ on the lower half circle of ${C_\infty }$. From these initial studies, it is generally admitted nowadays that completeness is achieved inside the minimal sphere that surrounds the resonator, outside which the permittivity is uniform. However, a rigorous proof of the vanishing of the integral over ${C_\infty }$ inside the resonators with arbitrary geometries is still waiting to be done. More details can be found in the Section 4 and in the Appendix C in [21].

Disappointingly, the line integral over ${C_\infty }$ is not equal to zero in general and Eq. (6) is rarely met in practice. Indeed, most systems [51] involve a boundary between two media, defining a substrate and a superstrate, e.g., a metal particle deposited on a substrate or a photonic-crystal cavity in a semiconductor membrane. The outgoing-wave condition has to be satisfied into two different half spaces and two branch-cut singularities appear. The circle ${C_\infty }$ should then be distorted to avoid crossing the branches, as illustrated in Fig. 1, and the contour integral is not null. One may argue that the contribution of the branch cut may also be computed numerically by integration in the complex plane along the distorted contour ${C_\infty }$, but in practice this is not performed easily; it has only been achieved for very simple geometries where the branch-cut contribution can be analytically obtained or when all the poles can be computed [52].

We note that a precise knowledge of the spatial region for which the QNM basis is complete is not critical nowadays, since this knowledge is useful only for 1D or 3D resonators [53] in uniform backgrounds, which are more of academic than practical interest. We also stress that, in recent works for 1D Fabry-Perot resonators [54] and spheres [31], completeness is ensured inside and outside the resonator, see also [4] for gravitational waves. These results are obtained through a regularization based on a real coordinate transformation of both the spatial and temporal coordinates. The spatio-temporal regularization has not received a generalization for higher dimensions so far (note that, owing to the spherical-harmonic orthogonality, the Mie theory of spheres is similar to that of slabs illuminated with oblique incidences), but further studies would be highly appreciated.

2.4 Conclusion

The first line in Table 1 summarizes the discussion on the completeness for the open space. Completeness is guaranteed inside the resonator for geometries with uniform backgrounds only. The presence of a substrate results in branch cuts that are problematic. When the function is not known analytically, it is not easy to compute the branch-cut contribution along the deformed contours. Note that the contribution depends on the scattered field itself. As a result, even if the geometry is fixed, it changes as the driving source is modified. We are not aware of earlier works that accurately compute the branch-cut contribution, except for the specific case of 2D resonators in free space (without substrate) [52], for which the branch cut arises from the singularity of the Green tensor. This case is conveniently much simpler than the case of branch cuts arising from the presence of a substrate as will be analyzed in Fig. 3.

Table 1. Completeness of the QNMs in unregularized (open) and regularized spaces.

View Table | View all tables in this article

So, what to do in practice? The literature, as a whole, proposes to use a restricted finite subset of QNMs, and to complement this subset by another contribution that might be called a non-resonant contribution. The expansion may then be written as

Several approaches have been experimented to evaluate the non-resonant contribution. In the Riesz projection method [55,56], a finite contour, similar to ${C_0}$ in Fig. 1, is considered. The non-resonant contribution is accessed by computing $\textbf{E}({\textbf{r},z} )$ for complex frequencies z along this contour, which requires solving Maxwell equations rigorously typically one hundred times [55]. High accuracy is achieved in general; note however that the case with branch cuts has not been tested so far [57]. Alternatively, the non-resonant contribution may also be computed by directly solving Maxwell equations at a few real frequencies, and since it is expected to gently vary with the frequency, it can be easily interpolated with the help of the real-frequency data [58,59]. High accuracy has been demonstrated with this hybrid method for complex geometries, a dielectric grating [58] and a silicon nanodisk on a metal thin layer [59].

The previous approaches have an obvious practical interest, as they may enable fast, intuitive, and accurate methods for analyzing the responses of resonators. However, they lack the reassuring consistency provided by completeness. For instance, coupled-QNM theories for the coupling of a resonator with power-carrying channels [60] or with other resonators [61,62], take advantage of the convergence guaranteed by increasing the number of QNMs in the expansion. The same could be said also for the resonant-state expansion [28,29,47,63,64], a rigorous high-order perturbation method for simply computing the QNMs of any resonator by using the complete set of QNMs of a simpler geometry as a basis.

In the last few years, it has become clear that alternative formalisms relying on expansions that are complete inside and outside the resonators may be developed [33,34,65]. In this case, completeness is achieved by QNM regularization and by complementing an infinite set of QNMs by another infinite set of numerical modes (often called PML modes [66]). The expansion takes the following form

In Eq. (8), the numerical modes ${[{\tilde{\textbf{E}}_p^{num}(\textbf{r} ),\tilde{\textbf{H}}_p^{num}(\textbf{r} )} ]^T}$ are denoted just like the actual QNMs and their excitation coefficients also (we simply add a superscript ‘$num$’ for the sake of differentiation). The reason is that the numerical modes are computed exactly as the QNMs and that their excitation coefficients have a closed form expression, which is exactly the same as that of the QNMs.

Importantly, we note that, since the completeness of the expansion is guaranteed inside or outside the resonators, the formalism is particularly suited to various coupled mode theories, including coupling between scattering bodies and waveguide channels [60], hybridization of nanoresonators [61,62] and rigorous high-order perturbation theory for large shape deformation [67–69].

The present review is not concerned about formulas for the $\alpha $’s [3,65]. Let us mention that there is a single valid expression of the $\alpha $‘s for resonators made of non-dispersive materials for a given inner product [70,71], since the latter allows us to project ${[{{\textbf{E}_S},{\textbf{H}_S}} ]^T}$ onto the basis vectors ${[{{{\tilde{\textbf{E}}}_m},{{\tilde{\textbf{H}}}_m}} ]^T}$ [3,32]. The case of dispersive resonators will be briefly covered in Section 3.4.3.

The expansion in Eq. (8) fundamentally relies on the regularization of QNMs. We now examine its implication.

3. Regularized QNMs

Let us start by recalling that the concept of the inner product constitutes a fundamental building block of standard (Hermitian) electrodynamics or quantum mechanics in classical textbooks on conservative systems, see Chapter 10 in [50]. With inner products, one may conveniently normalize the vectors, construct an orthogonal basis and perform projections onto one vector of the basis.

Indeed, one anticipates that, if we may define a sort of inner product for non-Hermitian systems, one may then accommodate the tools of conventional Hermitian electrodynamics and considerably simplify the treatment of non-Hermitian system. This is exactly the approach that is highlighted in this section, which will end up with the definition of regularized QNMs and numerical modes, which all obey the same orthogonality product and together form a complete basis of a regularized Hilbert space of square-integrable functions.

The concept of QNM regularization has been initially developed in a series of theoretical studies in quantum mechanics, which will be reviewed in the Section 3.1. We will then focus on QNMs that are regularized with perfectly matched layers, first with an infinite thickness in Section 3.2, and then with a finite one in Section 3.3. Finite thicknesses correspond to the current practice in computational electrodynamics. This will give us the opportunity to well differentiate the spectra of continuous operators of infinite space with those of discrete operators of discretized and finite space. We will then conclude the Section by commenting on the completeness of the basis formed by regularized QNMs and numerical modes, and the different tests that are convincingly supporting the completeness.

3.1 Overview from the quantum mechanics literature

From the beginning of theoretical study on the spectral decomposition of scattered waves into a superposition of QNMs, a fundamental question has drawn considerable attention: How may we properly “map” the divergent QNMs of the open space onto a Hilbert space and force the QNMs to be square-integrable. The general recipe, widely developed in quantum mechanics and gravitational physics, consists in carrying out a suitable mathematical mapping that converts the original exponentially divergent QNMs into bounded square-integrable QNMs, which can be considered as new vectors of a new Hilbert space that shares many similarities with the normal modes of conservative systems.

In this overview section, we adopt the presentation of Chapter 5 of the monograph in [25], which emphasizes mappings with similarity transformations. For the sake of easy reading by an interested reader, we symbolically recast the continuous Maxwellian operator of Eq. (1) into the classical Hamiltonian format $\textbf{H}\tilde{\boldsymbol{\Psi}} = \tilde{E}\tilde{\boldsymbol{\Psi}}$ of quantum mechanics, with $\tilde{E} \equiv \hbar \tilde{\omega }$ the complex energy and $\tilde{\boldsymbol{\Psi}} \equiv [{\tilde{\textbf{E}};\tilde{\textbf{H}}} ]$ the wavefunction. The similarity transformation

The key point is how to find such a linear mapping. The literature on non-Hermitian quantum mechanics contains various approaches in which such a transformation is introduced, starting with the Gaussian regularization introduced by Zel’dovich in the 1960s [8,72] and recently revisited in [73] for Maxwellian operators. Most of the approaches are in fact rather mutually equivalent. They are often applicable for simple analytical potentials (scalar waves in 1D) and are inconveniently transferable to contemporary electromagnetic resonators encountered in nanophotonics.

The motivation for applying the similarity transformation is obvious. Under this transformation, the resonance wavefunctions initially exponentially growing up in space become ‘regular’ functions that are square-integrable. One immediately anticipates that normalization and orthogonality may then be defined in the new Hilbert space of regularized functions.

In electromagnetism, the preferred path to build regularized QNMs is that provided by perfectly matched layers (PMLs). Owing to the importance of scaled spaces for deriving the orthonormality relation, the PML machinery will be carefully documented in the next subsection. We simply note, at this stage, that a PML does not provide a one-to-one mapping between the divergent and regularized QNMs, as Eq. (9) simplistically suggests.

As an aside, let us note that another regularization approach has been recently introduced in the literature on QNMs [74,75] and further utilized in a series of applications, including the important case of QNM quantization [76,77]. The regularization relies on a dual scheme, with the embodiment of a QNM field evaluated at the real-valued resonance frequency of the QNM outside the resonator while keeping the actual QNM field evaluated at the complex-valued frequency inside the resonator. The divergence is therefore removed by construction. The heuristic is simple and removes all divergence issues with a wave of a magic wand. For dispersive media, the regularized QNM no longer satisfies the interface conditions for tangent electromagnetic fields at the resonator boundary. More fundamentally, the heuristic does not preserve the machinery of analytic continuation, which is central in QNM theory, and fundamental notions such as completeness and orthogonality are not preserved. In our opinion, it is a quick fix that requires deeper considerations.

3.2 PML-regularization in infinite space

Usually, PMLs are used for damping oscillating waves corresponding to scattered fields at real frequencies. Here we review how they can be used for the regularization and normalization of QNMs. A convenient and efficient method to divert the divergent QNM into the physical domain of square-integrable wavefunctions is to rotate the coordinate along which the divergence occurs into the complex plane by an angle $\theta $.

To guide the analysis, Fig. 2(a) shows a 2D resonator on a complicated substrate. On purpose, the “substrate” is composed of two semi-infinite uniform half-planes with different permittivities, see the pink and yellow domains. Importantly, we note that for large enough values of x, say, for $|x |$ larger than some positive coordinate ${x_0}$, the permittivity is invariant in x; it depends only on y. Conversely, the same remark can be made for the y coordinate. The key fact about these geometries that already encompass a great variety of practical cases is that the solution $\tilde{\boldsymbol{\Psi}} = [{\tilde{\textbf{E}};\; \tilde{\textbf{H}}} ]$ of Maxwell equations is an analytic function of the spatial coordinate $\textbf{r}$ (see for instance [78]), which can be freely analytically continued and rotated at complex values of $\textbf{r}$ with PMLs [79,80].

 

Fig. 2. Complex mappings with PMLs. (a) Illustration of a PML implementation in 2D for a resonator with a complicated “substrate” composed of two semi-infinite uniform orthogonal half planes with different permittivities. Importantly, we note that the PMLs are implemented by complex mappings that can all be factorized as the product of a function of x or y only. (b) Example of an appropriate trajectory in the complex coordinate plane that is classically used in the PML machinery. The key point is that the trajectory encompasses a linearly growing imaginary part for $\tilde{X} \to \infty + i\infty $. Note that no mapping is applied for $|x |< {x_0}$, implying that the regularized QNMs defined by Eq. (11) and the initial diverging ones defined by Eq. (1) coincide in this bounded region of space. (b) is adapted from [16].

Download Full Size | PDF

Figure 2(a) illustrates a possible PML implementation composed of three different regions. In the center part, no coordinate transform is applied; the physical space is unmapped. In the regions along x or y axis, only a single direction is mapped; finally, at the corners, both the x and y coordinates are mapped, independently. This complex rotation that keeps the coordinates unscaled in the resonator region where the permittivity is complicated is referred to as the exterior scaling transformation in quantum mechanics studies; it has been introduced in the late 70s, see the Section 5.3 in [25].

The following analysis can be conducted for the 2D geometry of Fig. 2(a) by detailing all the PML regions [78], but we will adopt a simplified analysis, by just highlighting the main idea. First, we note that one moves from any rectangular region to a neighbor one in Fig. 2(a), by applying or removing a single coordinate transform in a single direction, either x or y. Therefore, for the following analysis, we consider that the mapping is performed along a single coordinate, denoted x, hereafter. Furthermore, we consider a mapping for $x \to + \infty $ only, the other Cartesian coordinates and $x < 0$ being treated similarly. By hypothesis [see Fig. 2(a)], the permittivity and permeability are both independent of x for $x > {x_0}$. Thus, the QNM in infinite space may be decomposed as a superposition of waves of the form ${\boldsymbol{\Psi}_k}({y,\; z} )\textrm{exp}({i{{\tilde{k}}_x}x - i\tilde{\omega }t} )$ (modal expansion with radiation and guided waves) [78] or simply planes waves if the half space $x > {x_0}$ is uniform for any $y,\; z$ (Helmholtz–Weyl decomposition).

Let us first consider the field $\boldsymbol{\Psi}({x,y,z} )$ scattered by the resonator under illumination by a driving field at real frequency $\omega $. Since $\boldsymbol{\Psi}({x,y,z} )$ satisfies the outgoing wave condition, $\boldsymbol{\Psi}({x,y,z} )$ gently decays as $1/r$ and oscillates at large distances r from the resonator. Thus, $\boldsymbol{\Psi}$ admits an analytical continuation $\boldsymbol{\Psi}({\tilde{X},y,z} )$ in the complex domain defined by $\textrm{Re}({\tilde{X}} )> {x_0}$ (${x_0}$ being the inner boundary of the scaled exterior) and $\textrm{Im}({\tilde{X}} )> 0$, with $\boldsymbol{\Psi} \to 0$ as $\tilde{X} \to \infty + i\infty $. The Maxwell equations, initially restricted to the real axis $x = \textrm{Re}({\tilde{X}} )$, can be solved in this infinite complex space along an appropriate trajectory with, for instance, a linearly growing imaginary part, see Fig. 2. This is the essence of the PML concept [79,80]. We adopt the presentation and notation of [78] hereafter.

The key point is that, owing to the choice $\textrm{Im}({\tilde{X}} )> 0$, only the outgoing waves admits an analytical continuation, since all propagative or evanescent waves of the form ${\boldsymbol{\Psi}_k}({y,\; z} )\textrm{exp}({ - ikx - i\omega t} )$ exponentially diverge: $\textrm{exp}({ - ik\tilde{X} - i\omega t} )\to \infty $ for $\tilde{X} \to \infty + i\infty $. Thus, by construction, PMLs automatically implement the outgoing wave condition. Additionally, since the analytical continuation guarantees that no reflection occurs at $\textrm{Re}({\tilde{X}} )= {x_0}$, PMLs (with an infinite thickness as we investigate in this Section; finite-thickness effects are considered in the next Section) perfectly satisfy the outgoing wave conditions by damping the field without reflection, at least for many geometries of wide interest [81]. They are often called reflection-less absorbing materials and this is why they have revolutionized boundary conditions in electromagnetism and other wave equations.

The integration along a trajectory with a complex coordinate is inconvenient. Thus, we conveniently change variables and map the complex trajectory into a new one depending only on a real variable, $x^{\prime} = \textrm{Re}({\tilde{X}} )$ [82]. The mapping requires to know the relationship between $\textrm{Im}({\tilde{X}} )$ and $\textrm{Re}({\tilde{X}} )$, i.e., the function $\tilde{X} = F({x^{\prime}} )$. Owing to the equivalence between material parameters and spatial coordinate transforms [83,84], the mapped Maxwell equations can be either written with transformed operators by replacing $\frac{{\partial ({\cdot} )}}{{\partial \tilde{X}}}$ by $\frac{{\partial x}}{{\partial \tilde{X}}}\frac{{\partial ({\cdot} )}}{{\partial x}}$ in Eq. (1) or by using transformed permittivity and permeability tensors, or by a combination of both approaches [78].

Figure 2(b) shows the simplest possible mapping corresponding to $\tilde{X} = x^{\prime}$ for $x^{\prime} < {x_0}$ and $\tilde{X} = x^{\prime} + i\tan (\theta )({x^{\prime} - {x_0}} )$ for $x^{\prime} > {x_0}$ ($x^{\prime}$ and x coincides in this simple case). This transform defines a rotational angle $\theta $ of the coordinates in the complex plane, $\tan (\theta )= \textrm{Im}({\tilde{X}} )/({x^{\prime} - {x_0}} )$. Other convenient mappings involve complex coordinate stretching, equivalent to an additional change of the length units. Dispersive coordinate transforms [80] with a frequency-dependent stretching factor are discussed later on.

Let us come back to the QNM solution of Eq. (1) with a complex frequency $\tilde{\omega }$. Just like for real frequencies, the PML machinery can be applied. Using the transformed permittivity and permeability tensors, we may symbolically rewrite the regularized spectral problem of Eq. (1)

There is however a minor, albeit essential, difference with the usual real-frequency case: QNM fields diverge as $\textrm{exp}({i\tilde{k}x - i\tilde{\omega }t} )\propto \textrm{exp}({ - \textrm{Im}({\tilde{k}} )x} )$ at large distances in real spaces. Thus, only the QNMs with an exponentially damped field along the trajectory in the complex plane are transformed into square integrable wavefunctions. This happens only if $\textrm{exp}({i\tilde{k}\tilde{X}} )\to 0$ when $\tilde{X} \to \infty + i\infty $. For a vacuum, $\tilde{k} = \tilde{\omega }/c$, the necessary condition is met only for QNMs with a quality factor $Q ={-} \frac{1}{2}\textrm{Re}({\tilde{\omega }} )/\textrm{Im}({\tilde{\omega }} )$ that verifies

Therefore, in the mapped space of regularized QNMs, only a subset of the initial complex QNMs emerges [33]. Figure 3(a), which is directly inspired from Fig. 2 in [33], illustrates our purpose. Only the QNMs with large quality factors emerge in the new complex frequency plane of the mapped space. They have complex frequencies that are identical to the initial ones obtained with an unscaled calculation, see Fig. 1. In the next Section 3.3, we will clarify what is happening with the veiled QNMs in the hatched half plane.

 

Fig. 3. Spectrum of regularized spaces. (a) Complex frequency plane after regularization by a non-dispersive PML with an infinite thickness. Owing to the rotation in the complex plane, the continuous spectrum composed of propagative or evanescent radiation modes (the real axis in the initial unmapped space) is rotated by an angle $\theta $. Its rotated version is shown with the dashed green line. The latter divides the plane in two half planes. Above the line, the initial QNMs are unveiled by the regularization and the spectrum is unchanged. Below the line in the shaded half plane, the initial QNMs are veiled and the initial spectrum is largely impacted. Inspired by Fig. 2 in [33]. (b) Complex-frequency plane of a discretized operator for a quite complicated geometry, a silver-Drude nanorod (plasma frequency ${\lambda _p} \equiv 2\pi c/{\omega _p} = $138 nm, 60-nm diameter, 90-nm height) lying on a 138-nm-thick silicon slab in air. Several branch cuts are discretized; two of them are marked with black arrows. Modes corresponding to QNMs are shown with a circle. They are independent of the PML. Numerical modes are shown with square marks. The red color indicates the excitation rate of the modes (here upon illumination by a plane wave), dark and light red corresponding to high and weak excitations, respectively. The spectrum was computed with the QNMEig solver using a non-dispersive PML. After Fig. SI.7(a) in [34].

Download Full Size | PDF

Importantly, note that the regularized QNM fields $\tilde{\boldsymbol{\Psi}} = [{\tilde{\textbf{E}};\; \tilde{\textbf{H}}} ]$ are different from the original diverging fields that are solutions of Eq. (1). They only coincide in the unscaled interior segment $|x |< {x_0}$. Another important remark concerns the fact that the QNM field in the real (physical) space $|x |> {x_0}$ cannot be recovered from the QNM field calculated in the PML for $|{x^{\prime}} |> {x_0}$. That would be possible only if a one-to-one mapping would exist between x and $x^{\prime}$. And this does not happen since a PML is a branch cut in the complex plane, which helps satisfying the outgoing wave condition. De facto, the one-to-one mapping only exists between the complex trajectory $\tilde{X} = F({x^{\prime}} )$ and $x^{\prime}$. It is thus important to realize that the PML regularization does not implement a similarity transform between the initial divergent QNMs and the regularized ones (the operator $\textbf{U}$ of Eq. (9) is meaningless in the present context). In our opinion, the PML machinery that generally applies to a wide range of geometries and not just to analytical functions is a much more abstract concept than the simplified similarity transform approach introduced in Section 3.1.

We know from Section 2 that the entire set of QNMs provides completeness only on rare occasions. Now we see that only an infinite subset of QNMs is unveiled by the regularization. Thus, the question of the completeness of regularized QNMs arises. That is where the numerical modes come in, as we discuss in Section 3.4.3. Let us first illustrate what are these numerical modes.

3.3 Practical PML-regularization with finite and discretized spaces

Except for some analytical cases, such as 1D Fabry-Perot resonators or spherical Mie potentials, the complex mapping with an infinite spatial extend is truly inconvenient. In practice, QNMs are computed with numerical methods [85]. The numerical space is always finite and owing to the ‘discretization’, the operators become matrices with a finite size. Numerical PMLs with finite thicknesses are similar to the PMLs introduced in the previous section. They additionally feature perfectly conducting magnetic or electric walls on the outer boundaries of the PML space (Neumann or Dirichlet boundary condition for a finite value $x^{\prime} > {x_0}$), so that the mapped space becomes bounded. To account for this change, we symbolically rewrite the regularized spectral problem of Eq. (11) as

Indeed, since the regularized QNMs are exponentially decaying in the PML region, most of them are weakly impacted by truncating the mapped space by hard wall: only their tiny exponential tails “see” the wall and reflect off it, before being attenuated again on the way back towards the region of interest $x^{\prime} < {x_0}$. However, some QNMs have small quality factors and are diverging extremely fast. Their regularized fields are not damped enough at the outer PML boundary and are quite impacted by the hard-wall boundaries. Thus, the spectrum of the new bounded operator with hard walls in Eq. (13) may be quite different from the one with infinite PMLs in Eq. (11). However, importantly, the actual QNMs that are relevant to understanding the physics of the resonators in some frequency range of interest are hopefully preserved by the regularization and the space truncation.

Figure 3 illustrates the significant impact on the spectrum of the passage from a continuous operator in an unbounded space to a discretized one in a bounded space. The spectrum corresponds to a silver nanorod on a thin silicon layer. It was computed with the solver QNMEig operating with the COMSOL environment with non-dispersive PMLs, see details in [34].

The spectrum encompasses QNMs shown with blue circles and numerical modes shown with squares. For this complex geometry, the numerical modes and the QNMs sit close to one another in the complex-frequency plane. To discriminate them, the authors in [34] have exploited the fact that the QNMs are insensitive to PML-parameter variations (as long as they are revealed by the complex coordinate transform), whereas numerical modes are sensitive [86]. They have thus performed two computations for two different PML parameters and have identified the eigenvalues that are independent of the PML. Further details are provided in the Supplemental Section 3.4 in [34].

The second salient feature is that the spectrum includes several nearly vertical aligned rows of numerical modes that are regularly spaced along the real axis. The PML implements the outgoing wave conditions for the guided modes of the slab, and just like for free-space plane waves in uniform media, differentiating ingoing and outgoing guided modes requires a cut in the complex plane. This cut is automatically implemented by the PML (at least for some specific geometries where the phase velocity is opposite to the group velocity [81]) and the aligned rows of numerical modes are thus interpreted as discretized branch cuts operated by the different waveguide modes [87]. Note that similar (discretized) branch patterns are also observed for gratings due to the passing-off of diffraction orders at cutoff frequencies determined by the grating equation, see Fig. 11 in [33] and Fig. 1 in [88].

Some numerical modes, especially those with a negative frequency ($\textrm{Re}({\tilde{\omega }} )< 0)$, have positive imaginary parts ($\textrm{Im}({\tilde{\omega }} )> 0)$, implying that they exponentially diverge as $t \to \infty $. The existence of these unphysical modes (we are not considering amplifying media) is due to the fact that the non-dispersive PML used for the computation in Fig. 3 does not satisfy the causality property ${\boldsymbol{\varepsilon }^\ast }({\textbf{r},\omega } )= \boldsymbol{\varepsilon }({\textbf{r}, - {\omega^\ast }} )$ and ${\boldsymbol{\mathrm{\mu}}^\ast }({\textbf{r},\omega } )= \boldsymbol{\mathrm{\mu}}({\textbf{r}, - {\omega^\ast }} )$. Non-dispersive PMLs do not prevent the obtention of a high accuracy when the field scattered by a driving pulse is reconstructed in the temporal domain. Indeed as $t \to \infty $, the reconstruction becomes false, but energy has leaked away from the resonator well before divergence effects become observable, see the discussion related to Fig. 2(b) in [34].

To obtain poles only with negative imaginary parts, one may use dispersive PMLs that implement a complex coordinate stretching with a stretching coefficient that is inversely proportional to the frequency [79,80]. All frequencies are treated on equal footing, the permittivity and permeability satisfy ${\boldsymbol{\varepsilon }^\ast }({\textbf{r},\omega } )= \boldsymbol{\varepsilon }({\textbf{r}, - {\omega^\ast }} )$ and ${\boldsymbol{\mathrm{\mu}}^\ast }({\textbf{r},\omega } )= \boldsymbol{\mathrm{\mu}}({\textbf{r}, - {\omega^\ast }} )$, and poles are obtained pairwise with frequencies $\tilde{\omega }$ and $- {\tilde{\omega }^\ast }$. Dispersive PMLs translate (instead of rotating) the continuum spectrum, which is useful for computing low-frequency modes (such as static modes). They are thus expected to offer important benefits compared to non-dispersive PMLs. However, according to the experience of some of the authors, the quality and convergence of the reconstruction are not improved. Some special modes, not observed with non-dispersive PMLs, appear and densely fill some frequency regions of the complex plane. In addition, the numerical implementation for computing the QNMs with dispersive PMLs is more complex; thus, the expected advantages seem not to be met in practice. We let the readers try both PMLs and form their own opinion.

In Fig. 3, all the modes are colored and sorted by increasing order of their excitation coefficient ${\alpha _m}$ in Eq. (8). Dark red coloring corresponds to a strong impact. QNMs are mostly impacting, but we also note some numerical modes arising from the branch cuts which significantly contribute to the reconstruction. This shows that the branch cuts cannot be neglected in general, and that leakage into the membrane mode is a dominant channel in this example.

Further note that some numerical modes shown with a square, like those labelled by the letter $A^{\prime}$, are nearly independent of the PML but are not considered as QNMs since the imaginary part of their frequency is slightly positive owing to the use of non-dispersive PMLs. These are numerical versions of the $- {\omega ^\ast }$ modes associated to the QNMs labelled by the letter A.

3.4 Fundamental properties of regularized QNMs

3.4.1 Norm of regularized QNMs

The norm of regularized QNMs relies on a volume integral of the form $\smallint \left[ {\tilde{\textbf{E}}\cdot \frac{{\partial \omega \boldsymbol{\varepsilon }}}{{\partial \omega }}\tilde{\textbf{E}} - \tilde{\textbf{H}}\cdot \frac{{\partial \omega \boldsymbol{\mathrm{\mu}}}}{{\partial \omega }}\tilde{\textbf{H}}} \right]{d^3}\textbf{r}$, see Section 4 and [32]. The integral is undefined for the unscaled real space owing to the QNM divergence, but it is defined for the regularized QNMs that are square integrable since they are exponentially damped in the mapped space for large $\textbf{r}$ (we should say ${\textbf{r}^B}$ but remember that we dropped the superscript). The question arises as to whether the integral over the infinite regularized space depends on the complex mapping (or PML) used?

The key point is that the integral, and therefore the norm, is independent of the trajectory used, as long as the inequality in Eq. (12) is satisfied. Since the integral over the complex trajectory is equal to the integral over the mapped space $x^{\prime}$ (the demonstration that relies on the equivalence between material parameters and spatial coordinate transforms [83] is not repeated here but can be found in the Supplementary Section 2 in [32]), the PML norm computed over the entire space is independent of the PML considered. Appendix 2 provides another demonstration of this key property. The demonstration is very simple and can be followed by all the readership.

The norm is thus a unique quantity, which is intrinsic to the QNM and leads to an effective and quite general method to compute normalized QNMs in practice [32,34,85], see Section 4 for a comparison of the PML-norm with other rigorous norms.

3.4.2 Biorthogonality

The concept of orthogonality is essential, and its demonstration is quite trivial for regularized QNMs. We will consider the QNMs of nanoresonators built with non-dispersive materials for the sake of simplicity. For dispersive materials, the permittivity depends on the frequency and Eq. (11) corresponds to a nonlinear eigenvalue problem. To derive the orthogonality relation, one needs to linearize the operator [3]. This can be performed with auxiliary fields [34,89]. In general, the auxiliary fields are physical quantities that represent the material polarizations $\textbf{P}$ and currents $\textbf{J}$ and the biorthogonality relation holds between QNMs enlarged with the auxiliary fields, $\tilde{\boldsymbol{\Psi}} = [{\tilde{\textbf{E}};\; \tilde{\textbf{H}};\; \tilde{\textbf{P}};\; \tilde{\textbf{J}}} ]$, which combine geometrical and material resonances [90].

The following derivation holds for non-dispersive materials [32]. It is repeated here to show that regularized spaces and their regularized eigenvectors can be treated in much the same way as the normal modes of Hermitian system. In order to realize this, please compare the following derivation with that of Section entitled “General Properties of the Harmonic Modes” in [2].

Let us consider two regularized QNMs ${\tilde{\boldsymbol{\Psi}}_n} = [{{{\tilde{\textbf{E}}}_n};\; {{\tilde{\textbf{H}}}_n}} ]$ and ${\tilde{\boldsymbol{\Psi}}_m} = [{{{\tilde{\textbf{E}}}_m};\; {{\tilde{\textbf{H}}}_m}} ]$, therein satisfying the same source-free Maxwellian operator, Eq. (11). We simply apply the Lorentz reciprocity theorem, taking the QNM ${\tilde{\boldsymbol{\Psi}}_n}$ with a complex-valued frequency ${\tilde{\omega }_n}$ as solution 1 and the QNM ${\tilde{\boldsymbol{\Psi}}_m}$ with a complex-valued frequency ${\tilde{\omega }_m}$ as solution 2. Since there is no source, we obtain

If two QNMs have equal frequencies ${\tilde{\omega }_n} = {\tilde{\omega }_m}$, then they are degenerate and are not necessarily orthogonal. For two modes to be degenerate requires that two different field patterns happen to have precisely the same frequency. Usually, there is a symmetry that is responsible for the astonishing coincidence. For example, if the permittivity and permeability distributions are invariant under a 90° rotation, QNMs that differ only by a 90° rotation are expected to have the same frequency. Such QNMs are degenerate and are not necessarily orthogonal. However, since the Maxwellian operator is linear, any linear combination of these degenerate QNMs is itself a mode with that same frequency. We can always choose to work with linear combinations that are orthogonal [65].

3.4.3 Completeness of expansions including numerical modes

So far, we have examined the completeness of expansions composed of ‘true’ QNMs and argued that completeness is achieved for very limited geometries only. In this Section, we examine the critical issue of the completeness of expansions involving regularized QNMs and numerical modes in ‘PMLized’ spaces.

The theoretical case of an unbounded space with infinitely thick PML has been scarcely considered in the literature to our knowledge. Few studies have been carried out to understand what happens to the veiled QNMs below the rotated continuous spectrum, whether they are replaced by new poles, and what is the nature of these poles that are not ‘true’ QNMs. It would be interesting to determine the entire spectrum of mapped operators in infinite spaces and possibly their branch cuts [33,92].

In fact, most recent studies have focused on practical cases, i.e., on mapped spaces with finite-thickness PMLs bounded by perfectly conducting electric or magnetic walls and discretized operators (i.e., non-Hermitian matrices). The theoretical background, essentially linear algebra, is much simpler than for continuous operators defined on open spaces. In contrast to a Hermitian matrix, a non-Hermitian matrix does not have an orthogonal set of eigenvectors; in other words, a non-Hermitian matrix A can in general not be transformed by an orthogonal matrix Q to diagonal form $D = {Q^\mathrm{\ast }}AQ$, where the star is the transpose conjugate. However, most non-Hermitian matrices can be transformed by a nonorthogonal matrix X to diagonal form $D = {X^{ - 1}}AX$ through their right and left eigenvectors [93]. There exist rare exceptions, for which X converges to a singular operator. This case corresponds to defective matrices is often referred to as non-Hermitian degeneracy or self-orthogonality, see for instance the Chapter 9 in [25] and [94].

Except for these exceptional cases for which two or more complex eigenvalues cross and the associated eigenvectors coalesce, the literature generally assumes that the whole set of QNMs and numerical modes (also termed as PML modes or Bérenger modes) form a complete set. This has been demonstrated for problems for which the Green’s tensor is analytically known [92]. However, although no rigorous mathematical proof formally exists, the conjecture is supported by many numerical examples.

Though they play a key role in the analysis of non-Hermitian matrices [50], left eigenvectors (or left QNMs) have not been considered so far in this review. For instance, only right QNMs appear in the biorthogonality product of Eq. (15). The reason is that we have considered only reciprocal materials ($\boldsymbol{\varepsilon } = {\boldsymbol{\varepsilon }^\textrm{T}}$ and $\boldsymbol{\mathrm{\mu}} = {\boldsymbol{\mathrm{\mu}}^\textrm{T}}$ where the superscript ‘$\textrm{T}$’ stands for transpose) so far. In this case of broad practical interest, the right and left QNMs (or eigenvectors) are equal and it is therefore not appropriate to distinguish them. For further reading, please refer to [60,95] and Appendix C, where the normalization of QNMs of resonators with non-reciprocal materials is considered.

Convincingly, accurate reconstructions of the scattered field with QNM expansions augmented by numerical modes were first provided for a 2D dielectric triangular rod, for which convergence of the absorption cross-section spectrum was reported by increasing the total number of QNMs and numerical modes retained in the expansion, see Fig. 7 in [33].

For resonators made of dispersive materials, the eigenvalue problem defined in Eq. (13) is no longer linear, since the operator itself depends on the frequency via $\boldsymbol{\varepsilon }({\tilde{\omega }} )$ or $\boldsymbol{\mathrm{\mu}}({\tilde{\omega }} )$. A possible approach to recover linearity is to incorporate auxiliary-fields, for instance, the material polarizations $\textbf{P}$ and currents $\textbf{J}$. It is then possible to recover a biorthogonal relation (see Eq. (4) in [34]) for both QNMs and numerical modes. The biorthogonality can then be used for projection on the basis elements and an expression for the excitation coefficients, the $\alpha $’s, can be further derived. Note that, like for any Hermitian matrix, once the source term is defined, there is a unique closed-form expression for the $\alpha $’s per biorthogonal relation [34]. However, as shown in [65], new formulas can be found by choosing different auxiliary fields for linearization or different methods for splitting the source term [96].

Figure 4 shows the convergence of the extinction cross-section spectra of the silver nanorod for a broad range of frequency of six octaves. Despite the series of branch cuts and accumulation points (at the poles of the Drude permittivity and at the surface plasmon frequency), convergence is evidenced as the total number M of numerical modes and QNMs retained in the expansion of Eq. (8) is increased. Reference [34] contains other numerical evidence of the convergence for complicated geometries, such as the accurate reconstruction of light scattered by bowtie antennas in the temporal domain with short pulses or quenched emission at frequency close to the surface plasmon frequency and its accumulation point.

 

Fig. 4. Numerical evidence of the completeness of expansions with QNMs and PML modes. (a) Scattering cross-section spectrum of the silver nanobullet of Fig. 3(b). The solid-blue and dashed-black curves are computed by retaining 20 QNMs only and 20 QNMs plus 200 numerical modes, respectively, and the shadowed pink curves are the data obtained with the frequency-domain solver of COMSOL Multiphysics. After [34]. (b) Convergence rate as the total number M of QNMs and numerical modes retained in the expansion in Eq. (8) is increased. The vertical axis represents the spectrally-averaged ($0.1{\omega _p} < \omega < 0.6{\omega _p}$) error between the predictions of the extinction cross section obtained with the truncated modal expansion of Eq. (8) and fully-vectorial data computed with COMSOL Multiphysics. The results are obtained with discretized spaces regularized with finite PML thicknesses and the modes are sorted by increasing order of their impact on the reconstruction. Adapted from Fig. SI.7 in [34]. Other relevant examples illustrating completeness can be found in [33,65,88,95].

Download Full Size | PDF

Other convincing reconstructions have also been obtained for other complicated geometries (e.g., metallic gratings that also have accumulations points and multiple branch cuts at grating anomalies) [88]. In [95], the convergence of the expansion is tested for a 2D resonator with a prismatic shape and a Lorentz permittivity for an even broader spectral interval, though a small inaccuracy is reported close to the pole of the permittivity with $M = 200$ modes. In [65], the authors have considered simpler geometries, namely 2D cylinders and spheres in free space. Using dispersive PMLs, they have numerically verified that the expansion with QNMs and numerical modes is complete, inside or outside the scatterer, i.e., everywhere in the physical domain that corresponds to the unmapped subspace [97]. The authors have also checked that an expansion based on QNMs only is incomplete even inside the resonator, owing to branch-cut of the 2D Green function. See Fig. 2 in [65] for more details.

Table 1 summarizes the discussion.

3.5 Multiphysics

The regularization approach with mapped spaces can be applied to problems combining electromagnetism and other areas of physics. First, it is worth noting that, for dispersive materials, multiphysics is already enrolled in the regularization approach with auxiliary fields. This is because the dispersions of $\boldsymbol{\varepsilon }$ and $\boldsymbol{\mathrm{\mu}}$ are determined by the microscopic interaction between electromagnetic fields and other elementary excitations, such as phonons, polaritons, magnons, plasmons [100].

For dispersive materials, $\boldsymbol{\varepsilon }$ and $\boldsymbol{\mathrm{\mu}}$ depend on the frequency and Eq. (1) corresponds to a non-linear eigenvalue problem. One then encounters a difficulty in deriving an orthogonality relation, which prevents us from obtaining a closed-form expression for the modal excitation coefficient of the scattered fields (${\textbf{E}_s}$, ${\textbf{H}_s}$), as the excitation coefficient of one QNM depends on the other QNMs [32]. This issue can be elegantly solved by introducing auxiliary fields that allow the reformulation of the non-linear eigenvalue into a linear one [34,89,101]. For example, considering a medium with a dispersive permittivity described by a single-pole Lorentz model, $\varepsilon (\omega )= {\varepsilon _\infty } - {\varepsilon _\infty }\omega _p^2{({{\omega^2} - \omega_0^2 + i\omega \gamma } )^{ - 1}}$, two auxiliary fields, the polarization $\textbf{P} ={-} {\varepsilon _\infty }\omega _p^2{({{\omega^2} - \omega_0^2 + i\omega \gamma } )^{ - 1}}\textbf{E}$ and the current density $\textbf{J} ={-} i\omega \textbf{P}$, can be introduced to reformulate Eq. (1) into a higher dimensional linear eigenvalue problem. An orthogonality relation, $\int\!\!\!\int\!\!\!\int_{\mathrm{\Omega } \cup {\mathrm{\Omega }_{PML}}} {{\varepsilon _\infty }{{\tilde{\textbf{E}}}_n} \cdot {{\tilde{\textbf{E}}}_m} - {\mu _0}{{\tilde{\textbf{H}}}_n} \cdot {{\tilde{\textbf{H}}}_m} + \omega _0^2/({{\varepsilon_\infty }\omega_p^2} ){{\tilde{\textbf{P}}}_n} \cdot {{\tilde{\textbf{P}}}_m} - 1/({{\varepsilon_\infty }\omega_p^2} ){{\tilde{\textbf{J}}}_n} \cdot {{\tilde{\textbf{J}}}_m}{d^3}\textbf{r} = 0}$ for $m \ne n$, can then be obtained by incorporating the auxiliary fields $\textbf{P}$ and $\textbf{J}$ into the Lorentz reciprocity theorem. With this orthogonality relationship, it is possible to project the scattered fields (${\textbf{E}_s}$, ${\textbf{H}_s}$, ${\textbf{P}_s}$, ${\textbf{J}_s}$) into the enlarged QNM vector (${\tilde{\textbf{E}}_m}$, ${\tilde{\textbf{H}}_m}$, ${\tilde{\textbf{P}}_m}$, ${\tilde{\textbf{J}}_m}$), which encompasses electromagnetic fields and material properties. This projection leads to the derivation of a closed-form expression for the QNM excitation coefficient [34].

For more complicated cases involving cross-action of Maxwell resonant field with other equations of physics (e.g., optomechanical cooling [102], lasing [103], photonic switching [104], electron nonlocality [105]) the physics of the system can no longer be fully taken into account with the medium response functions $\boldsymbol{\varepsilon }$ and $\boldsymbol{\mathrm{\mu}}$. The norms and orthogonality relations need to be rederived. In general, the derivation can be performed by strictly following the method in [32,34]. The most important thing one needs to care about is to fully account for the contributions of the extra current distributions, modified polarization fields, or boundary conditions, when using the Lorentz reciprocity theorem.

An illustrative example is found with QNMs incorporating non-classical response functions to accurately model electronic spill-out, nonlocality, or surface-enabled Landau damping for deep subwavelength confinements. This can be achieved thanks to an effective nonclassical surface polarization, ${\textbf{P}_{\textrm{surf}}} = \boldsymbol{\pi } + i{\mathrm{\omega }^{ - 1}}\textbf{K}$ with $\boldsymbol{\pi }$ being an out-of-plane dipole moment and $\textbf{K}$ an in-plane surface current density, treated with the Feibelman $d$-parameters formalism [106]. Interestingly, it has been found that, in addition to a volume integral term, the normalization formulas acquire a surface contribution arising from the nonclassical surface polarization [107]. Another important thing worth mentioning is that, for dispersive d parameters, the self-consistent Maxwell equation defines a non-linear eigenvalue problem, even when $\boldsymbol{\varepsilon }$ and $\boldsymbol{\mathrm{\mu}}$ are frequency independent. In analogy to the case of dispersive permittivities mentioned above, the Maxwell equations can be again mapped into a higher dimensional linear eigenvalue problem by introducing a surface auxiliary field, and thus an orthogonality relationship can be readily obtained [107].

4. Normalization of QNMs

It is interesting to go through the history of QNM normalization with the theoretical concepts of Section 3 in mind. Two distinct periods can be conveniently considered. In the 1990’s, P. T. Leung, K. Young, and coworkers at the Chinese University of Hong Kong played an important role in the development of QNM theory. They studied scalar fields in 1D cavities as well as vector fields in 3D systems with a spherical symmetry. Except for a few works that extended these results to 1D photonic bandgap structures [108,109], the domain remained dormant for 10 years. Then, the 2010’s witnessed a second active period where different groups addressed the issue of QNM normalization for arbitrary 3D geometries for which no analytical solution of Maxwell equations are available.

Figure 5 summarizes the main achievements of these two periods and places them along a timeline. Four different normalization methods have been elaborated in the 2010’s for 3D optical resonators made of dispersive materials. The methods have already been compared [40–44,74] and a few numerical tests have also been performed [40–42]. However, no clear consensus has been reached [40,42–44]. Some works consider that three integral expressions of the norm nearly provide the same result [21,40] and that a difference is noticeable only on rare occasions. On the contrary, some other works have shown that significant differences exist [41–43] and have argued that the normalization inherited from the work by the Hong-Kong group (referred to as LK normalization hereafter) is incorrect [42–43]. In addition, some arguments raised in [40,74] concerning the PML normalization are questionable [110].

 

Fig. 5. Milestones in the evolution of ideas on QNM normalization in electromagnetism. Two different periods can be identified. Normalization of scalar fields in 1D cavities and vector fields in 3D systems with a spherical symmetry have mainly been studied during the 1990’s. The 2010’s witnessed a second active period on QNM normalization for arbitrary 3D geometries. The emphasis was also placed during this period on the numerical implementation. The four normalization methods presented in Section 4 are distinguished by different colors: blue (LK norm), green (M norm), red (PML norm), and magenta (pole response). The adjectives “non-dispersive”, “dispersive”, “non-magnetic”, “magnetic” and “bi-anisotropic” refer to the material properties considered in the corresponding works.

Download Full Size | PDF

For all these reasons, we believe that it is important to compare in detail the different norms, so that future works can be built on firm ground. Hereafter, we start with a historical perspective presenting the evolution of ideas, then we review previous arguments and results exposed in [40,42–44,74], to further provide new mathematical demonstrations and numerical evidences. We hope that the whole clarifies the issue.

As a side note, we emphasize that normalization is not always required. Imagine that a QNM is not normalized and that its field is, say, three times larger than the field of its normalized version. However, compared to the normalized QNM, the unnormalized one will be excited three times less by any incident field. Overall, the scattered field will be the same. Therefore, as long as one is interested in computing only scattering matrices that relate the amplitudes of outgoing and ingoing waves, normalization is not required and unnormalized QNMs form a suitable basis to decompose the scattering spectra of nanoresonators [111,112]. However, QNM normalization has an inherent physical meaning: it teaches us the QNM ability to respond to a driving field and allows us to compare this ability for different QNMs of the same resonator or of different resonators [22]. When this driving field corresponds to the light emitted by a dipole source, the inverse of the normalized field squared represents the mode volume [22,32].

4.1 Historical perspective on QNM normalization

4.1.1 1990s: 1D and spherical cavities

The 1990s witnessed important development in QNM theory, essentially with the works from the Chinese University of Hong Kong [15,16,19,20,113]. The Hong Kong group derived important theorems on the orthogonality and the completeness of QNM expansions for scalar fields in 1D cavities. In [16], the authors regularized the QNMs by considering an analytic continuation for complex coordinates, exactly as discussed in Section 3. The regularization allowed them to define the norm of 1D QNMs and their orthogonality. In addition, the authors showed that QNMs form a complete set inside the cavity, essentially with the arguments developed in Section 2.3. However, the PML machinery was not yet invented at that time and the QNM regularization approach has not been further developed for geometries other than the 1D Fabry-Perot cavity.

During the course of their studies, the Hong Kong researchers have also developed another approach to normalize the QNMs of dielectric microspheres [15]. They presented their approach as an extension of the Zel’dovich 1D regularization formalism [8] to vector fields in 3D. However, no regularization is implemented, and the norm depends on a volume integral and a surface integral. We refer to this norm as the Leung-Kristensen (LK) norm, using the same terminology as in [42,43]. For dispersive but non-magnetic ($\mu = {\mu_0}$) media, it can be written as

Let us emphasize that the limit $\mathrm{\Omega } \to \infty $ was not present in the initial proposal [15] and that the authors were aware of a lack of rigor; they recognized in Appendix A that the norm is accurate only for QNMs with a large quality factor and for a given range of sizes of the domain $\mathrm{\Omega }$ surrounding the resonator.

The limit $\mathrm{\Omega } \to \infty $ has been introduced later on in [19], and was popularized in a series of works in the 2010s [40,74,114,115] that studied cavities with non-spherical shapes, e.g., photonic-crystal cavities. Taking the limit suggests that the sum of the volume and surface integrals, which depend on the size of $\mathrm{\Omega }$, becomes constant when the size of the integration domain $\mathrm{\Omega }$ becomes asymptotically large. In agreement with [42,43], we show in Sections 4.2.2 and 4.2.3 that this suggestion is fully incorrect.

To understand the 2010’s developments and debates, it is worth coming back on the two proposals of the Hong Kong group: the norm with the complex-coordinate regularization (integration along a path in the complex plane) and the LK norm with a surface integral term. For 1D resonators, the scalar field outside the cavity is simply given by an outgoing plane wave with a complex propagation constant. The integral in the complex plane starting at $x = R$ can thus be calculated analytically and it is simply proportional to $E_m^2(R )$, the QNM field squared at $x = R$, see the Section IV.A in [116]. Both norms are thus strictly equal for 1D cavities. However, in 2D and 3D, the field outside the cavity is no longer a plane wave and the equivalence between the two norms is broken. Therefore, the surface term used in [15,19,20] is incorrect for 3D systems (even spheres). This point was already mentioned as a tricky issue in Appendix A of [15] but, unfortunately, it has been forgotten in subsequent works.

4.1.2 2010s: towards arbitrary 3D geometries

The turn of the twenty-first century witnessed major progress of nanotechnologies and it became possible to fabricate photonic structures with increasingly complex geometries. Consistently, high demand prompted the development of QNM theories for 3D arbitrary geometries for which no analytical solution of Maxwell equations is known. The 2010s have been particularly active and fruitful with the contributions of different groups.

LK norm. Let us start by the work of P. T. Kristensen and colleagues who used the normalization proposed in the 1990s, initially for computing the mode volume of high-Q photonic-crystal cavities [114] and then plasmonic antennas [74,115].

M norm. A different norm was available since the beginning of the decade. In 2010, E. A. Muljarov, W. Langbein, and R. Zimmermann were developing a Brillouin-Wigner perturbation theory for open electromagnetic systems, which rapidly became known as the resonant-state expansion (RSE) [28]. They derived a QNM norm in the form of a sum of a volume integral and a surface integral, initially for non-dispersive and non-magnetic materials [28], which was later generalized for dispersive materials [42], and finally for magnetic (and even bi-anisotropic) materials [63].

In the following, we refer to this normalization method as the M norm. In the case of dispersive materials without bi-anisotropy, it can be written as [63]

Note that $\mathrm{\Omega }$ is usually taken as the minimal convex volume containing the resonator, but it can also be larger.

The derivation of Eqs. (18)–(19) relies on the divergence theorem, an analytic continuation of the QNM field to a function of complex frequencies, and a Taylor expansion close to the QNM eigenfrequency [29,42,63]. First-order spectral derivatives are then replaced by spatial derivatives of the QNM field under the assumption that the resonator is embedded in a uniform background. Note that slightly different expressions of the M norm have been proposed over the years. The initial expression in [28] contains first and second-order spatial derivatives of the electric field. The M norm can also be formulated with only first-order spatial derivatives of the electric field, as shown in Appendix B of [42]. Finally, the most general expression of the M norm given by Eqs. (18)–(19) contains first-order spatial derivatives, only, of both the electric and magnetic fields. Albeit technical, these last remarks are important to further understand the domain of validity, the ease of numerical implementation (see Appendix D), and the accuracy of the computed norms (see Section 4.2.4).

At first sight, the M norm, Eqs. (18)–(19), and the LK norm, Eqs. (16)–(17), appear similar: they both encompass a surface integral to compensate the divergence of the volume integral as the size of the integration domain $\mathrm{\Omega }$ is enlarged. However, the surface integrals are different. It has been recently shown that the surface term $I_{\textrm{surf}}^{\textrm{LK}}$ is incorrect, whereas the surface term $I_{\textrm{surf}}^\textrm{M}$ provides correct normalization for resonators embedded in a homogeneous background (no substrate) [42,43]. These key points are confirmed by simulations results and mathematical evidence in Section 4.2.

PML norm. In 2013, C. Sauvan and his colleagues proposed to work with PML-regularized QNMs [32]. As discussed in Section 3, the PML normalization is an incarnation of the generalized version (for vector fields in 3D) of the complex coordinate transform introduced in quantum mechanics for simple systems and scalar fields. We refer to this normalization method as the PML norm. The latter relies on the sum of two volume integrals [32]

The permittivity and permeability in Eq. (21) are those inside the PML area, see Section 3. Note that, in Eqs. (16), (18), and (20), we have used the same notation $\mathrm{\Omega }$ for denoting an unmapped domain around the resonator. For comparison purposes, these convex volumes can be chosen to be identical (see for instance Sections 4.2.2, 4.2.3, and 4.2.4). In practical situations, they can also be chosen to be different.

It is interesting to note that the first term of the PML norm, Eq. (20), is identical to the first term of the M norm, Eq. (18). Therefore, the PML norm and the M norm are identical if, and only if, ${I_{\textrm{PML}}} = I_{\textrm{surf}}^\textrm{M}$. In Section 4.2.4, we prove this equality for resonators embedded in a free space (no substrate).

Pole-response norm. Finally, let us introduce a fourth normalization method [3,117]. It is quite different from the three others since it does not rely on the calculation of integrals of the QNM field. In contrast, it relies on the fact that, for complex frequencies $\omega $ sufficiently close to the complex QNM eigenfrequency ${\tilde{\omega }_m}$, the field scattered by the resonator is proportional to the normalized QNM field with a known excitation coefficient: $\mathop {\textrm{lim}}\nolimits_{\omega \to {{\tilde{\omega }}_m}} ({\omega - {{\tilde{\omega }}_m}} ){\textbf{E}_s}({\textbf{r},\omega } )= \mathop {\lim }\nolimits_{\omega \to {{\tilde{\omega }}_m}} ({\omega - {{\tilde{\omega }}_m}} ){\alpha _m}(\omega ){\tilde{\textbf{E}}_m}(\textbf{r} )$.

For example, if the driving source is a point dipole $\textbf{p}$ located at ${\textbf{r}_0}$, the excitation coefficient is given by ${\alpha _m}(\omega )={-} \omega \textbf{p} \cdot {\tilde{\textbf{E}}_m}({{\textbf{r}_0}} )/({\omega - {{\tilde{\omega }}_m}} )$ [32]. One readily obtains the expression for the normalized QNM ${\tilde{\textbf{E}}_m}(\textbf{r} )$ from the knowledge of the scattered field ${\textbf{E}_s}({\textbf{r},\omega } )$

Note that the accuracy increases as $\omega $ is getting closer to ${\tilde{\omega }_m}$, and then decreases as numerical inaccuracies take place owing to the singularity of ${\textbf{E}_s}$. In practice, we keep away from this regime and excellent precision is obtained, see the documentation of the QNMPole solver in [38] and Fig. 6. As in [3], we refer to this method [117] as the pole-response normalization.

 

Fig. 6. Comparison between the PML normalization and the pole-response normalization for a 2D resonator on a metal substrate. (a) Sketch of the nanopatch resonator together with a definition of the integration domains $\mathrm{\Omega }$ and ${\mathrm{\Omega }_{\textrm{PML}}}$ used to normalize the mode with Eq. (20). (b) Distribution of the modulus of the normalized electric field of the lowest-order magnetic dipolar mode (${\tilde{\lambda }_m} = 725 + i54$ nm). The dashed lines represent the PML inner boundaries. Note the strong damping of the field inside the PMLs. (c) Relative error for the real parts of the PML and pole-response norms, $\textrm{Re}(\langle{\tilde{\boldsymbol{\mathrm{\phi}} }\textrm{|}{{\tilde{\boldsymbol{\mathrm{\phi}} }}\rangle_{\textrm{PML}}} - \langle\tilde{\boldsymbol{\mathrm{\phi}} }\textrm{|}{{\tilde{\boldsymbol{\mathrm{\phi}} }}\rangle_{\textrm{pole}}}} )/\textrm{Re}(\langle{\tilde{\boldsymbol{\mathrm{\phi}} }\textrm{|}{{\tilde{\boldsymbol{\mathrm{\phi}} }}\rangle_{\textrm{pole}}}} )$. The size R of the domain $\mathrm{\Omega }$ is normalized by $\lambda = \textrm{Re}({{{\tilde{\lambda }}_m}} )$. Adapted from [41].

Download Full Size | PDF

The pole-response normalization is a very general method. It led to the freeware QNMPole [118] launched in 2013. It can be used with any frequency-domain electromagnetic solver (using PMLs or not), with any dispersive materials (as long as analytical expressions of the frequency-dependent permittivity and permeability are known), and for “any” geometries including microcavities that leak into semi-infinite periodic waveguides. The latter situation, which can be often encountered in photonic-crystal integrated optics, was successfully tested in [99].

4.2 Exactness and domain of validity of the different norms

In this Section, four different approaches to normalize the QNMs are tested. The tests summarize available results from the literature and incorporate new ones to consolidate the analysis. They are accompanied by several theoretical analyses that highlight the differences between the norms, their domain of validity, or their mathematical equivalence. In particular, we demonstrate the equivalence between the M norm and the PML norm for the restrictive case of resonators in a free space.

4.2.1 Two methods are available for almost any geometries

Let us recall that the LK and M norms only apply to resonators in free space. They are not valid for resonators deposited on substrates (e.g., nanoparticles on a mirror [119,120], photonic crystal cavities [121]), embedded in thin films [122], or coupled to waveguides [99,123]. Thus, we start by considering the PML and pole-response normalizations, which are valid for the general case. Note, however, that the PML norm cannot be used for cavities coupled to periodic waveguides, in contrast with the pole-response normalization [99].

Figure 6 compares the PML normalization and the pole-response normalization for a resonator deposited on a metal substrate. The 2D geometry, depicted in Fig. 6(a), consists of a thin dielectric patch (refractive index 1.5, thickness 30 nm and length 100 nm) sandwiched between a gold substrate and a gold patch (thickness 30 nm and length 100 nm). We consider the lowest-order magnetic-dipolar mode for transverse magnetic (TM) polarization (magnetic field polarized along the z-direction). Its electric-field modulus is displayed in Fig. 6(b).

Maxwell equations are solved with the aperiodic Fourier modal method [124], by retaining 2001 harmonics in the Fourier expansion of the electromagnetic fields along the $x$-direction. For the evaluation of the PML norm, the PMLs are implemented with uniform layers of anisotropic and absorbing materials. The corresponding permittivity and permeability tensors (${\mu _z},{\varepsilon _x},{\varepsilon _y}$ in TM polarization) are equal to the permeability and permittivity of the medium in which the PML is added (air or metallic substrate, see Fig. 6) multiplied by the complex factor ${f_{\textrm{PML}}}$ or $1/{f_{\textrm{PML}}}$ depending on the tensor component [78]. The PML thickness is ${h_{\textrm{PML}}} = 2\textrm{Re}(\tilde{\lambda}_m)/3$ and ${f_{\textrm{PML}}} = 15 + 15i$; these parameters guaranty that the ratio between the modal field at the outer PML boundaries and the maximum field is smaller than ${10^{ - 5}}$. The same PML parameters are used for all values of R. We have checked that different parameters provide the same result (within the numerical accuracy) if the field is very weak at the outer PML boundaries.

Figure 6(c) shows the relative error between the PML and pole-response norms as the size R of the computational domain $\mathrm{\Omega }$ is increased. Since the pole-response norm does not depend on R, we have chosen to calculate it for the smallest computational domain for the sake of numerical accuracy. The PML norm is also independent of discretization, theoretically. However, as R is increased, the discretization grain decreases since the number of Fourier harmonics is constant, and numerical errors increase. An irrelevant tiny variation ($< $ 0.04%) of the norm, due to numerical inaccuracies, is observed as R is increased by a factor 30, from $R = 0.13\lambda $ to $4\lambda $. To our knowledge, this calculation is the only comparison between different normalization methods for the important practical case of resonators that are not in free space.

4.2.2 Stringent difference between the M and LK norms

In this Section, we focus on two normalization approaches (the LK and M norms) that are restricted to resonators in free space. They are not mathematically equivalent, as previously demonstrated with an analytical treatment for dielectric and metal spheres in air [42,43]. Here, using the PML norm as a refence, we show that they numerically result in substantially different norms, even for simple geometries (silver sphere in air), see Fig. 7.

 

Fig. 7. Comparison between the PML, LK and M norms for lowest-order electric dipolar mode of a silver sphere of radius 40 nm (${\tilde{\lambda }_m} = 390 + i32$ nm). (a) The PML norm is used as a reference value to benchmark the relative error on the real parts of the L norm, $\textrm{Re}(\langle{\tilde{\boldsymbol{\mathrm{\phi}} }\textrm{|}{{\tilde{\boldsymbol{\mathrm{\phi}} }}\rangle_{\textrm{LK}}} - \langle\tilde{\boldsymbol{\mathrm{\phi}} }\textrm{|}{{\tilde{\boldsymbol{\mathrm{\phi}} }}\rangle_{\textrm{PML}}}} )/\textrm{Re}(\langle{\tilde{\boldsymbol{\mathrm{\phi}} }\textrm{|}{{\tilde{\boldsymbol{\mathrm{\phi}} }}\rangle_{\textrm{PML}}}} )$ (solid blue curve), and the M norm, $\textrm{Re}(\langle{\tilde{\boldsymbol{\mathrm{\phi}} }\textrm{|}{{\tilde{\boldsymbol{\mathrm{\phi}} }}\rangle_\textrm{M}} - \langle\tilde{\boldsymbol{\mathrm{\phi}} }\textrm{|}{{\tilde{\boldsymbol{\mathrm{\phi}} }}\rangle_{\textrm{PML}}}} )/\textrm{Re}(\langle{\tilde{\boldsymbol{\mathrm{\phi}} }\textrm{|}{{\tilde{\boldsymbol{\mathrm{\phi}} }}\rangle_{\textrm{PML}}}} )$ (dashed red curve). The size R of the domain $\mathrm{\Omega }$ is normalized by $\lambda = \textrm{Re}({{{\tilde{\lambda }}_m}} )$. (b) provides an enlarged view of the dashed rectangle in (a). The sphere (radius 40 nm) is embedded in air. The permittivity of silver is given by a Drude model taken from [125]. The PML in spherical coordinates is implemented as a finite-size complex coordinate transform between R and $R + 4$ µm, $r \to R + ({1 + i/2} )({r - R} )$. For the M norm, the first-order derivatives in $I_{\textrm{surf}}^\textrm{M}$ are computed with a centered finite-difference scheme.

Download Full Size | PDF

We consider the electric dipolar QNM of lowest order supported by a silver sphere of radius 40 nm. We use an analytical expression of the QNM fields [29]. The volume and surface integrals required to compute the PML, LK and M norms are computed numerically with the Gaussian quadrature rule. The integration domain $\mathrm{\Omega }$ is a sphere of radius R and the calculation is repeated for increasing values of R. Since the QNM field takes a simple analytical expression in spherical coordinates, the PML is implemented as a radial complex mapping inside a 4-µm-thick shell surrounding the unmapped domain $\mathrm{\Omega }$ defined by $r < R$, $r \to R + ({1 + i/2} )({r - R} )$ for $R < r < R + 4 \mathrm{\mu}\textrm{m}$.

As is increased from 60 nm to 4 µm, the PML norm remains constant with an excellent accuracy (variation smaller than ${510^{ - 13}}$). This value is further used as a reference to compare the performance of the M and LK norms. Figure 7 displays the relative errors (in $\%$) on the real parts of the LK norm, $\textrm{Re}(\langle{\tilde{\boldsymbol{\mathrm{\phi}} }\textrm{|}{{\tilde{\boldsymbol{\mathrm{\phi}} }}\rangle_{\textrm{LK}}} - \langle\tilde{\boldsymbol{\mathrm{\phi}} }\textrm{|}{{\tilde{\boldsymbol{\mathrm{\phi}} }}\rangle_{\textrm{PML}}}} )/\textrm{Re}(\langle{\tilde{\boldsymbol{\mathrm{\phi}} }\textrm{|}{{\tilde{\boldsymbol{\mathrm{\phi}} }}\rangle_{\textrm{PML}}}} )$ (solid blue curve), and M norm, $\textrm{Re}(\langle{\tilde{\boldsymbol{\mathrm{\phi}} }\textrm{|}{{\tilde{\boldsymbol{\mathrm{\phi}} }}\rangle_\textrm{M}} - \langle\tilde{\boldsymbol{\mathrm{\phi}} }\textrm{|}{{\tilde{\boldsymbol{\mathrm{\phi}} }}\rangle_{\textrm{PML}}}} )/\textrm{Re}(\langle{\tilde{\boldsymbol{\mathrm{\phi}} }\textrm{|}{{\tilde{\boldsymbol{\mathrm{\phi}} }}\rangle_{\textrm{PML}}}} )$ (dashed red curve).

Two conclusions can be drawn: (1) the LK norm is not independent of the size R of the calculation domain; it oscillates around the value of the PML norm for $R > \lambda /2$, with $\lambda = \textrm{Re}({{{\tilde{\lambda }}_m}} )$, and the oscillations exponentially grow up as R increases. The relative error is quite large even for small domains: it is $> 1\%$ for $R < \lambda /2$ and even $> 10\%$ for $R < \lambda /5$. The reason for the divergence of the LK norm is analyzed in Section 4.2.3 and in [42,43]. (2) The M and PML norms provide almost identical results; the relative difference is $< {10^{ - 4}}\%$. In Section 4.2.4, we mathematically prove that the M and PML norms are equivalent for resonators embedded in a free space. We thus conclude that the slight difference is due to difficulties encountered for the numerical implementation of the M norm, see Appendix D.

4.2.3 Divergence of the LK norm

The exponentially divergent oscillations of the LK norm in Fig. 7 does not come as a surprise [40,42,43]. A similar behavior can be found in Fig. 3 in [40]. In fact, the LK norm can be derived analytically for the QNMs of spheres in free space. Equations (24)–(26) in [40] and Eq. (C(1)5) in [42] show that the LK norm diverges as $\textrm{exp}({\textrm{Re}({\tilde{\omega }/c} )R/Q} )$, Q being the quality factor of the mode.

We now show that the LK norm corresponds, in fact, to the first term of the PML norm in Eq. (20). This implies that the LK norm completely misses the regularization contribution due to the volume integral in the PML. The surface term $I_{\textrm{surf}}^{\textrm{LK}}$ in Eq. (16) is actually a decoy: it does not compensate the divergence of the volume integral over the domain $\mathrm{\Omega }$, but simply accounts for the magnetic-field contribution to the volume integral.

By using the divergence theorem and radiation conditions on the border $\mathrm{\Sigma }$ of a large spherical domain $\mathrm{\Omega }$, we show in Appendix E that

The ${\approx} $ sign simply comes from the use of radiation conditions, which are not rigorous for a finite-size domain $\mathrm{\Omega }$. Then injecting Eq. (23) into Eq. (16), it follows that

This is exactly the expression of the first term of the PML norm in Eq. (20) for non-dispersive magnetic materials. Equation (24) clearly shows that $\langle{\tilde{\boldsymbol{\mathrm{\phi}} }_m}\textrm{|}{\tilde{\boldsymbol{\mathrm{\phi}} }_m}\rangle_{\textrm{LK}}$ is actually an integral of the QNM fields over the unmapped domain $\mathrm{\Omega }$ only. The absence of any regularization simply explains why the LK norm exponentially grows as $R \to \infty $.

We believe that it is important to be aware that the LK norm, intrinsically, does not contain any regularization. This conclusion removes and clarifies a long-standing mistake of the early 1990s which has been passed on until recently [15,19,74,114,115].

4.2.4 Mathematical equivalence of the PML and M norms for resonators in free space

Taking the same volume $\mathrm{\Omega }$, the first terms of the PML and M norms are identical, see Eqs. (18) and (20). Thus, the two norms are equivalent if, and only if, ${I_{\textrm{PML}}}$ in Eq. (21) and $I_{\textrm{surf}}^\textrm{M}$ in Eq. (19) are also identical. Hereafter, we show that ${I_{\textrm{PML}}} = I_{\textrm{surf}}^\textrm{M}$ for resonators in a free space (no substrate).

For the demonstration, we follow the same approach as in [42,63]. We consider the field radiated by a source emitting at a frequency $\omega$ that tends towards the QNM frequency ${\tilde{\omega }_m}$ and further use the Lorentz reciprocity theorem to transform the volume integral in ${I_{\textrm{PML}}}$ into the surface integral in $I_{\textrm{surf}}^\textrm{M}$.

Let us consider the electromagnetic field $\boldsymbol{\mathrm{\phi}} = ({\textbf{E},\textbf{H}} )$ solution of Maxwell equations at the frequency $\omega$ for an electric current source $\textbf{j}({\textbf{r},\omega} )= ({\omega- {{\tilde{\omega }}_m}} )\textbf{S}(\textbf{r} )$ located in the domain $\mathrm{\Omega }$. Because of the frequency dependence of the current, $({\textbf{E},\textbf{H}} )\to ({{{\tilde{\textbf{E}}}_m},{{\tilde{\textbf{H}}}_m}} )$ as $\omega\to {\tilde{\omega }_m}$. Note that the mode ${\tilde{\boldsymbol{\mathrm{\phi}} }_m}$ is not necessarily normalized.

We apply the Lorentz reciprocity theorem (see Appendix A) in the domain ${\mathrm{\Omega }_{\textrm{PML}}}$ to the field $\boldsymbol{\mathrm{\phi}} $ at the frequency $\omega$ and to the mode ${\tilde{\boldsymbol{\mathrm{\phi}} }_m}$ at ${\tilde{\omega }_m}$. Since it is located inside $\mathrm{\Omega }$, the source $\textbf{j}({\textbf{r},\omega} )$ is located outside ${\mathrm{\Omega }_{\textrm{PML}}}$, the source term in Eq. (31) cancels out. Additionally, since ${\tilde{\boldsymbol{\mathrm{\phi}} }_m}$ and $\boldsymbol{\mathrm{\phi}} $ vanish at the outer border of ${\mathrm{\Omega }_{\textrm{PML}}}$ for ideal PMLs, the surface term reduces to an integral over the closed surface $\mathrm{\Sigma }$. After dividing each term by $({\omega- {{\tilde{\omega }}_m}} )$, we obtain

Finally, taking the limit $\omega\to {\tilde{\omega }_m}$ in Eq. (25) provides

So far, we have not made any assumption on the medium surrounding the resonator and Eq. (27) shows that it is possible to replace the volume integral over the PML domain ${\mathrm{\Omega }_{\textrm{PML}}}$ in Eq. (20) by a surface integral over its inner boundary $\mathrm{\Sigma }$. However, without any further assumption that we detail hereafter, the right-hand term, and thus the M norm, does not only rely on the QNM itself but additionally encompasses the derivatives with respect to the frequency, $\frac{{\partial \textbf{E}}}{{\partial \omega}}$ and $\frac{{\partial \textbf{H}}}{{\partial \omega}}$, of the analytic continuation of the QNM.

This issue can be solved if we assume that the resonator is surrounded by a homogeneous medium. In that case, the fields over the surface $\mathrm{\Sigma }$ (i.e., outside the resonator and at the inner PML boundary) consist of only outgoing waves; they are functions of the variable $\omega\textbf{r}$ and the spectral derivatives can be rewritten as [42]

In conclusion, for resonators embedded in a homogeneous medium, ${I_{\textrm{PML}}} = I_{\textrm{surf}}^\textrm{M}$ and the PML and M normalizations, published respectively in [32] and in [63], are strictly equivalent. Table 2 summarises the domain of validity of the available normalization methods.

Table 2. Summary of normalization methods used in the 2010s for reciprocal materials.

View Table | View all tables in this article

5. Conclusion

We have reviewed key concepts attached to the quasinormal modes of open electromagnetic systems: normalization, orthogonality, and completeness. Normalization allows us to directly quantify the strength of the interaction of light with a resonance. Orthogonality enables the projection of any field on the QNM basis. Completeness guaranties the convergence and exactness of QNM expansions of the scattered field.

The normalization of QNMs in electromagnetism has not been achieved straightforwardly and an error or an approximation introduced in the 1990s has been unfortunately repeated and amplified in the 2010s, leading to some confusion in the literature. The difficulty is due to the divergence of the QNM field in the open space. We again emphasize that the divergence is not unphysical. The reverse is actually true. We have clarified the confusion, putting forth the errors made and the domain of validity of the available correct methods.

Nowadays the normalization is no longer a theoretical issue. Methods and even free software to normalize them for (virtually) any system have been available for almost ten years. Perhaps, a relevant question is to derive numerical methods that compute normalized QNMs in a very effective and accurate way, although the existing general methods already appear quite effective, see Section 4.2.1.

We have recalled a worthwhile result, quite well known for a long time [9,16,26,27], on the completeness of expansions based on true QNMs. Completeness is restrictively warranted inside resonators surrounded by free space (no substrate). We have highlighted a natural and quite general way to complement the incomplete QNM basis to achieve completeness: regularization [25]. We have described recent results in electromagnetism, aiming at setting up effective formalisms and numerical tools based on complex coordinate mappings, viewed as PMLs. We have presented data to show that an extended basis composed of regularized QNMs and additional numerical modes forms a complete basis of the regularized Hilbert space. Interestingly, we may use orthonormality relations that are quite similar to those used in “Hermitian electrodynamics” [1] and many geometries of current interest in plasmonics and nanophotonics can be faithfully analyzed.

After 30 years or more, research on QNMs in electromagnetism has continued unabated. This reflects the underlying importance of resonances in the contemporary physics and applications in photonics. We expect that this review will contribute to contextualizing and clarifying the achievements of the last decade.

Appendix A. Unconjugated form of Lorentz reciprocity theorem for absorbing and dispersive media

Consider a system characterized by the position- and frequency-dependent permittivity and permeability tensors $\boldsymbol{\varepsilon}({\textbf{r},\omega} )$ and $\boldsymbol{\mu}({\textbf{r},\omega})$. We assume reciprocal materials, $\boldsymbol{\mu} = {\boldsymbol{\mu}^T}$ and $\boldsymbol{\varepsilon} = {\boldsymbol{\varepsilon}^T}$, where the superscript T denotes matrix transposition. The electromagnetic field $({\textbf{E},\textbf{H}} )$ radiated by a current source density $\textbf{j}$ at the frequency $\omega$ is given by the two curl Maxwell equations,

(30)$$\nabla \times \textbf{E} = i\omega\boldsymbol{\mu}({\boldsymbol{r},\omega} )\textbf{H}{\; \; \; \; \; }\textrm{and}{\; \; \; \; \; }\nabla \times \textbf{H} ={-} i\omega\boldsymbol{\varepsilon}({\textbf{r},\omega} )\textbf{E} + \textbf{j}\; .$$

Lorentz reciprocity theorem relates two different solutions of Maxwell equations, $({{\textbf{E}_1},{\textbf{H}_1},{\omega_1},{\textbf{j}_1}} )$ and $({{\textbf{E}_2},{\textbf{H}_2},{\omega_2},{\textbf{j}_2}} )$ labelled by the index 1 and 2. It is derived by applying the divergence theorem to the vector ${\textbf{E}_2} \times {\textbf{H}_1} - {\textbf{E}_1} \times {\textbf{H}_2}$ and by using Eq. (30),

(31)$$\begin{aligned}&\mathop {\iint}\nolimits_\mathrm{\Sigma}^{} ({{\tilde{\textbf{E}}_2}\times{\tilde{\textbf{H}}_1} - {\tilde{\textbf{E}}_1}\times{\tilde{\textbf{H}}_2}} )\cdot \textbf{n}dS \\ &= i\iiint_\mathrm{\Omega } {\{{{\textbf{E}_1} \cdot [{{\mathrm{\omega }_1}\boldsymbol{\varepsilon}({{\mathrm{\omega }_1}} )- {\mathrm{\omega }_2}\boldsymbol{\varepsilon}({{\mathrm{\omega }_2}} )} ]{\textbf{E}_2} - {\textbf{H}_1} \cdot [{{\mathrm{\omega }_1}\boldsymbol{\mu}({{\mathrm{\omega }_1}} )- {\mathrm{\omega }_2}\boldsymbol{\mu}({{\mathrm{\omega }_2}} )} ]{\textbf{H}_2}} \}dV} \\ &- \iiint_\mathrm{\Omega } {({{\textbf{j}_1} \cdot {\textbf{E}_2} - {\textbf{j}_2} \cdot {\textbf{E}_1}} )dV}, \end{aligned}$$

where $\mathrm{\Sigma }$ is an arbitrary closed surface defining a volume $\mathrm{\Omega }$ and $\textbf{n}$ is the unit vector normal to the surface $\mathrm{\Sigma }$ (oriented outward).

Appendix B. Direct proof of the PML-norm invariance with the PML parameters

The following demonstration is very general and quite straightforward. Let us consider a resonator composed of non-magnetic materials and illuminated by a source at frequency $\omega $. In the scattered field formulation, the scattered field [${\textbf{E}_s},{\textbf{H}_s}]$ satisfies the Maxwell equations with a current distribution inside the resonator $\textbf{J}(\omega )={-} i\omega \mathrm{\Delta }\boldsymbol{\varepsilon }({\textbf{r},\omega } ){\textbf{E}_b}({\omega ,\boldsymbol{r}} )$, where $\mathrm{\Delta }\boldsymbol{\varepsilon }({\textbf{r},\omega } )$ is null everywhere except in the resonator and ${\textbf{E}_b}({\omega ,\boldsymbol{r}} )$ is the driving field. All details with the same notations are found in Annex 2 in [3].

We consider a first regularized version of the scattering problem, in which a PML is introduced outside the resonator that necessarily has a finite spatial extend ($\mathrm{\Delta }\boldsymbol{\varepsilon }$ is null outside a compact support ${\Omega _{\textrm{res}}}$). Referring to Fig. 2, this defines an unmapped compact domain bounded by the inner boundary of the PML domain, a rotational angle $\theta $ that satisfies Eq. (12) and a stretching coefficient if any.

We assume that the ${m^{th}}$ QNM is revealed by the PML and denote its regularized field by [${\tilde{\textbf{E}}_m}$,..]. By applying the Lorentz reciprocity theorem, see Appendix A, we find

(32)$$\begin{aligned} &\int_{\mathrm{\Omega } \cup {\mathrm{\Omega }_{\textrm{PML}}}} {[{{\textbf{E}_s} \cdot ({\omega \boldsymbol{\varepsilon }(\omega )- {{\tilde{\omega }}_m}\boldsymbol{\varepsilon }({{{\tilde{\omega }}_m}} )} ){{\tilde{\textbf{E}}}_m} - {\textbf{H}_s} \cdot({\omega \boldsymbol{\mathrm{\mu}}(\omega )- {{\tilde{\omega }}_m}\boldsymbol{\mathrm{\mu}}({{{\tilde{\omega }}_m}} )} ){{\tilde{\textbf{H}}}_m}} ]dV} \\ \qquad &=- i\int_\mathrm{\Omega } {\textbf{J} \cdot{{\tilde{\textbf{E}}}_m}dV,}\end{aligned}$$

where the integral is performed over the entire space, including the unmapped compact domain $\mathrm{\Omega } \supseteq {\mathrm{\Omega }_{\textrm{res}}}$ and the PML domain ${\mathrm{\Omega }_{\textrm{PML}}}$. Note that the surface-integral term vanishes and does not show up in Eq. (32) because all fields, ${\textbf{E}_s}$, ${\tilde{\textbf{E}}_m}$,… exponentially vanish for $|\textbf{r} |\to \infty $. QNMs are poles and thus, the scattered field diverges for any as the source frequency $\omega \to {\tilde{\omega }_m}$. For a pole of order 1, the residue can be written as

(33)$$\mathop {\lim }\limits_{\omega \to {{\tilde{\omega }}_m}} ({\omega - {{\tilde{\omega }}_m}} )[{{\textbf{E}_s}({\textbf{r},\omega } ),{\textbf{H}_s}({\textbf{r},\omega } )} ]= {\gamma _m}[{{{\tilde{\textbf{E}}}_m}(\textbf{r} ),{{\tilde{\textbf{H}}}_m}(\textbf{r} )} ].$$

Similar considerations apply to poles of higher orders. ${\gamma _m}$ is a complex coefficient that may be simply determined by evaluating the limit of Eq. (32) as $\omega \to {\tilde{\omega }_m}$. By using Eq. (33) and the definition of the derivative at $\omega = {\tilde{\omega }_m}$, $\mathop {\lim }\nolimits_{\omega \to {{\tilde{\omega }}_m}} [{\omega \boldsymbol{\varepsilon }(\omega )- {{\tilde{\omega }}_m}\boldsymbol{\varepsilon }({{{\tilde{\omega }}_m}} )} ]/({\omega - {{\tilde{\omega }}_m}} )\equiv \frac{{\partial \omega \boldsymbol{\varepsilon }}}{{\partial \omega }}$, we obtain [32]

(34)$${\gamma _m} = \frac{{ - i\mathop \smallint \nolimits_{{\mathrm{\Omega }_{\textrm{res}}}}^{} \textbf{J}({{{\tilde{\omega }}_m}} )\cdot {{\tilde{\textbf{E}}}_m}dV}}{{\mathop \smallint \nolimits_{\mathrm{\Omega } \cup {\mathrm{\Omega }_{\textrm{PML}}}}^{} \left[ {{{\tilde{\textbf{E}}}_m} \cdot \frac{{\partial \omega \boldsymbol{\varepsilon }}}{{\partial \omega }}{{\tilde{\textbf{E}}}_m} - {\mu_0}{{\tilde{\textbf{H}}}_m} \cdot {{\tilde{\textbf{H}}}_m}} \right]dV}}\; .$$

Note that the source term is given by $\textbf{J}({{{\tilde{\omega }}_m}} )={-} i{\tilde{\omega }_m}\mathrm{\Delta }\boldsymbol{\varepsilon }({\textbf{r},{{\tilde{\omega }}_m}} ){\textbf{E}_b}({{{\tilde{\omega }}_m},\boldsymbol{r}} )$.

Now let us consider another regularization performed with a new PML that is different from the previous one, possibly because the new unmapped compact domain $\mathrm{\Omega}^{\prime}$ or the new rotational angle are different. To demonstrate the invariance of the PML norm, we assume that the QNM [${\tilde{\textbf{E}}_m}$, ${\tilde{\textbf{H}}_m}$] is again revealed by the new PML, i.e., the new rotational angle $\theta^{\prime}$ satisfies Eq. (12). The new QNM, denoted [$\tilde{\textbf{E}}^{\prime}_m$, $\tilde{\textbf{H}}^{\prime}_m$], is identical to the previous one in the unmapped space: [$\tilde{\textbf{E}}^{\prime}_m$, $\tilde{\textbf{H}}^{\prime}_m$] $ = $ [${\tilde{\textbf{E}}_m}$, ${\tilde{\textbf{H}}_m}$] for $\textbf{r} \in \mathrm{\Omega } \cap \mathrm{\Omega}^{\prime}$. It differs from the previous one only in the PML domain. Owing to the fact that Maxwell equations have unique solutions, $\mathop {\lim }\nolimits_{\omega \to {{\tilde{\omega }}_m}} ({\omega - {{\tilde{\omega }}_m}} )[{{\textbf{E}_s}({\textbf{r},\omega } ),{\textbf{H}_s}({\textbf{r},\omega } )} ]= \mathop {\lim }\nolimits_{\omega \to {{\tilde{\omega }}_m}} ({\omega - {{\tilde{\omega }}_m}} )[{\textbf{E}^{\prime}_s({\textbf{r},\omega } ),\textbf{H}^{\prime}_s({\textbf{r},\omega } )} ]\;\forall \;\textbf{r} \in \mathrm{\Omega } \cap \mathrm{\Omega} ^{\prime}$. Therefore, the ${\gamma _m}$ coefficient is unchanged and, since the numerator in Eq. (34) is independent of the PML, the denominator, i.e., the PML norm, is shown to be also independent of the PML. QED.

Appendix C. Normalization for non-reciprocal materials

In this Appendix, we focus on the case where the reciprocity of the system is broken, namely, $\boldsymbol{\mathrm{\mu}} \ne {\boldsymbol{\mathrm{\mu}}^\textrm{T}}$ or $\boldsymbol{\varepsilon } \ne {\boldsymbol{\varepsilon }^\textrm{T}}$ and derive the orthonormality condition for regularized QNMs.

Non-reciprocal materials are usually found in magneto-optical materials, such as yttrium iron garnet (YIG), iron garnets (SrFe₁₂O₁₉), biased by an external or static magnetic field to break the time-reversal symmetry [127]. Optical resonators made of non-reciprocal materials are essential for many applications, such as one-way propagation, optical isolating, broadband field localization, and so on. To provide an intuitive understanding of these systems, it would be of interest to develop a method to normalize the QNMs when the reciprocity is broken.

The norm and orthogonality relationship for the QNMs of non-reciprocal resonators can be easily generalized from those of the reciprocal systems mentioned above. The generalization relies on the introduction of a left QNM, $\tilde{\textbf{E}}_m^{(L)}$, which can be found by solving the source-free Maxwell operator with transposed permittivity and permeability, ${\boldsymbol{\varepsilon}^\textrm{T}}$ and ${\boldsymbol{\mathrm{\mu}}^\textrm{T}}$ [95]. Plugging the left QNM and another QNM of the concerned system $\tilde{\textbf{E}}_n^{(R)}$ (i.e., the mode of the system with $\boldsymbol{\varepsilon }$ and $\boldsymbol{\mathrm{\mu}}$, we refer as the right QNM) into the Lorentz reciprocity formula, and noticing that they satisfy the source-free Maxwell equations, we obtain the relationship

(35)$$ \iiint_{\Omega \cup \Omega_{\mathrm{PML}}} \tilde{\mathbf{H}}_{m}^{(L)} \cdot\left[\tilde{\omega}_{n} \boldsymbol{\mu}\left(\tilde{\omega}_{n}\right)-\tilde{\omega}_{m} \boldsymbol{\mu}\left(\tilde{\omega}_{m}\right)\right] \tilde{\mathbf{H}}_{n}^{(R)}-\tilde{\mathbf{E}}_{m}^{(L)} \cdot\left[\tilde{\omega}_{n} \boldsymbol{\varepsilon}\left(\tilde{\omega}_{n}\right)-\tilde{\omega}_{m} \boldsymbol{\varepsilon}\left(\tilde{\omega}_{m}\right)\right] \tilde{\mathbf{E}}_{n}^{(R)} d V=0. $$

Note that, for nonreciprocal systems, $\tilde{\textbf{H}}_n^{(L )} \ne \tilde{\textbf{H}}_n^{(R )}$ and $\tilde{\textbf{E}}_n^{(L )} \ne \tilde{\textbf{E}}_n^{(R )}$, but they share the same eigenvalues, that is $\tilde{\omega }_n^{(R )} = \tilde{\omega }_n^{(L )} = {\tilde{\omega }_n}$. The volume integral is performed over the whole space including the PML regions.

Further applying the Lorentz reciprocity to the left QNM and the scattered field of the resonator when it is driven by a dipole source emitting close to the resonance frequency, strictly following the approach in [32] and using Eq. (35), the QNM norm for non-reciprocal systems is derived

(36)$$ Q N_{n}=\iiint_{\Omega \cup \Omega_{\mathrm{PML}}} \tilde{\mathbf{E}}_{n}^{(L)} \cdot \frac{\partial(\omega \boldsymbol{\varepsilon})}{\partial \omega} \tilde{\mathbf{E}}_{n}^{(R)}-\tilde{\mathbf{H}}_{n}^{(L)} \cdot \frac{\partial(\omega \mu)}{\partial \omega} \tilde{\mathbf{H}}_{n}^{(R)} d V. $$

For compactness, the detailed derivation of this formula is not shown here, and we refer the interested reader to Section 1.3.3 of [128].

The introduction of left QNMs in the normalization and orthogonality relations impacts the definition of some important physical quantities, compared to reciprocal cases. For example, the mode volume for non-reciprocal systems takes the following form [22]

(37)$${\tilde{V}_n}({\textbf{r},\textbf{u}} )= {[{2{\varepsilon_0}({\tilde{\textbf{E}}_n^{(L )} \cdot \textbf{u}} )({\tilde{\textbf{E}}_n^{(R )} \cdot \textbf{u}} )/Q{N_n}} ]^{ - 1}},$$

where $\textbf{u}$ is the unit vector defining the polarization direction.

Figure 8 shows the field and mode volume maps of a QNM of a 2D yttrium-iron garnet wire in a homogenous air background. Under an external dc magnetic field along the z direction, the wire exhibits a strong gyromagnetic anisotropy, for a relative permeability tensor taking the form $\boldsymbol{\mathrm{\mu}} = [{{\mu_r}, - i{k_r},0;\;i{k_r},{\mu_r},0;0,0,{\mu_\infty }} ]$. In contrast to the reciprocal resonators, the fields for the left and right QNMs are different, as can be seen from Fig. 8(a). Therefore, to compute the mode volume in Fig. 8(b), two independent computations are required.

 

Fig. 8. (a) QNM field maps of a 2D yttrium-iron garnet wire in an air background. The first and second panels correspond to the maps for the real part of the electric field $z$-component for the left ($\tilde{\textbf{E}}_n^{(L )}$) and right ($\tilde{\textbf{E}}_n^{(R )}$) QNMs, respectively. The maps are computed for a QNM with its electric field parallel to the z direction and an eigenfrequency of $7.90 + 0.04i$ GHz. The radius of the rod is $9.1$ mm and it has an isotropic relative permittivity of $\varepsilon = 15$. (b) The map of $\textrm{Re}({\tilde{V}_n^{ - 1}} )$ of the QNM with $\textbf{u} = {\textbf{e}_z}$. Note that because here we consider a 2D resonator, the mode volume takes the unit of ${m^2}$ instead of ${m^3}$ in 3D cases. This COMSOL model is freely available in the freeware package MAN [38].

Download Full Size | PDF

Appendix D. Discussion on the numerical implementation of the M norm

As discussed in Section 4.1.2, several equivalent expressions of the M norm have been published over the years. The first expression involves only the electric field, with first and second-order spatial derivatives [28,42,47]. The M norm can also be formulated with only first-order spatial derivatives of the electric field, as shown in Appendix B of [42]. Finally, a more general expression valid for chiral and bi-anisotropic materials requires only the calculation of first-order spatial derivatives of the electric and the magnetic field, see Eq. (18) or [63]. In the case of a sphere, the electromagnetic field is known analytically, and the spatial derivatives can be exactly calculated. For a different geometry, the QNM field, and thus its derivatives, must be calculated numerically. The numerical evaluation of the derivatives introduces a small numerical error in the surface integral. As a result, the divergence of the volume integral is not perfectly compensated, and a small error is introduced in the norm. The magnitude of the error depends on the implementation of the M norm. This issue has been investigated in Appendix G of [42].

Figure 9 displays a comparison between both implementations of the M norm if the derivatives are calculated numerically with a centered finite-difference scheme. The relative error between the PML norm and the M norm oscillates around zero with an amplitude that increases exponentially with R. An implementation of the M norm with only first-order derivatives [Eq. (19), dashed red curve] reduces the error by three orders of magnitude compared to the implementation with second-order derivatives (solid blue curve), see Fig. 9. The error could be further reduced by using a different numerical method to calculate the spatial derivative. However, methods more accurate than the centered finite-difference scheme usually increase the computational load since they require a more precise knowledge of the field behavior in the vicinity of the derivation point, i.e., multiple (larger than 2) numerical calculations of the field [129]. Note that, for both implementations, the size of the domain $\mathrm{\Omega }$ that provides the smallest error is in the range $\lambda < R < 2\lambda $, with $\lambda = \textrm{Re}({{{\tilde{\lambda }}_m}} )$. The use of smaller integration domains results in larger errors.

 

Fig. 9. Comparison between two different implementations of the M normalization for the same sphere as in Fig. 7. The spatial derivatives in the surface term $I_{\textrm{surf}}^\textrm{M}$ are calculated with a centered finite-difference scheme. We plot the same relative error as in Fig. 7. The solid blue curve shows the result of an implementation with second-order derivatives and electric field only (Eq. (4) in [47]) and the dashed red curve shows the result of an implementation with only first-order derivatives of the electric and magnetic fields contributing on equal footing, see Eqs. (18)–(19) and [63]. The latter reduces the relative error by three orders of magnitude. (b) provides an enlarged view of (a). The insets display a zoom of the areas inside the rectangular boxes. Note that the dashed red curve is the same as in Fig. 7.

Download Full Size | PDF

Appendix E. Demonstration of Eq. (23)

In this Appendix, we demonstrate Eq. (23), i.e., a relation between the surface integral $I_{\textrm{surf}}^{\textrm{LK}}$ appearing in the LK norm and a volume integral of the electric and magnetic fields over a spherical domain $\mathrm{\Omega }$. The expression of the surface integral $I_{\textrm{surf}}^{\textrm{LK}}$ is given by Eq. (17).

Let us consider a QNM $({\tilde{\textbf{E}}_m},\;{\tilde{\textbf{H}}_m},{\tilde{\omega }_m})$ solution of Maxwell equations in the absence of source. We first apply the divergence theorem to the vector ${\tilde{\textbf{E}}_m} \times {\tilde{\textbf{H}}_m}$ and the closed spherical surface $\mathrm{\Sigma }$ defining the volume $\mathrm{\Omega }$. Simple manipulations similar to the derivation of Lorentz reciprocity theorem in Appendix A lead to

(38)$$\mathop {\iint }\nolimits_\mathrm{\Sigma }^{} {\tilde{\textbf{E}}_m} \times {\tilde{\textbf{H}}_m} \cdot \textbf{u}dS = i{\tilde{\omega }_m}\iiint_\mathrm{\Omega } {({{{\tilde{\textbf{E}}}_m} \cdot \boldsymbol{\varepsilon }{{\tilde{\textbf{E}}}_m} + {{\tilde{\textbf{H}}}_m} \cdot \boldsymbol{\mathrm{\mu}}{{\tilde{\textbf{H}}}_m}} ){d^3}\textbf{r}\; ,} $$

with $\textbf{u}$ is the unit vector normal to the sphere $\mathrm{\Sigma }$ (oriented outward).

For a resonator surrounded by a homogeneous medium and a large domain $\mathrm{\Omega }$, one can use radiation conditions on the surface $\mathrm{\Sigma }$ that provide a simple relation between the electric and magnetic fields. In the far field, the electromagnetic field locally behaves as an outgoing plane wave propagating along $\textbf{u}$. Therefore, for positions $\textbf{r}$ on the surface $\mathrm{\Sigma }$, one can write

(39)$${\tilde{\textbf{H}}_m} \approx {\varepsilon _0}nc\textbf{u} \times {\tilde{\textbf{E}}_m}\; ,$$

where n is the refractive index of the surrounding medium and c is the speed of light. The ${\approx} $ sign simply comes from the fact that Eq. (39) is not rigorous for a finite-size domain $\mathrm{\Omega }$. It is only asymptotically valid as $r \to \infty $.

Inserting Eq. (39) into the left-hand side of Eq. (38) and using the fact that the electric field of a plane wave is transverse to its propagation direction leads to

(40)$$\mathop {\iint }\nolimits_\mathrm{\Sigma }^{} {\tilde{\textbf{E}}_m} \times {\tilde{\textbf{H}}_m} \cdot \textbf{u}dS \approx {\varepsilon _0}nc\mathop {\iint}\nolimits_\mathrm{\Sigma }^{} \tilde{\textbf{E}}_m^2dS.$$

According to the definition of the surface term $I_{\textrm{surf}}^{\textrm{LK}}$ in Eq. (17), Eq. (38) simply becomes

(41)$$I_{\textrm{surf}}^{\textrm{LK}} \approx{-} \iiint_\mathrm{\Omega } {({{{\tilde{\textbf{E}}}_m} \cdot \boldsymbol{\varepsilon }{{\tilde{\textbf{E}}}_m} + {{\tilde{\textbf{H}}}_m} \cdot \boldsymbol{\mathrm{\mu}}{{\tilde{\textbf{H}}}_m}} ){d^3}\textbf{r}\; .} $$

QED.

Acknowledgements

CS and PL thank Jean-Paul Hugonin for stimulating discussions and continuous help. PL and TW would like to thank Wei Yan. PL thanks Thomas Christopoulos, Odysseas Tsilipakos, Aristeidis Karalis, and Mohammed Benzaouia for helpful comments.

Disclosures

The authors declare no conflicts of interest.

Data availability

Data underlying the results presented in this paper are not publicly available at this time but may be obtained from the authors upon reasonable request.

References

1. R. F. Harrington, Time Harmonic Electromagnetic Fields (McGraw-Hill, 1961).

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

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

4. J. L. Jaramillo, R. P. Macedo, and L. Al Sheikh, “Pseudospectrum and Black Hole Quasinormal Mode Instability,” Phys. Rev. X 11(3), 031003 (2021). [CrossRef]  

5. Y. Ashida, Z. Gong, and M. Ueda, “Non-Hermitian Physics,” Adv. Phys. 69(3), 249–435 (2020). [CrossRef]  

6. G. Gamow, “Quantum Theory of the Atomic Nucleus,” Eur. Phys. J. A 51(3-4), 204–212 (1928). [CrossRef]  

7. A. J. F. Siegert, “On the Derivation of the Dispersion Formula for Nuclear Reactions,” Phys. Rev. 56(8), 750–752 (1939). [CrossRef]  

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

9. R. M. More, “Theory of decaying states,” Phys. Rev. A 4(5), 1782–1790 (1971). [CrossRef]  

10. R. M. More and E. Gerjuoy, “Properties of Resonance Wave Functions,” Phys. Rev. A 7(4), 1288–1303 (1973). [CrossRef]  

11. B. J. Hoenders, “On the completeness of the natural modes for quantum mechanical potential scattering,” J. Math. Phys. 20(2), 329–335 (1979). [CrossRef]  

12. W. H. Press, “Long wave trains of gravitational waves from a vibrating black hole,” Astrophys. J. 170, L105–L108 (1971). [CrossRef]  

13. K. D. Kokkotas and B. G. Schmidt, “Quasi-normal modes of stars and black holes,” Living Rev. Relativ. 2(1), 2 (1999). [CrossRef]  

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

15. 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(9), 5187–5198 (1990). [CrossRef]  

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

17. A. W. Snyder and J. D. Love, Optical Waveguide Theory (Chapman and Hall, 1983).

18. S. Hein, T. Hohage, and W. Koch, “On resonances in open systems,” J. Fluid Mech. 506, 255–284 (2004). [CrossRef]  

19. 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(5), 805–817 (1996). [CrossRef]  

20. 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(9), 1409–1417 (1999). [CrossRef]  

21. P. T. Kristensen, K. Herrmann, F. Intravaia, and K. Busch, “Modeling electromagnetic resonators using quasinormal modes,” Adv. Opt. Photonics 12(3), 612–708 (2020). [CrossRef]  

22. T. Wu, M. Gurioli, and P. Lalanne, “Nanoscale Light Confinement: the Q’s and V’s,” ACS Photonics 8(6), 1522–1538 (2021). [CrossRef]  

23. S. Both and T. Weiss, “Resonant states and their role in nanophotonics,” Semicond. Sci. Technol. 37(1), 013002 (2022). [CrossRef]  

24. N. Moiseyev, “Quantum theory of resonances: calculating energies, widths and cross-sections by complex scaling,” Phys. Rep. 302(5-6), 212–293 (1998). [CrossRef]  

25. N. Moiseyev, Non-Hermitian quantum mechanics, (Cambridge University, 2011).

26. C. E. Baum, Interaction Note 88, December 1971. See https://apps.dtic.mil/sti/citations/ADA066905

27. C. E. Baum, “The singularity expansion method,” p. 129-179 in Transient Electromagnetic Fields Springer Berlin Heidelberg, L. B. Felsen ed. (1976).

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

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

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

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

32. 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(23), 237401 (2013). [CrossRef]  

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

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

35. M. D. Truong, A. Nicolet, G. Demésy, and F. Zolla, “Continuous family of exact Dispersive Quasi-Normal Modal (DQNM) expansions for dispersive photonic structures,” Opt. Express 28(20), 29016–29032 (2020). [CrossRef]  

36. P. Y. Chen, D. J. Bergman, and Y. Sivan, “Generalizing normal mode expansion of electromagnetic Green’s tensor to open systems,” Phys. Rev. Appl. 11(4), 044018 (2019). [CrossRef]  

37. D. A. B. Miller, “Waves, modes, communications, and optics: a tutorial,” Adv. Opt. Photonics 11(3), 679–825 (2019). [CrossRef]  

38. MAN (Modal Analysis of Nanoresonators). QNM solvers and toolboxes are freely available at https://www.lp2n.institutoptique.fr/light-complex-nanostructures

39. 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,” Comput. Phys. Commun. 257, 107509 (2020). [CrossRef]  

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

41. C. Sauvan, J.-P. Hugonin, and P. Lalanne, “Photonic and plasmonic nanoresonators: a modal approach,” Proc. SPIE 9546, 95461C (2015). [CrossRef]  

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

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

44. 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(1), 017802 (2017). [CrossRef]  

45. T. Christopoulos, O. Tsilipakos, G. Sinatkas, and E. E. Kriezis, “On the calculation of the quality factor in contemporary photonic resonant structures,” Opt. Express 27(10), 14505–22 (2019). [CrossRef]  

46. Y. Brûlé, B. Gralak, and G. Demésy, “Calculation and analysis of the complex band structure of dispersive and dissipative two-dimensional photonic crystals,” J. Opt. Soc. Am. B 33(4), 691–702 (2016). [CrossRef]  

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

48. C. Sauvan, “Quasinormal modes expansions for nanoresonators made of absorbing dielectric materials: study of the role of static modes,” Opt. Express 29(6), 8268–8682 (2021). [CrossRef]  

49. M. Parto, Y. G. N. Liu, B. Bahari, M. Khajavikhan, and D. N. Christodoulides, “Non-Hermitian and topological photonics: optics at an exceptional point,” Nanophotonics 10(1), 403–423 (2020). [CrossRef]  

50. G. B. Arfken and H. J. Weber, Mathematical methods for physicists, 6th edition, (Elsevier Academic, 2005, Chap. 7,).

51. Exceptions occurs when one tries to lower substrate reflection by choosing the superstrate medium to have the same refractive index as the substrate.

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

53. In 2D, even for the ideal case of resonators in uniform backgrounds, there is always a singularity in 2D owing to the singularity of the 2D Hankel function at the origin r=0.

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

55. 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(4), 043806 (2018). [CrossRef]  

56. F. Betz, F. Binkowski, and S. Burger, “RPExpand: Software for Riesz projection expansion of resonance phenomena,” SoftwareX 15, 100763 (2021). [CrossRef]  

57. S. Burger, personal communication (2021).

58. 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(8), 085433 (2018). [CrossRef]  

59. T. Wu, D. Arrivault, M. Duruflé, A. Gras, F. Binkowski, S. Burger, W. Yan, and P. Lalanne, “Efficient hybrid method for the modal analysis of optical microcavities and nanoresonators,” J. Opt. Soc. Am. A 38(8), 1224–1231 (2021). [CrossRef]  

60. H. Zhang and O. D. Miller, “Quasinormal Coupled Mode Theory,” arXiv:2010.08650 (2020).

61. K. G. Cognée, Hybridization of open photonic resonators PhD dissertation, University of Bordeaux and University of Amsterdam (2020).

62. C. Tao, J. Zhu, Y. Zhong, and H. T. Liu, “Coupling theory of quasinormal modes for lossy and dispersive plasmonic nanoresonators,” Phys. Rev. B 102(4), 045430 (2020). [CrossRef]  

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

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

65. A. Gras, P. Lalanne, and M. Duruflé, “Non-uniqueness of the Quasinormal Mode Expansion of Electromagnetic Lorentz Dispersive Materials,” J. Opt. Soc. Am. A 37(7), 1219–1228 (2020). [CrossRef]  

66. The term ‘PML modes’ or ‘Bérenger modes’ is particularly used to describe the modes that are obtained with finite-thickness PMLs and are intended to replace the continuous spectrum of the operator resulting from the coupling with the open space. However, the PML modes are also obtained for discretized operator in practice. Discretization always induces errors especially at high frequency owing to the necessarily finite mesh elements. In addition, with non-dispersive PMLs that are quite often used, numerical errors might also occur at low frequencies. Therefore, the ideal PML modes are in practice contaminated by the unperfect discretization. To account for that point, we prefer to call the ‘PML modes’ ‘numerical modes’ in the review.

67. W. Yan, P. Lalanne, and M. Qiu, “Shape Deformation of Nanoresonator: A Quasinormal-Mode Perturbation Theory,” Phys. Rev. Lett. 125(1), 013901 (2020). [CrossRef]  

68. 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(6), 063811 (2019). [CrossRef]  

69. E. A. Muljarov, “Full electromagnetic Green's dyadic of spherically symmetric open optical systems and elimination of static modes from the resonant-state expansion,” Phys. Rev. A 101(5), 053854 (2020). [CrossRef]  

70. The unicity of the α‘s singularly contrasts with the case of expansions based on pole expansions, for which there is an infinity of mathematically-correct formulas owing to the fact that the QNM basis is overcomplete [71], see the Eqs. (2.23) and (2.24) in [16] and related discussions in [29,47,52].

71. L. W. Pearson, “Evidence that bears on the left half plane asymptotic behavior of the SEM expansion of surface currents,” Electromagnetics 1(4), 395–402 (1981). [CrossRef]  

72. A. I. Baz, Ya. B. Zel’dovich, and A. M. Perelomov, Scattering, Reactions, and Decay in Non-Relativistic Quantum Mechanics (IPST, 1969) (Moskva, Nauka, 1971).

73. B. Stout, R. Colom, N. Bonod, and R. C. McPhedran, “Spectral expansions of open and dispersive optical systems: Gaussian regularization and convergence,” New J. Phys. 23(8), 083004 (2021). [CrossRef]  

74. 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(11), 113048 (2014). [CrossRef]  

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

76. 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(21), 213901 (2019). [CrossRef]  

77. S. Franke, M. Richter, J. Ren, A. Knorr, and S. Hughes, “Quantized quasinormal mode description of non-linear cavity QED effects from coupled resonators with a Fano-like resonance,” Phys. Rev. Res. 2(3), 033456 (2020). [CrossRef]  

78. J.-P. Hugonin and P. Lalanne, “Perfectly-matched-layers as nonlinear coordinate transforms: a generalized formalization,” J. Opt. Soc. Am. A 22(9), 1844–1849 (2005). [CrossRef]  

79. W. C. Chew and W. H. Weedon, “A 3D perfectly matched medium from modified Maxwell’s equations with stretched coordinates,” Microw. Opt. Technol. Lett. 7(13), 599–604 (1994). [CrossRef]  

80. J. P. Berenger, Perfectly Matched Layer (PML) for Computational Electromagnetic (Morgan & Claypool, 2007).

81. Our presentation remains simplistic, and the PML concept is not fully general. It may completely fail for instance for geometries that are not spatially invariant, e.g., half spaces with periodic waveguides. Even translation invariant half spaces supporting backward waves with opposite phase and group velocities may lead to insurmountable difficulties, see for instance E. Bécache, S. Fauqueux, and P. Joly, “Stability of perfectly matched layers, group velocities and anisotropic waves,” J. Comput. Phys. 188(2), 399–433 (2003). [CrossRef]  

82. Note that x' if often denoted simply x, which leads to a confusion since this x is different from the original x coordinate in Eq. (1). In practice, in numerics, x' corresponds to the numerical space inside the PML.

83. A. Nicolet, J.-F. Remacle, B. Meys, A. Genon, and W. Legros, “Transformation methods in computational electromagnetics,” J. Appl. Phys. 75(10), 6036–6038 (1994). [CrossRef]  

84. U. Leonhardt and T. G. Philbin, “General relativity in electrical engineering,” New J. Phys. 8(10), 247 (2006). [CrossRef]  

85. P. Lalanne, W. Yan, A. Gras, C. Sauvan, J.-P. Hugonin, M. Besbes, G. Demesy, 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(4), 686–704 (2019). [CrossRef]  

86. The literature has taken the habit to distinguish two kinds of modes, the QNMs and the PML modes. To lift it and also to be more general and encompass methods that do not rely on PMLs, e.g., integral methods [see D. A. Powell, “Interference between the Modes of an All-Dielectric Meta-atom,” Phys. Rev. Appl. 7(3), 034006 (2017)], hereafter we will call the PML modes ‘numerical modes’. The latter are affected by the boundary condition (e.g., hard wall for finite PMLs), whereas the former are not. [CrossRef]  

87. D. Marcuse, Theory of Dielectric Optical Waveguides2nd ed. (Academic, 1991).

88. A. Gras, W. Yan, and P. Lalanne, “Quasinormal-mode analysis of grating spectra at fixed incidence angles,” Opt. Lett. 44(14), 3494–3497 (2019). [CrossRef]  

89. A. Raman and S. Fan, “Photonic Band Structure of Dispersive Metamaterials Formulated as a Hermitian Eigenvalue Problem,” Phys. Rev. Lett. 104(8), 087401 (2010). [CrossRef]  

90. The demonstration in [34] relies on recognizing that the permittivity of any linear dispersive material can be modeled with a N-pole Lorentz-pole permittivity, with increasing accuracy as the number of poles N increases.

91. The term biorthogonal (rather than orthogonal) product is used because the product involves right and left eigenvectors, which are identical for reciprocal systems, see Section3.4.3 and Appendix C, where the normalization of QNMs or resonators with non-reciprocal materials is considered.

92. F. Olyslager, “Discretization of Continuous Spectra Based on Perfectly Matched Layers,” SIAM J. Appl. Math. 64(4), 1408–1433 (2004). [CrossRef]  

93. By most matrices, we intend to mean that the set of diagonalizable matrices (in fact, even the set of matrices with distinct eigenvalues), is a set that is dense in the algebra of matrices with complex coefficients.

94. R. El-Ganainy, K. G. Makris, M. Khajavikhan, Z. H. Musslimani, S. Rotter, and D. N. Christodoulides, “Non-Hermitian physics and PT symmetry,” Nat. Phys. 14(1), 11–19 (2018). [CrossRef]  

95. F. Zolla, A. Nicolet, and G. Demésy, “Photonics in highly dispersive media: the exact modal expansion,” Opt. Lett. 43(23), 5813–5816 (2018). [CrossRef]  

96. The existence of an infinity of expressions for the excitation coefficient should not be confused with the so called ‘overcompleteness’ of QNM expansions inside resonators places in free space. Overcompleteness is due to the fact that only a compact subspace (the resonator inside) of the whole open space is used for the expansion. With regularization, the expansion holds over the whole space and the infinity of expressions arises from the infinity of auxiliary fields for dispersive materials only.

97. In fact, completeness is likely to be achieved in the whole computational domain, even in the PML regions.

98. 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(22), 6359–6362 (2014). [CrossRef]  

99. R. Faggiani, J. Yang, R. Hostein, and P. Lalanne, “Implementing structural slow light on short length scales: the photonic speed-bump,” Optica 4(4), 393–399 (2017). The normalization of cavities coupled to periodic waveguides has also been studied in [98] by regularizing the divergent field in the waveguide with advanced theorems of complex analysis. It would be interesting to compare both approaches. [CrossRef]  

100. N. Rivera and I. Kaminer, “Light-matter interactions with photonic quasiparticles,” Nat. Rev. Phys. 2(10), 538–561 (2020). [CrossRef]  

101. J. Zimmerling, L. Wei, P. Urbach, and R. Remis, “A Lanczos model-order reduction technique to efficiently simulate electromagnetic wave propagation in dispersive media,” J. Comput. Phys. 315, 348–362 (2016). [CrossRef]  

102. T. J. Kippenberg and K. J. Vahala, “Cavity Optomechanics: Back-Action at the Mesoscale,” Science 321(5893), 1172–1176 (2008). [CrossRef]  

103. 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(4), 043822 (2006). [CrossRef]  

104. A. M. Shaltout, V. M. Shalaev, and M. L. Brongersma, “Spatiotemporal light control with active metasurfaces,” Science 364(6441), eaat3100 (2019). [CrossRef]  

105. Q. Zhou, P. Zhang, and X.-W. Chen, “General Framework of Canonical Quasinormal Mode Analysis for Extreme Nano-optics,” Phys. Rev. Lett. 127(26), 267401 (2021). [CrossRef]  

106. P. J. Feibelman, “Surface electromagnetic fields,” Prog. Surf. Sci. 12(4), 287–407 (1982). [CrossRef]  

107. Q. Zhou, P. Zhang, and X.-W. Chen, “Quasinormal mode theory for nanoscale electromagnetism with quantum surface responses,” arXiv:2105.06328 (2021).

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

109. 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(4), 876–891 (2009). [CrossRef]  

110. P. Lalanne, “Mode volume of electromagnetic resonators: let us try giving credit where it is due,” arXiv:2011.00218 (2020).

111. 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(1), 011803 (2013). [CrossRef]  

112. M. Benzaouia, J. D. Joannopoulos, S. G. Johnson, and A. Karalis, “Quasi-normal mode theory of the scattering matrix, enforcing fundamental constraints for truncated expansions,” Phys. Rev. Res. 3(3), 033228 (2021). [CrossRef]  

113. 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(5), 3982–3989 (1994). [CrossRef]  

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

115. P.T. Kristensen and S. Hughes, “Modes and Mode Volumes of Leaky Optical Cavities and Plasmonic Nanoresonators,” ACS Photonics 1(1), 2–10 (2014). [CrossRef]  

116. S.-Y. Lee, “Decaying and growing eigenmodes in open quantum systems: Biorthogonality and the Petermann factor,” Phys. Rev. A 80(4), 042104 (2009). [CrossRef]  

117. 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(22), 27371–27382 (2013). [CrossRef]  

118. QNMPole is one of the two QNM solvers of the freeware MAN (Modal Analysis of Nanoresonators) freely available at https://www.lp2n.institutoptique.fr/light-complex-nanostructures.

119. G. M. Akselrod, C. Argyropoulos, T. B. Hoang, C. Ciracì, C. Fang, J. Huang, D. R. Smith, and M. H. Mikkelsen, “Probing the mechanisms of large Purcell enhancement in plasmonic nanoantennas,” Nat. Photonics 8(11), 835–840 (2014). [CrossRef]  

120. K. Santosh, O. Bitton, L. Chuntonov, and G. Haran, “Vacuum Rabi splitting in a plasmonic cavity at the single quantum emitter limit,” Nat. Commun. 7(1), ncomms11823 (2016). [CrossRef]  

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

122. S. Fahr, C. Rockstuhl, and F. Lederer, “Metallic nanoparticles as intermediate reflectors in tandem solar cells,” Appl. Phys. Lett. 95(12), 121105 (2009). [CrossRef]  

123. T. Tanabe, M. Notomi, E. Kuramochi, A. Shinya, and H. Taniyama, “Trapping and delaying photons for one nanosecond in an ultrasmall high-Q photonic-crystal nanocavity,” Nat. Photonics 1(1), 49–52 (2007). [CrossRef]  

124. E. Silberstein, P. Lalanne, J.P. Hugonin, and Q. Cao, “On the use of grating theory in integrated optics,” J. Opt. Soc. Am. A 18(11), 2865–2875 (2001). [CrossRef]  

125. A. Archambault, F. Marquier, J.-J. Greffet, and C. Arnold, “Quantum theory of spontaneous and stimulated emission of surface plasmons,” Phys. Rev. B 82(3), 035411 (2010). [CrossRef]  

126. T. Christopoulos, O. Tsilipakos, and E. Kriezis, “Perturbation theory for Kerr nonlinear leaky cavities,” Opt. Lett. 45(23), 6442–6445 (2020). [CrossRef]  

127. V. S. Asadchy, M. S. Mirmoosa, A. Díaz-Rubio, S. Fan, and S. A. Tretyakov, “Tutorial on electromagnetic nonreciprocity and its origins,” Proc. IEEE 108(10), 1684–1727 (2020). [CrossRef]  

128. T. Wu and P. Lalanne, “QNMnonreciprocal_resonators: an openly available toolbox for computing the Quasinormal Modes of nonreciprocal resonators,” arXiv:2106.05502 (2021).

129. W. H. Press, S. A. Teukolsky, W. T. Vetterling and, and B. P. Flannery, Numerical Recipes (Cambridge University, 2007).