Quasinormal mode approach to modelling light-emission and propagation in nanoplasmonics

We describe a powerful and intuitive technique for modeling light-matter interactions in classical and quantum nanoplasmonics. Our approach uses a quasinormal mode expansion of the Green function within a metal nanoresonator of arbitrary shape, together with a Dyson equation, to derive an expression for the spontaneous decay rate and far field propagator from dipole oscillators outside resonators. For a single quasinormal mode, at field positions outside the quasi-static coupling regime, we give a closed form solution for the Purcell factor and generalized effective mode volume. We augment this with an analytic expression for the divergent LDOS very near the metal surface, which allows us to derive a simple and highly accurate expression for the electric field outside the metal resonator at distances from a few nanometers to infinity. This intuitive formalism provides an enormous simplification over full numerical calculations and fixes several pending problems in quasinormal mode theory.


Introduction
High-index-contrast dielectric cavities and metallic nano-resonators both facilitate control of light-matter interaction by engineering the local density of optical states (LDOS). While this occurs on the scale of a wavelength for dielectric cavities [1], it can extend down to nm length scales for metallic nanoparticles (MNPs) [2]. Traditional quantum optical effects such as vacuum Rabi oscillations and resonance fluorescence become much richer with MNPs [3][4][5]. Localized surface plasmons dramatically increase the LDOS within a few nm of a MNP, which can enhance the spontaneous emission (SE) rate of a single photon emitter [6][7][8]. The pronounced coupling to surface plasmons and an increased LDOS have been used in a number of MNP configurations, such as single photon transistors [9], chemical detection/imaging [10,11], optical antennas [12], spasing [13,14], enhanced SE [7,12], and long range entanglement [15].
Accurate quantitative modelling of dipole emitters near arbitrarily shaped MNPs is challenging. In many cases, however, the LDOS enhancements in cavity structures are directly attributable to one or just a few local resonances which may be rigorously described mathematically as quasinormal modes (QNMs) of the open system [16]. In particular, it has been shown that the scattering resonance of the MNPs can be described by QNMs [17,18]; so intuitively, it should be possible to capture most of the physics of the dipole-MNP interaction around the local plasmon resonances, using just a few QMNs. Indeed, this is the typical approach in dielectric cavity geometries where, e.g. the Purcell factor [19] for dipole emitters within the cavity can be accurately determined in terms of the cavity mode quality factor (Q) and mode volume (V). However, as recently shown in [20], even in this relatively simple situation, the mode volume is in fact non-trivial to define in a rigorous way because of the leaky nature of the cavity modes which causes the field to diverge far from the cavity. In dealing with MNPs, the situation is further complicated because the dielectric constant of the MNP is complex and strongly dispersive. While some authors believe that it is unnatural to work with modes in a lossy system [21], we argue that QNMs have enormous intuitive appeal in MNP geometries, and that they help to accurately explain the underlying physics of light-matter coupling in a remarkably clear and transparent way. The use of QNMs in the field of nanoplasmonics is made difficult by a number of challenges: (i) techniques developed in quantum optics for lossy inhomogeneous structures suggest that traditional mode expansion techniques do not work [22]; (ii) proper calculations of localized surface plasmons as QNMs are non-trivial; (iii) because QNMs diverge in space, it is impossible to directly use them in calculations of the electromagnetic propagator to positions far away from the resonators; and (iv): the LDOS is known to diverge near a metal surface due to quasi-static coupling (e.g. causing Ohmic heating) so at these positions a single QNM expansion is not expected to work.
In this paper, we describe a new QNM expansion technique that can be used to evaluate the electric field from a dipole emitter at positions both nearby and far away from MNPs (i.e. outside the resonator). We first compute the effective mode volumes and Purcell factors for metal nanorods, and compare the resulting approximate results to full numerical calculations. We elaborate on the concept of mode volumes for the QNMs of MNPs [23] and show explicitly that a recently introduced mode volume for MNPs [18] is exactly the generalization of the result in [20] to the case of dispersive materials. We present a solution to the problem of using QNM expansions in scattering calculations at positions outside and far away from the resonator where the modes diverge [24,25]. Further, we introduce an analytic technique to account for Ohmic losses and quasi-static coupling which become important when the emitter is within just a few nm of a MNP.

