Modal Analysis of photonic and plasmonic resonators

Quasi-normal modes (QNMs) are ubiquitous throughout photonics and are utilized in a wide variety of applications, but determining these modes remains a formidable task in general. Here we show that by exploiting the structure of Maxwell's equations it is possible to effectively compute QNMs of photonic and plasmonic nanoresonators. The symmetry of Maxwell's equations allows for a reduction to a system of small order via a Lanczos reduction process through which dominant QNMs can be identified. A closed-form reduced-order model for the spontaneous decay (SD) rate of a quantum emitter is also obtained, which does not require an a priori QNM expansion of the fields. The model is parametric in wavelength and field expansions in dominant QNMs are determined a posteriori. We demonstrate and validate that QNMs of open resonators and the SD rate of a quantum emitter are accurately predicted.

Quasi-normal modes (QNMs) are ubiquitous throughout photonics and are utilized in a wide variety of applications, but determining these modes remains a formidable task in general. Here we show that by exploiting the structure of Maxwell's equations it is possible to effectively compute QNMs of photonic and plasmonic nanoresonators. The symmetry of Maxwell's equations allows for a reduction to a system of small order via a Lanczos reduction process through which dominant QNMs can be identified. A closed-form reduced-order model for the spontaneous decay (SD) rate of a quantum emitter is also obtained, which does not require an a priori QNM expansion of the fields. The model is parametric in wavelength and field expansions in dominant QNMs are determined a posteriori. We demonstrate and validate that QNMs of open resonators and the SD rate of a quantum emitter are accurately predicted.
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][2][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 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 γ denoting the decay rate of the emitter in the resonator configuration of interest and γ 0 the decay rate of the same emitter in a reference medium, γ/γ 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. Explicitly, for an emitter located at x = x S and an electric dipole of the form J ext = ∂ t p(t)δ(x−x S ) with dipole moment p(t) = p(t)n s , p(t) = |p(t)|, and n s a unit vector, we have in steady-stateĴ ext = −iωp(ω)δ(x − x S )n s and the time-averaged radiated power is given by To evaluate this power over a frequency or wavelength interval of interest, the electric field strength at the dipole location is required for all frequencies belonging to this interval.
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][7][8][9][10][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)).
In this letter we show that by exploiting the symmetry of the first-order Maxwell system, it is possible to efficiently determine QNMs of open resonating structures consisting of dispersive metallic nanoparticles. In addition, we show that the SD rate can be computed without any a priori mode selection, that is, the decay rate can be computed without an explicit mode expansion of the fields as is more commonly done in decay rate computations (see, e.g. [7]).
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 asD = εÊ +P = ε c (ω)Ê with ε = ε 0 ε ∞ , where ε ∞ is the instantaneous (high-frequency) permittivity and a polarization vectorP that is related to the electric field strength via the generic constitutive relation determine what type of relaxation is considered (Drude, Lorentz). For a Drude model, for example, we have β 0 = ε 0 ω 2 p , β 1 = 0, and β 2 = γ p , where ω p is the volume plasma frequency and γ p the collision frequency of the metal.
Introducing the auxiliary field variableÛ = iωP, we can write the above constitutive relation and Maxwell's equations in the consistent first-order form [12] which can be written as (D + S − iωM)F = −Q, where S and M are medium matrices, and the curl operators of Maxwell's equations are contained in the spatial differentiation operator D. The electromagnetic field quantities and external sources are collected in the field vectorF and source vectorQ, respectively. For most external sources used in practice (electric dipole, for example), the frequency dependence of the source can be factored out and we writeQ =p(ω)Q , wherep(ω) is the source wavelet and Q is frequency independent. We note that the partial differential operator in Eq. (2) can be symmetrized by scaling the second row with β 1 β −1 0 , the third row with −β −1 0 , and the fourth row by −1. The efficiency of our method is based upon this symmetry. Furthermore, measured (causal) material behavior can 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.
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 [12,13]. 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) [14] in which the spatial coordinates are stretched using frequency dependent stretching functions [15]. 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 [16,17], 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 where D contains the discretized curl operators, S and M are the discretized medium matrices, andf cs and q are the discretized field and source vector, respectively. The above system is not conjugate-symmetric with respect to frequency and its time-domain counterpart is unstable due to the application of a frequency-independent PML. However, conjugate-symmetric frequency-domain field approximations can be obtained from the above system as [16]f where is the complex Heaviside unit step function defined as χ(z) = 1 for Re(z) > 0 and χ(z) = 0 for Re(z) < 0. Note thatf(ω) is conjugate-symmetric, that is, it satisfieŝ f * (ω) =f(−ω), provided thatp is conjugate-symmetric. For practical three-dimensional problems direct evaluation of Eq. (4) is usually not feasible, since the order n of the Maxwell system matrix A is simply too large (in 3D, typically n = O(10 6−7 )). It can be shown, however, that matrix A satisfies a particular symmetry property that allows for efficient Lanczos model-order reduction. In particular, defining the bilinear form x|y WM = x T MWy = x T WMy, for vectors from C n , where W is a specific step size matrix [12,18], it can be shown that Ax|y WM = x|Ay WM for all vectors x, y ∈ C n . Moreover, the bilinear form f|f WM is a discrete approximation of the integral which in the literature is used to normalize QNMs [7].
In 1931, Krylov [19] used what are now called polynomial Krylov subspaces in his analysis of oscillations of mechanical systems (e.g. ships). Here, the symmetry of matrix A allows us to follow a similar approach. Specifically, the symmetry of A can be used to reduce this matrix to tridiagonal form using a three-term Lanczos-type recurrence relation [12,18]. Carrying out m steps of this reduction process, we obtain the decomposition where T m is a tridiagonal matrix of order m n containing the Lanczos recurrence coefficients and V m is a tall n-by-m matrix with a column partitioning V m = (v 1 , v 2 , ..., v m ). The columns of matrix V m are referred to as Lanczos vectors, which are taken to be quasiorthonormal i.e., v i |v j WM = δ ij , where δ ij is the Kronecker delta. Furthermore, β m+1 in Eq. (6) is a Lanczos recurrence coefficient and e m is the mth canonical basis vector. To find an approximate spectrum of the Maxwell system matrix A, the Lanczos reduction process can be started with any (randomly generated) starting vector v 1 satisfying v 1 |v 1 WM = 1. If, however, modes excited by a given external source are of interest (as in SD rate computations, for example) then we take v 1 = q|q −1/2 WM q as a starting vector in the reduction process.
The Lanczos decomposition of Eq. (6) serves as a starting point for our modal analysis and SD rate computations. First, as is well known [20], the decomposition can be used to find approximate QNMs of the open resonator system. Specifically, if (θ Second, for a given external source q the decomposition can be used to construct the reduced-order model (ROM) [12,16] f m (ω) = iωp(ω) q|q which gives an approximation of the three-dimensional field of order m. In SD rate computations, however, only the projection of the electric field onto the direction of the dipole moment at the dipole location is of interest. For this projection, we haveÊ(x S , ω) · n s ≈ f m (ω)|q WM and substitution in Eq. (1) gives the ROM for the radiated power P m (ω) = P a Re e 1 |Ĝ(T m , ω)|e 1 (8) with P a = 0.5 ω 2 |p| 2 q|q WM . Only filtered resolvents of the reduced tridiagonal matrix T m need to be computed to evaluate this power over a complete frequency (wavelength) interval of interest and no a priori expansion of the fields in QNMs in required. Explicitly, assuming that T m can be diagonalized and arranging its eigenvectors as columns in matrix Z m = (z  where w k is the kth element of |Z T m e 1 . As mentioned above, converged QNMs can be identified by computing the residual of the approximate QNMs and their contribution to the radiated power P m (ω) can be determined using the spectral expansion of Eq. (9). With I QNM denoting the index set of converged QNMs that contribute to the radiated power we then arrive at a low order QNM expansion by replacing the sum in Eq. (9) by a sum over all k ∈ I QNM . In other words, the Lanczos decomposition allows us to determine a posteriori which converged QNMs actually contribute to the radiated power and ultimately the SD rate of the quantum emitter.
To validate the presented approach, we compute the Purcell factor of a golden nanorod that has been considered in the literature before [7]. The configuration consists of a vertically oriented dipole centered 10 nm above a 30 nm×100 nm golden nanorod, embedded in a dispersionless background material with relative permittivity ε r = 2.25. A Drude model is used as a dispersion model for gold with a plasma frequency ω p = 1.26 · 10 16 Hz and a collision frequency γ p = 1.41 · 10 14 Hz. This dispersion model is used throughout this letter. 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 FEM-RCWA method [7], while the dashed line shows the converged reduced-order model response obtained via Lanczos reduction. 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 mil- lion, while the order of the converged reduced system is m = 4500. Dominant QNMs can be identified from the spectrum of the reduction matrix T 4500 . For this configuration, it turns out that essentially only a single QNM with a complex resonance wavelength of λ = 926+47i 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 Re(Ê z ) and Re(Ê x ) of the dominant QNM as computed via Lanczos reduction are shown in Figs. 1(b) and 1(c), respectively, where the red and blue surfaces have opposite signs. The isosurface has been chosen to best visualize the field distribution.
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 nm × 40 nm × 20 nm nanoplate as illustrated in Fig. 2(c). The wavelength of interest now runs from 0.5 µm to 1.2 µm so that the contribution of higher order QNMs can be investigated. The Purcell factor is computed using the ROM of Eq. (8) and the converged model is shown in Fig. 2 (solid line). 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 (θ 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 along with the contribution of each QNM separately. The real parts of theÊ 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.
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. 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 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 λ = 1034 + 34i nm and λ = 891 + 68i nm, respectively. Figure 3(a) shows an isosurface plot of Re(Ê x ) of the fundamental symmetric QNM, whereas isosurface plots of Re(Ê x/y/z ) of the anti-symmetric resonance are shown Figs. 3(c) -(e). Finally, a higher harmonic antisymmetric resonance is depicted in Fig. 3(b).
In conclusion, we have shown that 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 firstorder 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 Lanczostype reduction process and to find QNMs that are quasiorthonormal with respect to the bilinear form. Moreover, we have presented a new closed-form reduced-order model for the SD rate of a quantum emitter that 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 super- imposing the most contributing modes until a specified error criterion is met. 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.