Eigen-energy effects and non-orthogonality in the quasi-normal mode expansion of Maxwell equations

We derive a quasi-normal mode theory for three-dimensional scatterers, taking care to remove an hypothesis of weakly dispersive materials implicitely used in previous works. In our approach, the normalized modes remain unchanged, but the analytic expansion coefficients onto the set of QNM are modified. In particular, we take into account in a simple way the non-orthogonality of the modes, and we set up a rigourous frame, to treat the case where several QNMs are excited. Eventally, the complex concept of PML integration, previously introduced, becomes unnecessary, even to compute the QNM mode volume. Besides, crossover integrals of QNM fields over the whole space can now be written rigourously, as integrals on the finite volume of the scatterer, without surface terms. OCIS codes: (290.5825) Scattering theory; (290.2200) Extinction; (260.5740) Resonance. References and links 1. D. J. Bergman and M. I. Stockman, “Surface plasmon amplification by stimulated emission of radiation: quantum generation of coherent surface plasmons in nanosystems,” Phys. Rev. Lett. 90(2), 027402 (2003). 2. P. Berini and I. De Leon, “Surface plasmon-polariton amplifiers, and lasers,” Nat. Photonics 6(1), 16–24 (2011). 3. R. D. Artuso and G. W. Bryant, “Optical response of strongly coupled quantum dot-metal nanoparticle systems: double peaked Fano structure and bistability,” Nano Lett. 8(7), 2106–2111 (2008). 4. A. Manjavacas, F. J. García de Abajo, and P. Nordlander, “Quantum plexcitonics: strongly interacting plasmons and excitons,” Nano Lett. 11(6), 2318–2323 (2011). 5. J. Yang, M. Perrin, and P. Lalanne, “Analytical formalism for the interaction of two-level quantum systems with metal nanoresonators,” Phys. Rev. X 5(2), 021008 (2015). 6. L. Hayati, C. Lane, B. Barbiellini, A. Bansil, and H. Mosallaei, “Self-consistent scheme for optical response of large hybrid networks of semiconductor quantum dots and plasmonic metal nanoparticles,” Phys. Rev. B 93(24), 245411 (2016). 7. J. N. Anker, W. P. Hall, O. Lyandres, N. C. Shah, J. Zhao, and R. P. Van Duyne, “Biosensing with plasmonic nanosensors,” Nat. Mater. 7(6), 442–453 (2008). 8. J. Yang, H. Giessen, and P. Lalanne, “Simple analytical expression for the peak-frequency shifts of plasmonic resonances for sensing,” Nano Lett. 15(5), 3439–3444 (2015). 9. T. Weiss, M. Mesch, M. Schäferling, H. Giessen, W. Langbein, and E. A. Muljarov, “From dark to bright: firstorder perturbation theory with analytical mode normalization for plasmonic nanoantenna arrays applied to refractive index sensing,” Phys. Rev. Lett. 116(23), 237401 (2016). 10. G. D. Bernasconi, J. Butet, and O. J. F. Martin, “Mode analysis of second-harmonic generation in plasmonic anostructures,” J. Opt. Soc. Am. B 33(4), 768–779 (2016). 11. M. Perrin, J. Yang, and P. Lalanne, “Analytical treatment of the interaction between light, plasmonic and quantum resonances: quasi-normal mode expansion,” Proc. SPIE 9755, 97551J (2016). 12. 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). 13. 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). 14. 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). 15. P. T. Kristensen, C. Van Vlack, and S. Hughes, “Generalized effective mode volume for leaky optical cavities,” Opt. Lett. 37(10), 1649–1651 (2012). 16. 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).


