Nonlocal quasinormal modes for arbitrarily shaped three-dimensional plasmonic resonators

Nonlocal effects have been shown to be responsible for a variety of non-trivial optical effects in small-size plasmonic nanoparticles, beyond classical electrodynamics. However, it is not clear whether optical mode descriptions can be applied to such extreme confinement regimes. Here, we present a powerful and intuitive quasinormal mode description of the nonlocal optical response for three-dimensional plasmonic nanoresonators. The nonlocal hydrodynamical model and a generalized nonlocal optical response model for plasmonic nanoresonators are used to construct an intuitive modal theory and to compare to the local Drude model response theory. Using the example of a gold nanorod, we show how an efficient quasinormal mode picture is able to accurately capture the blueshift of the resonances, the higher damping rates in plasmonic nanoresonators, and the modified spatial profile of the plasmon quasinormal modes, even at the single mode level. We exemplify the use of this theory by calculating the Purcell factors of single quantum emitters, the electron energy-loss spectroscopy spatial maps, as well as the Mollow triplet spectra of field-driven quantum dots with and without nonlocal effects for different size nanoresonators. Our nonlocal quasinormal mode theory offers a reliable and efficient technique to study both classical and quantum optical problems in nanoplasmonics.


I. INTRODUCTION
Fundamental studies of light-matter interactions using plasmonic devices continue to make considerable progress and offering up a wide range of applications [1][2][3][4][5][6][7]. For spatial positions very close to metal resonators, the local Drude model is known to fail, which challenges many of the usual modeling techniques that use the classical Maxwell equations. In particular, charge density oscillations become relevant, causing frequency shifts of the localized surface plasmon (LSP) resonance as well as the appearance of additional resonances above the plasmon frequency [8][9][10][11][12][13]. Such investigations have been performed using both density functional theory (DFT) at the atomistic level [14], and using macroscopic nonlocal Maxwell's equations in the form of the hydrodynamical model (HDM) [13] and a generalized nonlocal optical response (GNOR) model [15]. However, so far, with the exception of the simple cases of spherical or cylindrical nanoparticles, nonlocal investigations have been primarily done using purely numerical simulations [16][17][18], which is not only computationally very extensive for arbitrary shaped plasmonic systems, but also can lack important physical insight; most of these calculations are also restricted to 2D geometries or simple particle shapes. Thus there is now need for more intuitive and efficient formalisms with nonlocal effects included, for arbitrarily shaped metal resonators in a numerically feasible way.
In optics and nanophotonics, one of the most successful analytical approaches to most resonator problems is to adopt a modal picture of the optical cavity (e.g., in cavity-QED and coupled mode theory). Recently, it has been also shown how quasinormal modes (QNMs) can quantitatively describe the dissipative modes of both dielectric cavities and LSP resonances [19], and even hybrid structures of metals and photonic crystals [20]. In contrast to "normal modes," which are solutions to Maxwell's equations subjected to (usually) fixed or periodic boundary conditions, QNMs are obtained with open boundary conditions [21], and they are associated with complex frequencies whose imaginary parts quantify the system losses. These QNMs require a more generalized normalization [21][22][23][24][25][26], allowing for accurate mode quantities to be obtained such as the effective mode volume or Purcell factor [27], i.e., the enhanced spontaneous emission (SE) factor of a dipole emitter. These QNMs are typically computed numerically from the Helmholtz equation with open boundary conditions, e.g., with perfectly matched layers (PMLs), whose solution can then be used to construct the full photon Green function (GF) of the medium-a function that is well known to connect to many useful quantities in classical and quantum optics [28][29][30][31][32][33][34][35][36][37][38]. The GF can also be used (and indeed is required) to compute electron energy loss spectroscopy (EELS) maps for plasmonic nanostructures [39][40][41][42][43][44][45][46][47][48][49][50], which is a notoriously difficult problem in computational electrodynamics, especially for nanoparticles or arbitrary shape. Despite these successes with QNMs, in the presence of nonlocal effects, it is not known if such a mode description even applies.
In this work, we show that, somewhat surprisingly, QNMs can indeed be obtained and used to construct the full system GF for complex 3D plasmonic nanoresonators with nonlocal effects, and even a single mode description is accurate over a wide range of frequencies and spatial positions. We start by redefining the Helmholtz equation that is usually solved to obtain the local QNMs [37], and then extend this approach to the case of nonlocal systems using a generalized Helmholtz equation, which is applicable to both HDM and GNOR models. A semi-analytical modal GF is then used to perform Purcell factor calculations of dipole emitters positioned nearby plasmonic gold nanorods (a structure for which there is no known analytical GF). We then show the accuracy of the modal Purcell factors against fully vectorial dipole calculations, also computed in the presence of the nonlocal corrections. The calculated QNMs are also used to accurately quantify the effective mode volume associated with coupling to quantum emitters, that can be used, e.g., for quantifying single photon source figures of merit [6,51]. Additionally, we examine the size dependence of the nonlocal behavior by investigating nanorods of different sizes, verifying the anticipated LSP blueshifts [11] and damping with decreasing nanoparticle size [15]. Next we use our QNM technique to efficiently calculate the EELS maps for different sizes of nanoparticles [40,47,48]. Finally, to more rigorously show the benefit of our nonlocal modal picture for use in quantum theory of light-matter interaction, we study the behavior of the Mollow triplets of field-driven quantum dots (QDs) coupled to plasmonic resonators [52], under the influence of nonlocal effects.

