## Abstract

Determining the electromagnetic field response of photonic and plasmonic resonators is a formidable task in general. Field expansions in terms of quasi-normal modes (QNMs) are often used, since only a few of these modes are typically required for an accurate field description. We show that by exploiting the structure of Maxwell’s equations, conjugate-symmetric frequency-domain field expansions can be efficiently computed via a Lanczos-type algorithm. Dominant QNMs can be identified *a posteriori* with error control and without *a priori* mode selection. Discrete QNM approximations of resonating nanostructures are presented and the spontaneous decay rate of a quantum emitter is also considered.

© 2020 Optical Society of America under the terms of the OSA Open Access Publishing Agreement

## 1. Introduction

Optical nanoresonators enable us to confine electromagnetic energy to subwavelength domains and give rise to locally enhanced fields that may stimulate various optical processes in a wide variety of applications and research areas such as biophotonics, optical antennas, and diffraction gratings [1–3]. Resonators consisting of metallic nanoparticles that are excited by femtosecond laser pulses are often of particular interest [4], since such resonators allow for the control of light-matter interactions with nanometer and subfemtosecond precision in space and time, respectively, thereby enabling new and exciting applications in cell biology and quantum optics, for example. Moreover, metallic nanoparticles are also often used in resonating structures designed to enhance the spontaneous decay (SD) rate of a quantum emitter that is embedded in such a structure, since this rate depends on the surroundings of the emitter and can be enhanced by an electromagnetic resonance (Purcell effect). The spontaneous decay of a quantum emitter is a purely quantum mechanical effect, but can be computed classically in the so-called weak-coupling regime [5]. Specifically, with $\gamma$ denoting the decay rate of the emitter in the resonator configuration of interest and $\gamma _0$ the decay rate of the same emitter in a reference medium, we have $\gamma /\gamma _0 = P/P_0$, where $P$ and $P_0$ are the time-averaged powers radiated by an electric dipole positioned at the location of the emitter in the resonator configuration and reference medium, respectively. For an emitter located at $\textbf {x}= \textbf {x}_{\textrm {S}}$ and an electric dipole of the form $\textbf {J}^{\textrm {ext}}=\partial _t \textbf {p}(t) \delta ( \textbf {x}- \textbf {x}_{\textrm {S}})$ with dipole moment $\textbf {p}(t)=p(t) \textbf {n}_{\textrm {s}}$, and $\textbf {n}_{\textrm {s}}$ a unit vector, we have in steady-state (time factor $\exp (-\textrm {i}\omega t$)) $\hat { \textbf {J}}^{\textrm {ext}}=-\textrm {i}\omega \hat {p}(\omega ) \delta ( \textbf {x}- \textbf {x}_{\textrm {S}}) \textbf {n}_{\textrm {s}}$ and the time-averaged radiated power is given by

To investigate what local field or decay rate enhancements can be realized, a modal analysis of a resonating structure is typically carried out. For open resonator structures these modes are called Quasi Normal Modes or QNMs and are characteristic of the structure at hand and independent of the excitation. An external source (or incident field) determines what resonant modes are actually excited, while the contribution of these excited modes to a measured field response is determined by the receiver. In open resonant structures, typically only a small number of QNMs are necessary to accurately model measured field responses [6–11] and in SD rate computations the source and receiver location actually coincide, since the electric field strength at the source (dipole) location is required to determine the radiated power (see Eq. (1)).

Computing QNMs for three-dimensional dispersive structures is not trivial and computationally demanding. In this paper, we first briefly review and expand upon our previous work and show that by exploiting the symmetry of the first-order Maxwell system, it is possible to efficiently determine QNMs via a Lanczos-type reduction algorithm. Furthermore, we show here how the symmetry of the Maxwell system is intimately related to QNM normalization. Subsequently, we show how the algorithm can be used to formulate frequency-domain modal field expansions that are complex-symmetric with respect to frequency, which corresponds to real-valued and causal field expansions in the time-domain. Furthermore, in this paper we demonstrate how this field expansion can be used to determine dominant QNMs *a posteriori* and with error control as opposed to the more common approach, which uses *a priori* mode selection (see, e.g. [7]). In other words, by using the proposed modal expansion technique it is possible to identify which modes are dominant and actually contribute to a measured field response after the Lanczos reduction process is finished.