Introduction
Nanoplasmonics is a blooming field where quantum, classical and non-linear physics meet, leading to promising applications, such as, new light sources [1,2], artificial materials [3][4][5][6], sensors [7][8][9], devices for frequency convertion [10].In all these domains, researchers are dealing with complex systems, which are dispersive, dissipative and often three dimensional.Obviously, there is an acute need for design rules and simple tools that would ease their understanding.This cannot be entirely provided by brute force computations, that can be lengthy, if one wants to determine in detail the optical response of a device (every angle of incidence, frequency, polarization of a plane wave excitation).
In this quest, the Quasi -Normal Mode expansion seems promising for many purposes both in quantum and classical physics [11].This theory was first proposed for non-dispersive, uni-dimensional systems [12], and then extended to three-dimensional, dispersive systems of arbitrary shape [13][14][15][16], to periodic non-dispersive [17], and dispersive structures [9,18], as well as to waveguides [19,20].Let us also mention that QNM theory can be helpful in the perturbative limit [8,11], as well as the Resonant State Expansion method [21,22], that permits to obtain the excitation coefficients for a perturbed system, using the basis of an unperturbed one.
The mathematical foundation of QNM theory for 3D dispersive material is complex, and takes its root in fundamental works on the resolution of generalized nonlinear eigenvalue problems for non-self-adjoint opertors [23][24][25].There, the electromagnetic field is developed as a function of a complex frequency in the vicinity of the poles of the scattering operator.Then, assuming no essential singularities are present, it can be casted under a Laurent serie [25], with two contributions: a resonant one -due to the poles, and described by an expansion on QNMs -, and a non-resonant one, given by an unknown [13], but holomorphic [25] function of frequency.In the available studies, the knowledge of the resonant contribution has proven to be sufficient to obtain a good accuracy, but it should not always be the case.
Once the field is expanded onto a set of eigenmodes [13,14], many important physical quantities, such as the Purcell factor [13,14] and cross sections [14] can be computed in a semi-analytical way.Let us mention as well applications to Electron Energy Loss Spectroscopy [26], and problems of quantum optics [5,27].It has been shown that all these physical quantities -and even the Purcell factor -can be computed with the sole knowledge of the QNM field inside the nanoparticle (NP), what helps circumventing the open problem of the QNM divergence far from the scatterer.
However, a serious difficulty in QNM theory concerns its normalization.A priori, it requires the integration over the whole space of a quantity which oscillates and diverges far from the particule [16].Even if the QNM norm (the result of the integral) is finite, this is a difficult task to regularize the integral to compute it.Two main methods, proved to be equivalent [16], can be followed to evaluate the norm by direct integration.First, one can evaluate the integral on a finite volume, taking care to consider a boundary (finite surface) integral as well.Second, one can use the fact the integration in the perfectly matched layers (PML) that surrounds the particule permits to regularize the integral [13,16].In a much simpler approach [14], the value of the norm is not computed but simply set to a constant.This leads to a third way to normalize QNMs, and to analytical expressions of the expansion coefficients onto the set of QNMs [13,14].The price to pay, is simply that one now also need to compute the field at one frequency in the vicinity of the QNM eigen-frequency, so as to normalize the mode [14] (see in particular Eq. ( 8) in [14]).
One other difficulty in QNM theory lies in the fact the modes are not orthogonal [13] in dispersive media.Examples in other fields of physics -surface plasmonics [28,29], or cavity optics [30] -teaches us that non-orthogonality can be responsible for strange behaviors, e.g. the total energy absorbed or scattered is not the sum of the energy absorbed [28] or scattered by each mode.Despite this non-orthogonality has already been pointed out [13], no simple way has been proposed until now (to our knowledge) to handle this effect -which has been largely overlooked -, except a complicated integration procedure inside the PML.This method however requires a good knowledge of PML modeling, and cannot be used if PML are absent [31].Note that it might also be possible to use the volume and surface integration method [16] to handle non-orthogonality.
Besides, in previous works, a Taylor expansion of permittivity has been used in the beginning of the derivation of QNM expansion coefficients [13,14].We shall show that this hypothesis can be removed, what improves the accuracy of the method, and renders the PML integration unnecessary, even if one wants to model rigorously a set of non-orthogonal QNMs, and compute rapidly the QNM mode volume.
In section 2, we set up the scattering problem, and detail the two possible point of view to define the modes -either with two first order either with one second order Maxwell equation.In section 3, we present a new normalization, which involves only the field inside the nanoresonator and show how this can be used to compute simply the QNM mode volume.We compare the obtained results with canonical examples from the literature.Then, we present a technique to express the integrals of QNM fields on the whole space into integrals on the sole volume of the nanoparticle(s), without surface terms.This method is general, and we use it to derive an analytical expression for the expansion coefficient that is more rigourous and accurate than what done in previous works [13,14].Several numerical tests are presented.In the fourth section, we generalize this method to a set of several QNMs, and obtain a closed form analytical expression to compute the expansion coefficients, that takes into account rigorously the non-orthogonality of the modes.