II. CAVITY MODE APPROACH TO NONLOCAL PLASMONICS
Without nonlocal corrections to the metal, the QNMs, f µ (r), can be defined as the solution to the Helmholtz equation with open boundary conditions (such as PMLs), where ε (r, ω) is the relative dielectric function of the system, andω µ = ω µ − iγ µ is the complex resonance frequency that can also be used to quantify the QNM quality factor, Q µ = ω µ /2γ µ . For metallic regions, the dielectric function can be described using the local Drude model, ε MNP (r, ω) = 1 − ω 2 p /ω (ω + iγ p ) , with ω p = 8.29 eV and γ p = 0.09 eV for the plasmon frequency and collision rate of gold [53], respectively. However, when considering the nonlocal nature of the plasmonic system, the electric field displacement relates to the electric field through an integral equation rather than a simple proportionality [12,54]. In this nonlocal case, a modified set of equations [13,15] can be used to define nonlocal QNMs,f nl µ (r), from where J µ is the induced current density and ξ is a phenomenological length scale associated with the nonlocal corrections [15]. Indeed, where β is the hydrodynamic parameter proportional to the electron Fermi velocity, v F = 1.39 × 10 6 m/s, and D = 2.9 × 10 −4 m 2 /s [55] is the diffusion constant associated with the short-range nonlocal response. While ξ in its full form represents the GNOR model, we can simply switch to the HDM by neglecting the diffusion. Traditionally in cavity physics, the concept of effective mode volume, V eff , plays a key role in characterizing the mode properties; historically, V eff quantifies the degree of light confinement in optical cavities, and is normally defined at the modal antinode where, e.g., a quantum emitter is typically placed. Even though for plasmonic dimers one can reasonably choose the gap center as the place to calculate the mode volume, for plasmonic resonators in general, this simple picture of mode volume is ambiguous. However, one can still quantify an effective inition holds for the local QNM, only one usesf µ ) [19], for rigorous use in Purcell's formula, which is associated with coupling to emitters at different locations outside (but typically near) the metal nanoparticle within a background medium of refractive index n b . Such a positiondependent mode volume can then be used in a generalized Purcell factor, to obtain the SE enhancement rate of a dipole emitter placed at r around a cavity with the resonance wavelength of λ c and quality factor of Q. The quantum emitter is assumed to be on resonance and aligned in polarization with the LSP mode. Recent work has shown that QNMs, when obtained in normalized form (as done in this work), accurately describe lossy plasmonic resonators using the local Drude model [20,37,56]. Here, we extend such an approach to include the nonlocal effects by considering the ansatz, for the scattered GF, which is extremely useful as it can be immediately used to obtain the full position and frequency dependence of the generalized Purcell factor (SE enhancement factor) for a dipole emitter polarized along n: [32] F (r; ω) = 1 + 6πc 3 ω 3 n b n · Im{G nl sc (r, r; ω)} · n.
Note that in a single mode regime, if at the peak of the resonance frequency, ω = 2πc/λ c , then (6) reduces to (4).