Finally, we illustrate the effectiveness of the proposed modal expansion approach by determining dominant QNMs of several three-dimensional metallic nanoresonators consisting of dispersive materials. The SD rate of a quantum emitter in a resonating structure is also considered. In this latter case, the Lanczos reduction process even provides us with a closed-form decay rate formula, since this quantity is determined by the projection of the electric field onto the direction of the dipole source at the source location (source and receiver coincide). Results for the SD rate are validated against other methods. Full Lanczos field expansions and dominant QNM expansions are compared as well.

## 2. Basic equations and symmetry

To describe the reaction of a metallic nanoparticle to the presence of an electromagnetic field, we write the electric displacement vector in Maxwell’s equations as $\hat { \textbf {D}}=\varepsilon \hat { \textbf {E}} + \hat { \textbf {P}}=\varepsilon _{\textrm {c}}(\omega ) \hat { \textbf {E}}$ with $\varepsilon =\varepsilon _0 \varepsilon _{\infty }$, where $\varepsilon _{\infty }$ is the instantaneous (high-frequency) permittivity and a polarization vector $\hat { \textbf {P}}$ that is related to the electric field strength via the generic constitutive relation

Introducing the auxiliary field variable $\hat { \textbf {U}} = \textrm {i}\omega \hat { \textbf {P}}$, we can write the above constitutive relation and Maxwell’s equations in the consistent first-order form [13]

We note that the partial differential operator in Eq. (3) can be symmetrized by scaling the second row with $\beta _1\beta _0^{-1}$, the third row with $-\beta _0^{-1}$, and the fourth row by $-1$. The efficiency of our method is based upon this symmetry, as will be explained later. Furthermore, measured (causal) material behavior can also be modeled using this formulation by fitting a rational function representation (i.e. a multipole expansion consisting of a superposition of Lorentz and Drude models) for the complex permittivity to permittivity measurements. This leads to the introduction of multiple auxiliary field variables and the resulting system can be symmetrized in a similar manner as described above [12].

To carry out a modal analysis of arbitrarily-shaped open resonators, we discretize the first-order Maxwell system in space using a staggered finite-difference Yee mesh. We discretize on such a mesh, since it can be shown that the discretization procedure is mimetic, that is, it is structure preserving and conservation laws and important physical symmetry properties of Maxwell’s equations (symmetry related to energy conservation or symmetry related to reciprocity, for example) have a counterpart after discretization [13,14]. Other discretization schemes (finite elements, for example) can also be used, of course, so long as these schemes are mimetic as well.

In addition, radiation towards infinity has to be taken into account, since we are interested in open nanoresonators. Typically, this is realized by surrounding the domain of interest by a so-called Perfectly Matched Layer (PML) [15] in which the spatial coordinates are stretched using frequency *dependent* stretching functions [16]. However, a disadvantage of such an approach is that in two- and three-dimensional problems this leads to nonlinear eigenvalue problems that need to be solved to find dominant QNMs. Therefore, our approach is to apply the PML technique of [17,18], which uses complex spatial step sizes to realize a perfectly matched layer, which do *not* explicitly depend on frequency and leads to linear eigenproblems. Incorporating this PML technique into our spatial discretization scheme then leads to the discretized first-order Maxwell system