Formulation of the scattering problem
In the following, we consider monochromatic fieds, and use the exp( ) Two ways are commonly used to set the modal problem.This can be through a first order representation [13,14]: ( , ) , and ( , ) , or a second order representation [16,21,22]: for non-magnetic materials.In Eqs. ( 1)-( 2), ( , ) ε ω r is the permittivity, that depends on the frequency -inside the nanoparticule -and is assumed to be constant outside the nanoparticle, surrounded by a medium of permittivity -unless otherwise stated.These local equations are complemented by the condition that the field should obey outgoing boundary conditions.
In the following, we assume that one or several nanoparticle(s) scatters an exciting field ( , ) r , and we write the scattered field at frequency ω, ( , )   ω s E r , as a sum of QNM, plus a non-resonant term [13,14], that can be assumed to be an holomorphic function of frequency [25]: Note that, the QNM field of eigenvector at a large distance r from the scatterer, whereas the scattered field should decay as 1/ r .Therefore, Eq. (3) is not valid in the whole space, but in some vicinity of the NP only, whose size would depend on Im( ) k  -for the QNM of lowest frequency of the NP we study, this size is about one wavelength, see [32].However, this point does not raise much problem in practice, as: (i) we shall show that the modes can be normalized with the sole knowledge of the field inside the NP.(ii) most of the physical quantities can be expressed as functions of the QNM field and excitation field inside the NP [14].Note that a good agreement between QNM expansion and exact results has been obtained for the Purcell factor of a dipole placed far from the NP, when the scattered field is expanded [14].
The scattered electric field excited in the system verifies where ( , ) r is the current source for the scattered field [14], and ε ∆ represents the permittivity change that models the nanoparticule(s).
Our goal is to find the expansion coefficients ( ) m β ω that will give the best approximation to the scattered field, so that we will consider ( , ) 0 . Then, once the ( ) m β ω are obtained, the absorption and extinction cross sections can be respectively computed as [14]: and Additionally, the Purcell factor of an electric dipole 0 J located at 0 r is obtained [13,14] as: [ ] where

Field normalisation
In this section, we assume only one QNM is excited.
The QNM normalization can performed by direct integration on the whole space of a function (depending on QNM fields) which is actually computed on a finite volume [13,16].If one were to compute such integral it could be useful to use the PML integration scheme proposed in [13], which is based on the fact that the integrand is invariant under PML coordinate stretching.But a direct calculation where integration is performed on a finite volume V also give the same result, provided a surface term (or conjunct integral, see Eq. ( 5).10, chapter 5.1 of [23])-is considered as well [16].In the following we denote the integral over the whole space, computed with PML integration as 3 ∫  to simplify the notation.
Instead of computing directly this integral, we impose a normalization condition, that fixes the value of the norm.For the first order scattering problem, this condition reads [14]: and comes naturally from Lorentz reciprocity [13,14].As shown previously [16,31], the integrand of the left hand side (l h.s) in Eq. ( 8) can be expressed as a function of the electric field only, using the equality . Then, using the vectorial identity A B , one obtains an equivalent normalization condition: Consequently, the QNM norm, as defined by [16] now reads: In order to define ( )

E r 
, let us first expand the scattered field in the vicinity of the pole, for the second order problem Eq. ( 2).We apply the divergence theorem (or Green -Ostrogradsky theorem) to the vectorial identity [16]: Then, taking ( , ) and ( ) , and integrating on the whole space (including PML), one eventually gets, using Eq. ( 4): Lets now consider a frequency Ω close to the pole, ω  , so that a Taylor expansion of the permittivity is valid, and Eq. ( 12) becomes: When ω Ω →  , the scattered field increases dramatically [14] and one can neglect all the contributions except the ω  pole in the QNM expansion, to write ( , ) ( ) ( ) . Multiplying both sides of the later equation by r and integrating, one gets: Inserting Eq. ( 14) in Eq. ( 15) and using that ω Ω ≈  , one gets: , and eventually one obtains the QNM definition: where Ω is a complex frequency, arbitrarily chosen, close to the pole, and the integral is carried out only on the nanoparticle volume.We checked that the normalization does not depend on the way Ω converges towards ω  .Besides, using Eq. ( 17) to define ( )