III. RESULTS AND EXAMPLE APPLICATIONS
In this section, a range of applications are presented to demonstrate the power and reliability of the QNM theory for optical and quantum optical investigations of plasmonic resonators, down to the nonlocal regime. We emphasize that, once the QNMs are calculated, as discussed in subsection III A, all the example studies are performed in matters of seconds owing to the analytical power of the technique.

A. Local versus nonlocal quasinormal modes
To obtain the system QNMs for the nonlocal HDM/GNOR model defined via (2) and (3), we employ the frequency domain technique discussed in Ref. [25] (used for the local Drude model), where an inverse GF approach is used to return the QNM in normalized form without having to carry out any spatial integral. We extend this method by incorporating nonlocal corrections, both for the HDM method [57] as well as the more complete GNOR model [15,58]. While the system GF using the notation of [25] finds a different form, for consistency, we follow the ansatz of (5) to briefly explain the technique.
The basic idea is to use a dipole excitation at the location of interest, r 0 , having a dipole moment of d, to numerically obtain the scattered Green function (as explained in more detailed later), G sc (r 0 , r 0 ; ω), and then reverse (5) to arrive at where a single QNM,f c (r), with the complex frequencyω c is considered, and we have defined A(ω) = ω 2 /2ω c (ω c − ω) for simplicity. The above quantity is in fact all one needs to perform an integration-free normalization for the QNM, and in practice it is calculated at frequencies very close to the QNM frequency (so is the QNM of (8)) to avoid divergent behavior. When inserted back into (5), one arrives at which provides the full spatial profile of the QNM, given that one also keeps track of the system response at all other locations, G sc (r, r 0 ; ω), within the same simulation. The numerical implementation is done using the frequency domain finite-element solver COMSOL [59], where an electric current dipole source is used to excite the system and iteratively search for the QNM frequencies by monitoring the strength of the system response [25]. To obtain the QNM, one obtains the scattered GF as the difference between two dipole simulation GFs at frequencies very close to the QNM frequency, with and without the metal nanoparticles [25], either in local or nonlocal case. The computed QNM can then provide the full spectral and spatial shape of the resonances involved. In our calculations, a computational domain of 0.5 µm 3 was used for all simulations with the maximum element size of 0.2 nm on the nanoparticle surface and 0.6 nm inside. The maximum element size elsewhere is set to 33 nm to ensure convergent results over a wide range of frequencies, and 10 layers of PML were used. We have checked that these parameters provide accurate numerical convergence for both local and nonlocal simulations done in this work.
Depicted in Fig. 1 are the computed QNMs for three different gold cylindrical nanorods with the same aspect ratio, varying from 20 nm to 5 nm in length (see figure caption for details). The left panels represent the local Drude model QNMs while the right panels show the QNMs using the nonlocal GNOR model. As seen, the main QNM shapes are similar but a redistribution of the localized field clearly occurs due to the inclusion of the nonlocal corrections. While the local Drude model predicts a similar mode shape for the different nanoparticle sizes, the nonlocal corrections introduce a pronounced degree of modal reshaping for smaller nanoparticles. Indeed, even for the largest nanoparticle shown in Figs. 1(a,b), higher field values are seen both inside as well as outside (but near) the metallic region. Figure 2 shows the computed QNM Purcell factors using the local Drude model as well as the two different nonlocal models for the h = 20 nm nanorod. As can be seen, both HDM and GNOR models predict the known blueshift of the plasmonic resonance [12][13][14]. However, the nonlocal prediction of the peak enhancement strongly depends on the model chosen. The GNOR model in particular predicts a considerably lower Purcell factor due to the inclusion of diffusion, which accounts for surfaceenhanced Landau damping [15]. Indeed, as will be discussed shortly, including the nonlocal effects modifies both the quality factor and the mode volume associated with QNMs. The inset also shows the validity of our Pur- for nanorods of different heights, h = 20 nm, h = 10 nm and h = 5 nm, where y0 = 0 is at the center of the nanorods. The same geometrical aspect ratio of 2 is used corresponding to a radius of r = 5 nm for the largest resonator. The double arrow in (a) shows the location of the dipole emitter at 10 nm away from the metallic surface, that is kept the same for all QNM calculations, and the green box represents the metallic border.