For practical three-dimensional problems direct evaluation of Eq. (5) is usually not feasible, since the order $n$ of the Maxwell system matrix $\mathsf {A}$ is simply too large (in 3D, typically $n=O(10^{6}) - O(10^{7}$)). It can be shown, however, that matrix $\mathsf {A}$ satisfies a particular symmetry property that allows for efficient Lanczos model-order reduction. In particular, in the $\mathsf {M} \mathsf {W}$ weighted bilinear form $\mathsf {x}^{T} \mathsf {M} \mathsf {W} \mathsf {y} = \mathsf {x}^{T} \mathsf {W} \mathsf {M} \mathsf {y}$, for vectors from $\mathbb {C}^{n}$, where $\mathsf {W}$ is a specific diagonal step size matrix [13,19], it can be shown that

## 3. Lanczos reduction and causal modal field expansions

In 1931, Krylov [20] used what are now called polynomial Krylov subspaces in his analysis of oscillations of mechanical systems (e.g. ships). Here, the symmetry of matrix $\mathsf {A}$ allows us to follow a similar approach. Specifically, the symmetry of $\mathsf {A}$ can be used to reduce this matrix to tridiagonal form using a three-term Lanczos-type recurrence relation [13,19,21]. Carrying out $m$ steps of this reduction process, we obtain the decomposition

The Lanczos decomposition of Eq. (10) serves as a starting point for our modal analysis and SD rate computations. First, as is well known [22], the decomposition can be used to find approximate QNMs of the open resonator system. Specifically, if the eigenvalue $\theta _j^{[m]}$ and corresponding eigenvector $\mathsf {z}_j^{[m]}$ are an eigenpair of the reduced matrix $\mathsf {T}_m$ then postmultiplication of (10) by $\mathsf {z}_j^{[m]}$ shows that $(\theta _j^{[m]}, \mathsf {V}_m \mathsf {z}_j^{[m]})$ is an approximate eigenpair of $\mathsf {A}$ with a residual vector given by $\beta _{m+1} ( \mathsf {e}_{m}^{T} \mathsf {z}^{[m]}) \mathsf {v}_{m+1}$. Converged QNMs

can be identified by computing the norm of this residual. Note that normalizing the eigenvectors $\mathsf {z}_j^{[m]}$ of $\mathsf {T}_m$ such that $ \left ( { \mathsf {z}_j^{[m]}} \right )^{T} { \mathsf {z}_i^{[m]}} =\delta _{ij}$ ensures that the approximate QNMs $\mathsf {y}_j$ are normalized with respect to the bilinear form (9), i.e. $ { \mathsf {y}_j}^{T} \mathsf {W} \mathsf {M} \, { \mathsf {y}_i}=\delta _{ij}$. Second, for a given external source $\mathsf {q}$ the decomposition can be used to construct the reduced-order model (ROM) [13,17]*a priori*expansion of the fields in QNMs in required. Explicitly, assuming that $\mathsf {T}_m$ can be diagonalized and arranging its eigenvectors as columns in matrix $\mathsf {Z}_m=( \mathsf {z}_1^{[m]}, \mathsf {z}_2^{[m]},\ldots , \mathsf {z}_m^{[m]})$, we have

*a posteriori*which converged QNMs actually contribute to the radiated power and ultimately the SD rate of the quantum emitter. The error of such a resonance expansion can be quantified by comparing it to Eq. (17) and leads to conjugate-symmetric and causal field responses, contrary to other approaches (e.g. [7]).

## 4. Results

#### 4.1 Validation of the reduction method

To validate the presented approach, we compute the Purcell factor of a golden nanorod that has been considered in the literature before [7,19]. The configuration consists of a vertically oriented dipole centered 10 nm above a $30~\textrm {nm}\times 100~\textrm {nm}$ golden nanorod, embedded in a dispersionless background material with relative permittivity $\varepsilon _{\textrm {r}}=2.25$. A Drude model is used as a dispersion model for gold with a plasma frequency $\omega _{\textrm {p}}=1.26 \cdot 10^{16}$ Hz and a collision frequency $\gamma _{\textrm {p}}=1.41\cdot 10^{14}$ Hz. This dispersion model is used throughout this paper. Finally, we mention that since our reduction framework is designed for arbitrarily-shaped nanoresonators, we do not make use of any rotational symmetry.