E r 
, we obtain the same results than in previous works [13,14] for the QNM profile, cross sections and Purcell factor.
Equation ( 17) permits to define the QNMs when any type of source is exciting the mode, what extends the specific normalization of [14], valid for a dipole source excitation, and now uses the scattered field only inside the nanoparticle.
This definition of QNMs is valid for any geometry, but is particularly useful for structures with cylindrical symmetry, excited, e.g., by a loop of current.Then, the q th azimuthal modes, whose eigen-fields vary in exp(iq ) ϕ are found, for a specified value of the integer q.Note that we used a finite element method (FEM) [14] with an axi-symmetric model to make the calculations.The computation time is decreased by one or two order of magnitude, compared to three dimensional calculations with same mesh parameters [14].In practice, the maximum size of the mesh elements is 2 nm in the metal and 10 nm outside (including in PML), the spherical PML internal and external radius where set respectively to 0.4 µm and 0.6 µm.
In order to compare this normalization to others, we compute the mode volume of QNMs in different canonical systems.Using the definition of [13], and taking advantage of Eq. ( 8), the mode volume for Purcell factor has a simple expression, using the present normalization (see also [34] for Mie resonators): Where 0 r is a specific point around the NP, and u the direction of the dipole moment of the point source.Note that one does not need to integrate in a PML to compute the mode volume Eq. ( 18) [11].
In Table 1, we compare the results from Eq. ( 18), computed with our normalization Eq. ( 17), to results from the literature.For meaningfull comparison, all the data have been converted if necessary into the exp( ) i t ω convention we use for harmonic fields.Nanorod A: 100nm length and 15nm radius gold nanorod in glass (n b = 1.5).0 r is placed on axis, 10 nm above the rod, u is oriented along the rod axis.Dimer B is made of two coaxial gold nanorods in glass.Their diameter are D1 = 20 nm and D2 = 85 nm, and their length are L1 = 80 nm and L2 = 145 nm; they are separated by a gap g = 45 nm.For both A and B, the mode volume is expressed in units of 3 4 / 10 where λ is the (real) resonance wavelength, and has been computed with Eq. (18).Dimer C is made of two identical and coaxial rods of 100 nm length and r = 15 nm radius, separated by a 20 nm gap.For both dimers, the dipole is located at the gap center, and oriented along the rod axis.
We find a good agreement between our approach Eqs. ( 17)-( 18) and previously published results, relying either on PML [13] or on a direct volume integration [16].The small differences between the predictions might be explained by the fact that different numerical methods have been used (FEM -present work -, FDTD [16], aperiodic Fourier modal method [13]).
Importantly, PML integration is not needed anymore to compute the mode volume if one uses Eq. ( 17) to define the QNM.

Conversion of PML integrals to finite volume integrals
Equation ( 14) raises upon a Taylor expansion, Eq. ( 13), so that one can question its validity when the excitation frequency gets farther from the pole frequency.
To obtain a more accurate and rigorous expansion than was done previously [13,14], one would need, at first sight, to do a PML-integration to compute in the l.h.s of Eq. (12).
To circumvent this difficulty, one cuts the l.h.s integral in Eq. ( 9) in two parts: an integral over the volume of the nanoparticle(s) and an integral over the background media, outside nanoparticle(s).This is the originality of our method.One obtains: where the symbol 3 V   means that integration is done on the whole space, including PML, except the nanoparticule volume, V , a region where permittivity is non dispersive, but can a priori depend on position: ( ) b ε r .The very important point in Eq. ( 19) is that it permits to connect an integral on 3   (that could be computed using PML domain) to a "regular" integral on physical fields in a finite volume, obtained under a normalization condition Eq. ( 9).