B. Purcell factors from coupled dipole emitters
cell factor calculations against full-dipole numerical calculations, only shown for the nonlocal GNOR response. However, a similar degree of agreement is observed for all other calculations both in Fig. 2 and what follows.
In Fig. 2, bottom panel, we additionally plot the corresponding effective mode volume for a range of dipole locations, from the nanorod surface (at z = 10 nm) up to 10 nm away. A comparison between the local Drude model (solid-blue) and the nonlocal GNOR (dashed-red) is shown where a nontrivial difference is observed. Closer to the metallic surface, smaller effective mode volumes are predicted by the nonlocal corrections while further away the opposite takes place. The difference at larger distances however, is mainly due to nonlocal corrected resonant wavelengths that are used to normalize the mode volumes.
We also consider enhanced SE from the three different gold nanoparticles discussed in Fig. 1. Plotted in Fig. 3 are the QNM calculations of the Purcell factors for dipole emitters located 10 nm away from nanorods, on the z-axis. In each case, the local Drude results are compared with the nonlocal GNOR results. Clearly, nonlocal corrections result in both larger resonance blueshifts and larger damping rates (lower quality factors).

C. Computing EELS spatial maps
For our next application of the nonlocal QNM theory, we calculate the spatial maps associated with EELS experiments that are obtained by nanometer-scale resolution in microscopy of LSP resonances [40,41,44,49,50]. Since the GF is available at all locations through QNM expansion of (5), the EELS spectra in the xz plane subjected to an electron beam propagating along the y axis can be easily obtained from [40,48], − 4e 2 v 2 ˆd tdt ′ Im e iω(t ′ −t) G yy (r e (t) , r e (t ′ ) ; ω) , where v is the speed of electrons, and the single mode expansion for our Green function-that is already confirmed to be very accurate-is used. The EELS calculations for all three nanoparticles (of Fig. 1) are shown in Fig. 4, all computed at the corresponding plasmonic peaks frequencies. Note that there are some noticeable numerical problems around the sharp corners of the metallic nanorod when using the conventional Drude model theory on the left (which is a known problem [60,61]). Using the same meshing scheme, however, the nonlocal description evidently helps to avoid such nonphysical predictions. More importantly, as the nanoparticle size decreases the EELS map becomes brighter at its maximum location which originates from the higher modal amplitudes of the QNMs discussed in Fig. 1. We stress that with the computed QNMs, such EELS maps are calculated instantaneously, which is a far cry from the usual brute force numerical solvers.