In Fig. 1 the computed Purcell factor over a complete wavelength interval of interest is shown. The solid line signifies the result obtained with the aperiodic Fourier modal method of [7] that takes the rotational symmetry into account, while the dashed line shows the converged reduced-order model response obtained via Lanczos reduction of the three dimensional model. The computed enhancement factors of both methods are in good agreement with each other. The unreduced Maxwell system has an order of $n=8.6$ million, while the order of the converged reduced system is $m=4500$. Dominant QNMs can be identified from the spectrum of the reduction matrix $\mathsf {T}_{4500}$. For this configuration, it turns out that essentially only a single QNM with a complex resonance wavelength of $\lambda ={{2 \pi c_0/\omega }=} 926 + 47\textrm {i}$ nm contributes to the SD rate over the considered wavelength interval. Higher order QNMs only contribute to the Purcell factor for wavelengths smaller than 600 nm. Finally, isosurface plots of $\textrm {Re}(\hat {E}_z)$ and $\textrm {Re}(\hat {E}_x)$ of the dominant QNM as computed via Lanczos reduction are shown in Figs. 1(b) and 1(c), respectively, where the red and black surfaces have opposite signs. The isosurface has been chosen to best visualize the field distribution.

#### 4.2 A posteriori mode expansion

To demonstrate that the Lanczos reduction technique can also handle configurations in which multiple QNMs contribute to the SD rate, we compute the Purcell factor of a quantum emitter that is placed 10 nm above a $102~\textrm {nm} \times 40~\textrm {nm} \times 20~\textrm {nm}$ nanoplate as illustrated in Fig. 2(c). The wavelength of interest now runs from $0.5~\mu$m to $1.2~\mu$m so that the contribution of higher order QNMs can be investigated. The background permittivity is $\varepsilon _r=2.25$, and the configuration is discretized on a Yee-grid with a stepsize of 2 nm and a PML with 9 primary and dual steps is used, which matches the free space impedance with a maximum error of $10^{-4}$ in the propagative interval of interest and a maximum error of $10^{-3}$ for evanescent waves [18].

The Purcell factor is computed using the ROM of Eq. (16) and the converged model is shown in Fig. 2 (solid line, bottom plot). Without any *a priori* mode selection, a low rank expansion in QNMs can now be obtained by ranking the individual contributions of the approximate eigenpairs $(\theta ^{[m]}_j, \mathsf {V}_m \mathsf {z}^{[m]}_{j})$ to the Purcell factor. For this configuration, we find that essentially only three QNMs are required to accurately describe the Purcell factor on the considered wavelength interval. The resulting three-term QNM expansion is shown in Fig. 2 (black dotted line) along with the contribution of each QNM separately (yellow, red, and green solid lines). This example illustrates that by using Lanczos reduction we are able to identify *a posteriori* which converged QNMs actually contribute to the Purcell factor on a wavelength interval of interest and no *a priori* mode selection for the fields is necessary. The real parts of the $\hat {E}_x$ fields of the contributing QNMs are shown in Figs. 2(a-b) for the higher order modes and in Fig. 2(d) for the fundamental QNM. In the top part of Fig. 2 the relative pointwise error of the QNM expansion is shown, i.e. $(\gamma ^{QNM}-\gamma )/\gamma$.

#### 4.3 Lanczos method as a resonance solver