Derivation of expansion coefficients
Equation ( 19) indeed permits to express rigourously the 3   integral in Eq. ( 12) in terms of integral on a finite volume.Finally, the expansion coefficient is expressed without using a Taylor expansion of the permittivity, and only as a function of integrals inside the volume of the nano object(s).
The expansion coefficicent for the second order problem is: Note that Eq. (19), is also valid to model a nanoparticle in a non-homogeneous background medium -e.g. a glass substrate, a stratified dielectric medium -, and so is Eq. ( 21).The same derivation can be used when the nanoparticle is surrounded by a non-homogeneous and dispersive medium -e.g. a metallic substrate -, provided that its permittivity is separable and can be casted under the form ( , ) ( ). ( ) . Then, a term depending on dh dω will appear in the expansion coefficient.
Redoing the same calculation, for the first order description, Eq. (1), using Eq. ( 8) instead of Eq. ( 9), one obtains the corrected expansion of the scattered field, at first order: Compared to previous theories [13,14], we notice the presence of an additive term in the denominator.For both expansions, this term is homogeneous to an energy, as indeed represents the electric energy of the mode stored inside the dispersive scatterer.The main difference by respect to standard definition of energy in dispersive media [35] is that energy now becomes a complex number, whose phase directly influences the real and imaginary part of the expansion coefficient.We denote this term: the eigen energy correction of QNM expansion.Finally, note that when ω is close enough to ω  , both the first and second order expression reduce to the previously derived relations for the analytic expansion coefficient [14].

Numerical study
In the following numerical tests, we consider one or two nanoparticle(s) made of gold and surrounded by a glass background of permittivity

Cross sections
Comparing numerical simulations made on the very same system than in [14] shows that the expansion coefficient we propose gives a better accuracy.Figure 1 shows the extinction and absorption cross section, computed from Eq. ( 5) and Eq. ( 6), where different models have been used to evaluate ( ) β ω .The nanoparticle is a cylinder of gold in glass, of length 100 nm and diameter 30 nm, excited by a plane wave polarized along the rod axis and propagating perpendicular to the rod axis (see inset).Magnifications are shown in inset, where an horizontal dot-dashed line helps to distinguish between positive and negative cross sections.The difference between first and second order expansion is of the order of the accuracy of the numerical method.
On the large wavelength aisle of the resonance, extinction cross section computed with previous models [14] becomes frankly negative.This problem is still present, but much less noticeable -see Fig. 1(a) -if one uses the new expansion coefficient Eq. ( 22), that actually gives a better estimation of the scattered field.Note that in this example, both first and second order give almost the same result.
To provide another example, we study a longer rod [11], for which we found three QNMs between 500nm and 2500 nm, at wavelength .We treat them as if they were orthogonal, assuming that the excitation coefficients can be computed using Eqs.( 21) - (22).It turns out that the mode at 2 λ is antisymmetric and cannot be excited by a plane wave polarized along the rod axis [11].Eventually, the results shown Fig. 2(a) display two resonances in the extinction spectrum.Apart from the improvement at long wavelength, see Fig. 2(c), we find that extinction at resonance is also described more accurately by the new model, with an improvement around 2%, see Fig. 2(b).We also note that the QNM method gives a good prediction of the weak resonance around 0.6 µm.Nevertheless, further investigation at smaller wavelength would necessitate to use a more elaborate dispersion relation for permittivity [36], and goes beyond the scope of this article.Fig. 2. Comparison of the different modal expansions for two QNMs treated as if they were orthogonal.The nanoparticle is a cylinder of gold in glass, of length 200 nm and diameter 30 nm, excited by a plane wave polarized along the rod axis and propagating perpendicular to the rod axis [11].Panel (a) shows the extinction cross section over a large spectrum, where two of the three QNMs are excited by the incoming plane wave.Panel (b) and (c) show an improvement respectively around the main resonance, and the long wavelength tail, when one uses our model (both first and second order give practically the same results, only the second order has been shown).An horizontal dot-dashed line helps to distinguish between positive and negative cross sections.Panel (d) shows the short wavelength extinction peak, where both the present and the former [14] models give the same result, slightly different from exact calculations.Legends are the same as in Fig. 1.