Green function expansion in terms of normalized QNMs
Consider a general, non-magnetic, inhomogeneous medium described by a complex permittivity ε ω . For a three-dimensional (3D) geometry, the total photon Green function satisfies the equation where ω = k c 0 and 1 is the unit dyadic. To describe quantum optical effects in lossy inhomogeneous structures [22,[26][27][28], two key Green functions are required, namely G r r ( , ) a a and G r r ( , ) a b ; the former can describe effects such as modified SE or the Lamb shift [28], while the latter includes the effects of photon propagation. Many observables in quantum optics, e.g. dipole-emitted spectra, can be described with the help of these two quantities [4,5,22]. For example, if we consider a dipole emitter at position r a with a dipole moment = d d n a (n a is a unit vector), the relative SE emission rate is [29] [29]. The Green functions may be computed in a number of ways [6,17,22], but in general they are rather expensive to calculate.
Localized surface plasmons may be directly understood as QNMs of the MNPs [17], defined as the frequency domain solutions to the wave equation with open boundary conditions [24,30] (the Silver-Müller radiation condition [31]). This makes QNMs the natural starting point for theoretical developments, although the radiation condition is not immediately compatible with typical mode solvers. For this reason, a common approach is the use of coordinate transforms, typically in the form of perfectly matched layers (PMLs), to model a system with no reflections from the simulation domain boundaries. The QNMs, μ f r ( ), have a discrete spectrum of complex resonance frequencies, ω ω γ = − μ μ μ i , from which the resonator quality factor is ω γ = μ μ Q 2 . An important consequence of the complex resonance frequency is that the QNMs diverge (exponentially) in the limit → ∞ r . In certain spatial regions, such as inside the resonator [24], the transverse part of the Green function can be expanded as [30] ω . This has been rigorously proven for spheres, but is also expected to be true for non-spherical scattering objects that have an abrupt discontinuity in the dielectric constant profile [30]. Thus for a single QNM, considering points inside the resonator, we define where μ f is normalized through [24] ∫ ∫ ω ε ω ω As a representative example, we first study a 2D MNP using a Drude model, We consider a metal rod with width (x axis) 10 nm, and length (y axis) 80 nm located in a homogeneous space with . To compute the QNM profile, f r ( ) c , we use the finite-difference timedomain (FDTD) technique [32] with PML boundary conditions and a run-time Fourier transform [33]. Specifically, the system is excited using a pulsed plane wave (incident along x, and polarized in y, which was chosen to maximally excite the QNM of interest) and a temporal window function is applied for computing the field. Figure 1 shows the spatial profiles of the QNM field as well as the scattered part of the field evaluated at the real part of the complex eigenfrequency. Only the QNM shows the expected spatial divergence [20]. Later we will introduce a regularized QNM field that coincides very well with the scattered field profile, as shown in figure 1(b) (dashed curve); this regularization is essential for modelling light-matter interactions outside the MNP.

Effective mode volume and Purcell factors for spatial locations beyond the quasistatic coupling regime
For most problems in light-matter interactions, one must describe how light behaves far away from the MNP. Examples include calculations of the spectrum that would be detected at some arbitrary spatial location [4,5,22], and the design of nanoantennas [29]. The QNMs diverge at large distances, and therefore it is impossible to use them directly for such calculations. However, since we do have a highly accurate approximation for the Green function inside the resonator, we can exploit a Dyson equation of the form ) to obtain a regularized Green function in the region far away from the scattering geometry (see appendix A). With this approach, we derive a corrected Green function expression for positions outside the resonator, where we have introduced a new regularized field, V B As shown in figure 1(b), this field coincides well with the scattered field because it is essentially the field that is scattered by the QNM at the (real) frequency ω c . In equation (6), we have neglected the influence of other mode contributions as we are currently interested in positions sufficiently far from the MNP, that these can be safely ignored. We now consider an oscillating dipole with (real) resonant frequency near the (complex) frequency of the QNM in figure 1, so that a single mode expansion (μ = c) of the Green function might be expected to provide a good approximation, at least for dipole-MNP separations beyond the quasi-static coupling regime. The next higher order plasmon mode is at least 500 meV away, so we can safely assume a single mode approximation. To describe the relative SE rate in terms of the Purcell factor, we write a a a a P where the Purcell factor has the familiar form 4  in agreement with previous generalized mode volume results for a lossless dielectric cavity [20]. We note that Sauvan et al [18] recently used the Lorentz reciprocity theorem to derive an expression for the effective mode volume which can be shown to be identical to equation (10) (see appendix C). Figures 2(a) and (b) show, respectively, the QNM spatial profile in more detail and the evaluation of V eff as the size of the calculation domain is increased. Although each term in equation (4) diverges as a function of domain size, we find that for the MNPs studied in this paper, the sum converges quickly to a finite value after ∼500 nm. We note, however, that in principle the convergence is non-trivial because of a small phase differences between the two terms, which may eventually show up as oscillations at relatively large distances from the resonator [20]. Although we have never found it to be a problem in practice, the issue can be handled by a regularization of the integral, for example by use of coordinate transforms [18].
In optical fiber geometries, this crossover region is referred to as the 'caustic radius' (r caustic ) [35]. Although one can define the SE rate in terms of the Purcell factor as in equation (8), it may in practice be more convenient to work directly with equations (2) and (5) as we do below. Figure 2(c) shows the relative SE rate as a function of frequency for a dipole aligned along the direction of the nanorod axis and located just above the top interface. The grey dashed line is the result from a full, independent numerical calculation using FDTD [6,32], i.e. with no mode expansions (namely, a full dipole calculation), and the green solid is the result from equation (8). Clearly, the analytic SE rate gives a very good fit to the full dipole calculation at this spatial location, including the entire non-Lorentzian lineshape. Figure 2 , as a function of distance (5-500 nm) from the MNP. For these distances, the single QNM expansion in equation (6) provides an excellent approximation to the full dipole results, even if one uses μ f˜in place of μ F (see figures 1(c) and (d) over this range). However, the approximation must fail at evaluating the SE rate for larger emitter-MNP separations, and for evaluating the propagated fields, since the rigorous QNMs, μ f˜, diverge in space. In figure 3(a) we show the on-resonance SE factor as a function of distance, spanning 5 nm-5 μm using G f in place of G F . Comparing to the full dipole solution, the single QNM approximation clearly fails to get the correct far field behavior, as anticipated from our discussions above. In contrast, the SE factor computed using G F in G far (see figure 3(b)) shows excellent agreement with the full dipole results and yields the correct behavior for large distances where it tends to unity. Importantly, in both cases we only include μ = c in the Green function summation, although G F accounts for the influence of light propagation. In figure 3(c) we show the two model results as a function of frequency at μ = r (5, 0) m, and G f again yields incorrect results. In figure 3(d), we show the absolute square of the propagator, ω G r r | ( , ; )| yy b a c 2 , as a function of r b , with selected full-dipole results shown with the circles; once again, we observe that only G F gives the correct behavior at large distances.