Finally, to show that QNMs in configurations consisting of multiple dispersive nanoresonators can be determined as well, we place two copies of the golden nanoplate next to each other such that the largest faces are parallel (see Fig. 3). The distance between the plates is 38 nm. This configuration supports anti-symmetric and symmetric resonances, where the wavelength of the anti-symmetric resonance is larger than the wavelength (in vacuum) of the symmetric resonance in accordance with the theory of electronic oscillators. In particular, the wavelengths of the fundamental anti-symmetric and symmetric resonances are $\lambda =1034+34\textrm {i}~$nm and $\lambda =891+68\textrm {i}~$nm, respectively. Figure 3(a) shows an isosurface plot of $\textrm {Re}(\hat {E}_x)$ of the fundamental symmetric QNM, whereas isosurface plots of $\textrm {Re}(\hat {E}_{x/y/z})$ of the anti-symmetric resonance are shown Figs. 3(c) – (e). A higher harmonic anti-symmetric resonance is depicted in Fig. 3(b).

## 5. Conclusion

The symmetry of Maxwell’s equations can be used to effectively compute QNMs of three-dimensional arbitrarily-shaped dispersive nanoresonators. A mimetic discretization of the first-order Maxwell equations for dispersive media leads to a large-scale discretized Maxwell system that is symmetric with respect to a particular bilinear form. This symmetry property allows us to reduce the large-scale Maxwell system to a system of much smaller order via a Lanczos-type reduction process and to find QNMs that are quasi-orthonormal with respect to the bilinear form. Moreover, we have presented resonance expansions for the fields and the SD rate of a quantum emitter that are conjugate-symmetric in frequency leading to causal and real-valued time-domain system responses. The ROM for the SD rate is parametric in wavelength meaning that a single model approximates the SD rate over a complete wavelength interval of interest, i.e. the model allows for wavelength sweeps. This feature is important in many applications in quantum optics, where the SD rate is controlled and optimized by modifying the background configuration of the quantum emitter. Specifically, for each background realization a single ROM provides an SD rate response over a complete wavelength interval of interest, which can significantly speed up the design and optimization of the resonating environment. Furthermore, the ROM does not require an *a priori* expansion of the electric field in terms of QNMs. It is not necessary to determine beforehand which QNMs contribute the most to the electric field at the dipole location. In fact, which modes actually contribute on a given wavelength interval can be determined *a posteriori* from the reduced Lanczos system and the corresponding converged ROM by ranking and superimposing the modes that contribute the most to the SD rate. In this manner, the ROM for the SD rate gives us control over the error that is introduced when a subset of QNMs is used to approximate the SD rate of a quantum emitter.

## Appendix A. Implementation of perfectly matched layers

Open domains are simulated using perfectly matched layers (PMLs). Specifically, coordinate stretching is applied in a layer (the PML) surrounding the domain of interest, which for Cartesian coordinates amounts to replacing the spatial derivatives $\partial _k$ in Maxwell’s equations by $\hat {\chi }^{-1}_k(k,-\textrm {i}\omega )\partial _k$, $k=x,y,z$, where $\hat {\chi }_k(k,-\textrm {i}\omega )$ is a frequency-dependent stretching function that depends of the $k$-coordinate only. Inside the domain of interest we have $\hat {\chi }_k(k,-\textrm {i}\omega )=1$, while within the PML one chooses the stretching functions such that waves entering the PML get attenuated or delayed or both. A typical choice for a stretching function within a PML is

where $a_k$ and $b_k$ are real-valued functions of $k$ only, $k=x,y,z$.Including the stretching functions in Maxwell’s equations leads to a system of the form

One obvious solution is to fix the PML frequency to $\omega =\omega _0$, say, where $\omega _0$ belongs to the frequency interval of interest. We then obtain the linearized system

To address these problems, in [17] a PML implementation is proposed in which optimal imaginary PML step sizes are determined to absorb propagating waves belonging to a frequency interval not necessarily close to a certain fixed frequency $\omega _0$. Specifically, for a PML consisting of $m$ layers, the discrete impedance function $\varphi _m(-\textrm {i}\omega )$ at the interface between the PML and the domain of interest is an $[m-1/m]$ rational function and can be written in pole-residue form as