Purcell factor
The Purcell factor of a dipole placed in the vicinity of a NP can also be evaluated from the QNM theory, expanding either the total field [13] or the scattered field [14] on the set of QNM.
An important difference in practice, is that the scattered field expansion, e.g.Equation ( 22), necessitates to compute crossover integrals between the field emitted by the dipole alone in the background, and the QNM, , at each frequency, whereas the total field expansion simply necessitates to evaluate the scalar product between the dipole current 0 J and the QNM field at the dipole location, 0 ( ) E r  , which has a constant value.This later is therefore more convenient to use, eventhough the field scattered by a dipole in an homogeneous media is known analytically [33].Following the derivation in [13], and using the method shown in section 3.1, one obtains: where the total field is well approximated by ( , ) . One should note that the very same eigen-energy corrective term is present in the denominator of Eq. ( 23) -for the total field expansion -and Eq. ( 22) -for the scattered field expansion.For the second order total field expansion, the correcting term takes the form of Eq. ( 21).
One can then use Eq. ( 7) to compute the Purcell factor, expanding either the total field, or the scattered field.At the difference of Fig. 2, the three QNMs can now be excited by the dipole near field, and a third resonance appears in the Purcell spectrum that was absent from cross sections spectra -see Fig. 3(a), and the discussion in [11].Looking at the long wavelength aisle, see Fig. 3(b), one observes that predictions are improved when one uses Eq. ( 23) to expand the total field, compared to former models [13,14].One also observes a fairly good agreement between exact results and those obtained by expanding the scattered field with Eq. ( 22), see Fig. 3(b).On the contrary, the computation of Purcell factor from scattered field, using the former model [14] leads to a very poor agreement, and even a negative Purcell factor, far enough from the resonance, see Fig. 3(a).Fig. 3. Comparison of the Purcell factor computed from different modal expansions for three QNMs treated as if they were orthogonal.The nanoparticle is a cylinder of gold in glass, of length 200 nm and diameter 30 nm, excited by an electric point dipole placed on axis, 10 nm above the rod, and oriented along its axis [11], see the inset in panel (a).Panel (a) shows the Purcell factor where three QNMs are excited.Panel (b) is a magnification of the large wavelength tail.Symbols represent exact FEM calculation.Equation (7) has been used to compute the Purcell factor, using either the scattered field expansion -thin dashed lineseither the total field expansion -thick solid lines.The new formulations Eq. ( 22), for scattered field and Eq. ( 23), for total field, are in red, and give better results than the former one [14], in blue.As both the first and second order model give practically the same results, only the first order has been represented.

Expansion coefficients
Numerical simulation tells us that two (or more) QNMs can sometimes be modeled independently -see Figs. 2 and 3 or [19].In this case, the relevant expansion coefficients are computed as if each QNMs were alone, for example using Eq. ( 21) and then the scattered field is obtained by summation on the set of QNM, using Eq. ( 3).
However, one may question this method.Indeed, as shown earlier [13], the excitation coefficient of the j 0 th QNM involves all the other QNMs as well.This effect is due to the fact that the QNMs are non-orthogonal [13], so that their excitation coefficients at a real frequency satisfy a matricial equation that reads: where M is a non-diagonal matrix.There, 2 ( ) m β ω is the vector of all expansion coefficients on m E  .Note the subscript '2' corresponds to the development of second order scattering problem, Eq. ( 2), and the subscript '1' to the first order problem Eq. ( 1).
An explicit expression of M can obtained without the need of PML integration: (i) using the eigen energy correction described in section 3 -for the diagonal terms.(ii) using an extra relation on mode orthogonality, for the off-diagonal terms.Once again, the latter is obtained by integrating Eq. ( 11) on 3   and then splitting the integral in two domains, 3 \ V  and V , to eventually obtain: Which is valid when m and p are different modes.Equation ( 25) complemented by Eq. ( 19) are an important result, that permit to expresses the crossover integral over 3   of two QNMs as an integral on a finite volume.Finally, Eq. ( 19) permits to obtain the diagonal terms: whereas the off-diagonal are, from Eq. ( 25): .
And for the first order scattering problem, we define in an analoguous way: Using mode orthogonality at frequencies m ω  and p ω  -see Eq. (S4) of [13] -, one gets: And then, the diagonal elements are: whereas the off-diagonal are: ( ) The matrices M gather the whole interaction of all QNMs, including the eigen-energy correcting term previously derived for a single QNM.Note that our definition of M do not use the conceptually complex PML integration [13], which appears to be unnecessary, as all the integrals on infinite domains can be rewritten as integrals on the finite volume of the nanoparticle(s).In the limit cancels out (both for first and second order description), however, for a real frequency there are in general non-zeros off-diagonal terms in M.
These expressions are valid for any number of coupled QNMs, and generalizes previous works that considered only one mode [13,14], or several orthogonal modes [13,19].Let us also mention that these expressions have been worked out for constant background permittivity, but the same expressions are also valid when the background is nonhomogeneous and non-dispersive.With little modifications, it can be generalized to the case of a non-homogeneous and dispersive medium, under an hypothesis of separation of the variable of ( , ) b ε ω r .
By looking directly at the formula, one understands that the strength of QNM coupling (the off-diagonal terms in M) depends on an analytic function that involves frequencies and permitivity, and on the crossover integral of the modal electric fields in the nanoparticule(s) volume.In particular, for particules made of an homogeneous dispersive material, when 0 , for m p ≠ , both modes are orthogonal (see Table 2).