D. Field-driven Mollow triplets and quantum optics regime
Finally, in addition to the previous discussions on Purcell factor and mode volume that are essential for building quantum optical models of light-matter interaction [62], we discuss the quantum regime of field-driven Mollow triplets for QDs coupled to plasmonic nanoparticles. In the dipole and rotating wave approximations, the to- tal Hamiltonian of the coupled system can be written as [52,63] where Ω = Ê pump (r d ) · d/ is the effective Rabi field, σ + , σ − are the Pauli operators of the two-level atom (or exciton), ω x is the resonance of the exciton, d is the dipole of the exciton, andf ,f † are the boson field operators. Following the approach in Ref. [52], and using the interaction picture at the laser frequency ω L , one can derive a self-consistent generalized master equation in the 2nd-order Born-Markov approximation: whereJ ph (τ ) =´∞ 0 dωJ ph (ω)e i(ωL−ω)τ , with the photon-reservoir spectral function given by J ph (ω) ≡ d·Im[G(r d ,r d ;ω)]·d π ε0 , and the time-dependent operators are defined through σ ± (−τ ) = e −iHS τ / σ ± e iHS τ / , with H S = (ω x −ω L )σ + σ − + Ω/2(σ + +σ − ), which results in a complex interplay between the values of the local density of states at the field-driven dressed states. Solving Detected spectra (S evaluated at rD = 200 nm) of a field-driven QD coupled to plasmonic nanoparticles, where the same ordering of the particle size and QD location as in Fig. 1 is followed, and we use an effective Rabi field of Ω = 50 meV. The plasmonic enhancement is also shown in dashed-gray in background. Nonlocal investigations on the right predict relatively stronger side peaks for the Mollow triplet with narrower linewidths.
the master equation, and exploiting the quantum regression theorem, one can compute the incoherent spectrum of the QD emission from [52] S 0 (ω) = lim as well as the detected spectrum, which includes quenching and propagation from QD at r 0 to a point detector at r D , from [52] S (r D ; ω) = 2 ε 0 |G(r D , r 0 ; ω) · d| 2 S 0 (ω).
For example calculations, we assume a QD with the dipole moment of |d| = 50 Debye at 10 nm away from the nanoparticle surface, at x=0. In particular, as with the calculations above, we compare the local Drude model versus nonlocal GNOR model as shown in Fig. 5. As can be recognized, including the nonlocal effects, in general, predicts narrower linewidths for the Mollow triplets (see e.g. Table I) where the relative strength of the side peaks are also increased. This is attributed to the modified plasmonic enhancment in the nonlocal description as confirmed in Fig. 3, which can be also rigorously confirmed using analytical equations for the linewdiths derived in h (nm) Drude FWHM (meV) GNOR FWHM (meV) 20 1 [52]. It should be noted that, in general, the detected spectrum S for the Mollow triplet problem can be different than the emitted S 0 as discussed in [52], however in our particular case, under the resonant excitation, we find that they have the same qualitative shape (and differ only quantitatively). We also stress that, these spectral calculations, at any detector position, can be trivially performed using a standard desktop through use of semianalytical GF of (5), demonstrating the power of our approach for carrying out complex problems in quantum optics.

IV. CONCLUSIONS
We have presented an efficient and accurate modal description of the nonlocal response of arbitrarily shaped metallic nanoparticles, using a fully 3D model. We have shown how analytical nonlocal QNMs can be used to accurately construct the system GF from which modal quantities of interest such as Purcell factor and effective mode volume can be derived. As anticipated, we first observe the blueshift as well as the larger damping rate for the LSP as a consequence of nonlocal effects. We further confirmed the validity of our approach for different nanoparticle sizes with full dipole solutions of the modified Maxwell equations, and we were able to successfully predict the size-dependent nonlocal modifications with ease. As example applications of the theory, we described how our nonlocal QNMs can be used to efficiently model Purcell factors of dipole emitters, EELS spatial maps as well as Mollow triplet spectra of field-driven QD. The presented model has many applications in both classical and quantum nanoplasmonics, offers considerable analytical insight into complex nonlocal problems, and could help pave the way for a quantum description of both light and matter.