The residues and poles can be obtained by minimizing the max-norm between this discrete impedance function and the continuous Neumann-to-Dirichlet map on a spectral interval of interest, which amounts to solving a Zolotarev optimal rational approximation problem for which a solution is available (see [17] and the references cited therein). Having found the residues and poles, the PML step sizes can be obtained using the Euclidean polynomial division algorithm. In [18] a PML approach has been proposed that uses complex step sizes with nonvanishing real and imaginary parts to optimally absorb evanescent and propagating waves simultaneously.Finally, having a complex stepsize PML realization available, it is shown in [17] that with a so-called stability-corrected wave function approach real-valued and stable field responses in the time-domain and conjugate-symmetric wave responses in the frequency-domain can be obtained based on these complex stepsize PML realizations.

## Funding

Nederlandse Organisatie voor Wetenschappelijk Onderzoek (14222).

## Disclosures

The authors declare no conflicts of interest.

## References

**1. **F. Vollmer, D. Braun, A. Libchaber, M. Khoshsima, I. Teraoka, and S. Arnold, “Protein detection by optical shift of a resonant microcavity,” Appl. Phys. Lett. **80**(21), 4057–4059 (2002). [CrossRef]

**2. **B. Rolly, B. Stout, and N. Bonod, “Boosting the directivity of optical antennas with magnetic and electric dipolar resonant particles,” Opt. Express **20**(18), 20376–20386 (2012). [CrossRef]

**3. **S. Peng and G. M. Morris, “Efficient implementation of rigorous coupled-wave analysis for surface-relief gratings,” J. Opt. Soc. Am. A **12**(5), 1087–1096 (1995). [CrossRef]

**4. **R. Faggiani, A. Losquin, J. Yang, E. Marsell, A. Mikkelsen, and P. Lalanne, “Modal analysis of the ultrafast dynamics of optical nanoresonators,” ACS Photonics **4**(4), 897–904 (2017). [CrossRef]

**5. **L. Novotny and B. Hecht, * Principles of Nano-Optics*, 2nd Ed. (Cambridge University, 2012).

**6. **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]

**7. **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]

**8. **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]

**9. **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]

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

**11. **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]

**12. **A. Taflove and S. C. Hagness, * Computational Electrodynamics: The Finite-Difference Time-Domain Method*, 2nd Ed. (Artech House, 2000).

**13. **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]

**14. **W. Chew, “Electromagnetic theory on a lattice,” J. Appl. Phys. **75**(10), 4843–4850 (1994). [CrossRef]

**15. **J.-P. Berenger, “A perfectly matched layer for the absorption of electromagnetic waves,” J. Comput. Phys. **114**(2), 185–200 (1994). [CrossRef]

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

**17. **V. Druskin and R. Remis, “A Krylov stability-corrected coordinate-stretching method to simulate wave propagation in unbounded domains,” SIAM J. on Sci. Comput. **35**(2), B376–B400 (2013). [CrossRef]

**18. **V. Druskin, S. Güttel, and L. Knizhnerman, “Near-optimal perfectly matched layers for indefinite Helmholtz problems,” SIAM Rev. **58**(1), 90–116 (2016). [CrossRef]

**19. **J. Zimmerling, L. Wei, P. Urbach, and R. Remis, “Efficient computation of the spontaneous decay rate of arbitrarily shaped 3D nanosized resonators: a Krylov model-order reduction approach,” Appl. Phys. A **122**(3), 158 (2016). [CrossRef]

**20. **A. N. Krylov, “On the numerical solution of the equation by which in technical questions frequencies of small oscillations of material systems are determined,” Otdelenie Matematicheskikh i Estestvennykh Nauk **7**, 491–539 (1931).

**21. **J. Liesen and Z. Strakos, * Krylov Subspace Methods* (Oxford University, 2013).

**22. **G. Golub and C. V. Loan, * Matrix Computations*, 4th Ed. (Johns Hopkins University, 2013).