Numerical simulations
In order to find a situation where non-orthogonality affects the modal expansion, one considers a dimer of coaxial gold nano cylinders coupled by a nanogap [13].Figure 4 shows how the modes of the system behaves in the complex plane, as the gap size is varied.One observes an anticrossing, sign of a strong interaction between the QNMs.Around the anticrossing (gap distance of 45 nm [13]), the crossover integral of the modes is of the order of the eigen-energy integral of each mode (see Table 1), and one can expect to see some effect of the non-orthogonality.and Non-orthogonal QNMs with order 2 expansion (black dashed).Note that the real part of QNM 2 and imaginary part of QNM 1, which are not represented, are predicted to be (almost) the same for the 3 models.Data on the two QNMs are summarized in Table 1.
Table 2. Coupling strongly two QNMs by a nanogap.To do so, we compute the excitation coefficient from Eq. ( 24) and Eq. ( 28).The results, see Fig. 4(b), show that the coefficients ( ) β ω are now clearly different whether one uses the first or second order Maxwell equation, or the simplified model [14], to obtain them.Going further, and computing the extinction spectrum of the dimer, see Fig. 5, one eventually observes a better prediction when non-orthogonality is taken into account and the second order model is used, Eq. ( 24).In this case, the disagreement between the exact result (symbols) and the model prediction is reduced by a factor ~2 at resonances, compared to predictions of the original model [14].Fig. 5. Extinction cross section for a gold dimer of geometry D 1 = 20 nm, L 1 = 80 nm ; D 2 = 85 nm, L 2 = 145 nm with a gap g = 45 nm, in a glass background.A plane wave polarized parallel to the dimer axis and propagating perpendicular to it is exciting the device.We show the predictions of different models: model with orthogonal QNMs [14] (red solid), Non-orthogonal QNMs with order 1 expansion (blue solid) and Non-orthogonal QNMs with order 2 expansion (black solid).The black dashed curve is obtained using Eq. ( 32), that adds a constant term to the second order expansion to model the scattered field (see text).Results from exact FEM calculation are represented with symbols.
Still, the agreement between models and exact result is not very satisfying -compare solid curves and symbols -, what was not the case previously (see Figs. 1-3).Additionally, the fact that first and second order model do not give the same results -compare blue and black solid lines -seems paradoxal.Indeed, the first and second order Maxwell equations should be equivalent for non-magnetic materials.
To solve this paradox, and improve the results, one has to recall that QNM expansion only models the resonant fields.A non-resonant, holomorphic function of frequency, should be added to get a more accurate description, see Eq. (3).Besides, whereas we considered two resonant poles in the wavelength range of interest, all the other poles would also give a contribution, provided their mode can be excited by the plane wave we considered.These higher order poles have their resonances at shorter wavelength, and their contribution could be evaluated accurately using e.g.Equation (24), if one knew all of them -what would necessitate another eigenmode solver than the iterative algorithm we use [14].
Eventually, to improve the model without adding much complexity, we assume that Then, by making a single additional computation, one computes ( ) cst f r as the difference between the exact scattered field at a given wavelength (we choose the local minimum of extinction, =0.9296 λ µm, in between the two extinction peaks), and the QNM expansion at this wavelength, using the second order expansion Eq. (24).Afterwards, ( ) cst f r is used in Eq. ( 32) to obtain a better estimation of the scattered field at all the wavelength.Eventually, one obtains a better description of the extinction cross section for the dimer (see the black dashed curve in Fig. 4).The discrepancy one can observe, going away from the resonances is attributed to our assumption that ( ) cst f r do not depend on frequency.Still, this simple method permits to improve the predictions of the models around the resonances, without much computational cost.And if one wants to keep a pure modal expansion (with f = 0), the second order model Eq. ( 24) gives the best results.