Enhanced SE factors for spatial locations very near the MNP resonator
Together, figures 2(d) and 3 demonstrate that the regularized QNM fields μ F , through G far , provide an accurate representation of G for all separations down to at least ∼5 nm. However, ultimately G far must fail at even shorter dipole-MNP separations because of the known quasistatic divergence of the LDOS. This divergent behavior cannot be accurately accounted for in a single QNM approximation. Our solution to this problem is to note that the dynamics of a dipole near a MNP surface is governed by essentially the same LDOS increase as a dipole near a metal half space. Moreover, this divergent term vastly dominates the LDOS, and therefore, at spatial points very near the resonator, we may simply add on this known quasi-static Green function to the single QNM approximation to the Green function, so that ω ω ω = + G r r G r r G r r ( , ; ) ( , ; ) ( , ; ), where the quasi-static term is given by [ , and ∓ is for s-or p-polarized dipoles, respectively. We find this formula to be quantitatively accurate for all positions that we have tried. For further details, see appendix A.
To demonstrate the accuracy of equation (11), we now consider a 3D metal nanorod (i.e. a cylinder) with radius 15 nm and length 100 nm, where the rotational axis coincides with the y axis (see figure 4(a)). The Drude parameters are the same as before, but we use γ = × − 1.41 10 rad s 14 1 to model gold. We find a frequency region dominated by a single QNM at ω π (2 ) c = 324.981−i16.584 THz (1.344−i0.0684 eV), and the next higher order mode is at least 600 meV away. The mode profile is shown in figure 4(a), and V eff is shown in figure 4(b). In figure 4(c) we show the SE spectrum for a dipole emitter at 2 nm (see arrow in figure 4(a)). The green/lower dashed curve uses equation (5) (with G f in place of G F ), while the grey/upper dashed curve is the full dipole result which is seen to be considerably different; in contrast, the proposed G out model (solid green), through equation (11), is seen to be in quantitative agreement with the full dipole result even in this spatial regime dominated by non-radiative coupling. In figure 4(d), we also show the propagator and confirm an excellent fit as a function of distance only if G F is used.

Conclusions
We have introduced a QNM expansion technique for arbitrarily shaped MNPs and numerically demonstrated its accuracy in evaluating the enhanced SE rate and far-field propagator associated with proximate dipole oscillators. In the limit of a single QNM resonance and spatial positions outside the regime of quasi-static coupling, we give a closed form solution for the Purcell factor. The corresponding effective mode volume is found to be identical to the effective mode volume recently introduced for dielectric cavities [20], with a generalization to account for material dispersion. Our results are obtained using a transparent analytic expression for the photon Green function evaluated outside of the MNP, qs , that is valid everywhere outside the particleʼs boundary. The G qs term is a simple, intuitive quasi-static dipole expression that accounts for quasi-static coupling within a few nanometers of the metal resonator, while F is the sum of the familiar Green function of a homogenous medium, and a resonant mode expansion. Our general approach can be applied to a wide variety of resonator geometries including hybrid metal-photonic-crystal structures [37] and coupled QD structures [38,39], which are otherwise very hard to understand and model.  It was shown in [24] for a one-dimensional lossy cavity that if the permittivity or any order of its derivative is discontinuous at the border of the cavity, then the QNMs inside the cavity form a complete basis; so one can use a mode expansion formulation for the Green function in terms of the QNMs. Later, [30] proved the same result for a 3D sphere, and they also discussed why the same argument could be applied to structures without spherical symmetry. Bergman and Stroud [40] used the same argument for scattering geometries made up of spherical structures. For our calculation of nanorods, we also find that a single QNM expansion works exceptionally well, and, moreover, we formally recover the same result as in [18] (which uses an entirely different method). Consequently we can assume that the mode expansion formulation approach is valid for points inside resonators of any shape of geometry as long as the permittivity (or its derivative) is discontinuous at the border of the resonator. To obtain the corrected Green function for points outside the scattering geometry (MNP), we utilize the Dyson equation MNP B , with ε B and G B the dielectric constant and Green function of the homogeneous background in which the MNP is embedded. For ′ r 1 inside the scattering geometry Substituting the QNM expansion of ω ′ ′ G r r ( , ; ) 2 1 , for which both points ′ r 2 and ′ r 1 are inside the MNP, into equation (A.2), and focusing on the single QNM of interest, we obtain in which r 1 and r 2 are two space points outside the resonator and we have introduced a new mode field, defined as In equation (A.4), G others formally includes the contributions from all other QNMs except μ f r ( ) which is typically by far the most dominant field at distances far enough away from the MNP that quasi-static effects may be neglected. As we show in the main text, equation (A.4) provides an excellent approximation to G at most points outside the resonator, when the G others term is neglected, apart from very close to the metal surface where it only fails for ⩽ x 5 nm. These two Green functions (equations (A.6) and (A.7)), which are both transverse, are used in our manuscript; and we demonstrate that only G F gets the correct far-field behaviour outside the effective mode volume region of the QNM (⩾ ∼ 500 nm), while G f in place of G F yields divergent propagators and divergent or/and negative enhanced emission rates for oscillating dipoles in this limit. The additional term in equation (A.4), namely G B , also explains exactly why the extra factor of unity is required to properly relate the relative SE rate to the Purcell factor; interestingly, we note that such a factor does not appear for spatial regimes inside the scattering geometry, such as with emitters inside photonic crystal cavities for which equation (A.7) is the appropriate single-mode approximation that can be used to derive the Purcell factor [20]. We stress that the only approximation is to use the QNM transverse Green function expansion within the MNP, with a Dyson equation theory that obtains the correct transverse Green function solution outside the scattering geometry.
Neither G far or G f are suitable for the extreme near field regime (i.e. a few nm from the MNP surface); however, the dominant response of a dipole near a metal surface is very similar to the behaviour of a dipole near a metal half space, and this exact quasi-static response is known analytically, through [36] ω where ∓ is for s-polarized or p-polarized dipoles, respectively. Moreover, this term vastly dominates the response at these distances, and therefore we propose to simply add this term in place of G others . Thus, in total we arrive at the approximate analytical Green function for use at all positions outside the MNP This expression is both simple and physically appealing in how it describes the various physical processes that occur beyond a single QNM expansion; yet, all that is required is a single QNM supplemented by the normalization procedure for a generalized effective mode volume and quasi-static coupling. In fact our formalism also makes it clear how to simplify complicated numerical calculations requiring large memory and small grids outside the MNP. Since all that is required is the QNM within the MNP, then the Green function can be computed to any arbitrary spatial point outside the MNP analytically. In figure A1, we show the predicted SE factor for a quantum dipole emitter oriented along the 3D MNP (gold nanorod with a Drude permittivity model-see main manuscript for more details) using the various Green function models. In figure A1(a) we show that calculations using G f (in place of G F in G far ) produce essentially identical results for the distances considered; and since G f is easier to compute, one can safely use G f for distances less than several hundred nm (see figure 2(d)). However, as we show explicitly in figure 3(a), for distances outside the caustic regime, ⩾ ∼ r 500 nm, the SE rates predicted by G f in place of G F grow exponentially (which is clearly unphysical), while G far consistently produces results that are in very good agreement with full dipole calculations.
To elaborate on the approximation made to equation (A.4), we note that G others includes also a first-order background term which is negligible above about 40 nm, and finite below this value (black solid curve). This term, however, is only a first-order Born approximation and turns out to be not reliable. To make this clear, in figure A1(b), we compare results with this term against the proposed G qs (i.e. the nonperturbative half space solution), and only the latter is seen to be in agreement with full dipole FDTD results [32]; indeed, the agreement is quantitative. Although we have chosen a single frequency here, in the main text we also show excellent agreement as a function of frequency at the spatial location of only 2 nm. It is interesting to note that Sauvan et al [18] also found that a single QNM expansion works well for evaluating the SE rate at a range of spatial distances outside the MNP using a dipole emitter position-dependent effective mode volume in the Purcell formula, and adding to it a factor of unity, seemingly 'by hand.' The formalism presented here: (i) defines a mode volume that is characteristic of the surface plasmon mode alone, and a separate factor η that characterizes the dipole coupling to that mode, (ii) provides an explanation for the factor of unity, (iii) provides an efficient and intuitive method for accurately calculating both SE rates and field propagation effects to arbitrarily large distances from the MNP, and (iv) also works accurately at short distances where quasistatic coupling is important.

Appendix B. Peak SE deviation factor
In the main text, we defined the relative SE rate in terms of the Purcell factor For spatial points near the resonator, one can simply replace μ F r ( ) a with μ f r ( ) a . Note, to derive these expressions we have used the Green function of a 3D homogeneous medium; expressions for G B for 2D and 3D are given in, e.g. [41] and [29], respectively, so the 2D expressions can be derived in the same way.

Appendix C. Effective mode volumes: correspondence with recently published results
In a recent paper, [18] Sauvan et al address the problem of normalization of QNMs in dispersive resonators. In particular, the authors use the Lorentz reciprocity theorem to derive a generalized mode volume for use in the Purcell formula. Below, we show that the normalization in [18] is a generalization of a previous result by Leung et al [24] to the case of vector fields and magnetic materials. For non-magnetic, dispersionless materials, the corresponding mode volume reduces to the same generalized mode volume that was previously introduced for leaky optical cavities [20]. Sauvan