Conclusion and perspectives
We have shown that QNM theory can be improved by removing an hypothesis of slowly varying permittivity in the frequency domain.This leads to a more accurate expansion without introducing much complexity.Besides, we are able to define properly and compute simply the crossover integrals of QNM fields as integrals over the finite volume of the nanoparticule(s).
The field-based normalization that stems from this theory permits to use 2D axisymmetric FEM models, and to obtain the normalized QNMs in a much faster way than previous general tools for 3D systems [14].We believe this is a good starting point to check and further improve the QNM theory.To test our finite volume normalization, we computed the mode volume (without PML integration), and we obtained results is in good agreement with previous works (with PML integration or with volume and surface integration).
In the case several modes are considered, we present a rigorous theory that encompasses the mode non-orthogonality, avoiding the conceptually complex PML integration.These expressions are also valid for a NP embedded in a non-dispersive layered media, and can be extended to the case of a nanoparticule surrounded by a non-homogeneous and dispersive medium (e.g. a nanoparticule on a metallic substrate) with little work.
In the future, the analytic formula we obtained for the expansion on several modes could be used to check whether or not the QNM theory is complete, provided one is able to find all the modes.But to our opinion, the most important perspective in QNM theory now concerns the regularization of the Volume Integral Equation.Indeed, due to its strong singularity, the scattering operator do not satisfy the necessary hypothesis of compactness [25] to be expanded in a sum of modes, what could be a source of discrepancy for nanoparticles of larger sizes.
field emited by the dipole in the absence of NP, inside the homogeneous background of index b n .

Fig. 1 .
Fig. 1.Comparison of the different modal expansions for a single QNM.The extinction cross section -panel (a) -and Absorption cross section -panel (b) -, obtained from order 1 (green solid), order 2 (black dashed) and former model without self energy correction [14] (red solid) are compared to exact FEM simulations (open circles).The nanoparticle is a cylinder of gold in glass, of length 100 nm and diameter 30 nm, excited by a plane wave polarized along the rod axis and propagating perpendicular to the rod axis (see inset).Magnifications are shown in inset, where an horizontal dot-dashed line helps to distinguish between positive and negative cross sections.The difference between first and second order expansion is of the order of the accuracy of the numerical method.

Fig. 4 .
Fig.4.Panel (a): Anti-crossing of the QNM eigen wavelength in the complex plane, when the gap size (g ) of a doublet of gold nanocylinders is varied.The rod diameters and length are[13]: D 1 = 20 nm, L 1 = 80 nm ; D 2 = 85 nm, L 2 = 145 nm.The expansion coefficients, when a plane wave polarized parallel to the dimer axis and propagating perpendicular to it is exciting the device are represented on panel (b), when g = 45 nm.There, the real part of QNM 1 and imaginary part of QNM 2 are shown, from the predictions of different models: model with orthogonal QNMs[14] (red solid), Non-orthogonal QNMs with order 1 expansion (green solid) and Non-orthogonal QNMs with order 2 expansion (black dashed).Note that the real part of QNM 2 and imaginary part of QNM 1, which are not represented, are predicted to be (almost) the same for the 3 models.Data on the two QNMs are summarized in Table1.

Table 1 . Comparing the mode volumes.
3 r