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

Mathias Perrin

► To cite this version:

Mathias Perrin. Eigen-energy effects and non-orthogonality in the quasi-normal mode expansion of Maxwell equations. Optics Express, 2016, 24 (24), pp.27137 (1-15). 10.1364/OE.24.027137. hal-01426676

HAL Id: hal-01426676

https://hal.science/hal-01426676

Submitted on 4 Jan 2017

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution - ShareAlike 4.0 International License
# Eigen-energy effects and non-orthogonality in the quasi-normal mode expansion of Maxwell equations

MATHIAS PERRIN

Laboratoire Ondes et Matière d'Aquitaine LOMA, Université Bordeaux, CNRS, UMR 5798, 33405 Talence - France

*mathias.perrin@u-bordeaux.fr

**Abstract:** 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. Eventually, 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: first-order 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).
17. B. Vial, F. Zolla, A. Nicolet, and M. Commandré, "Quasimodal expansion of electromagnetic fields in open two dimensional structures," Phys. Rev. A **89**(2), 023829 (2014).

18. B. Vial, G. Demesy, F. Zolla, A. Nicolet, M. Commandré, C. Hecquet, T. Begou, S. Tisserand, S. Gautier, and V. Sauget, "Resonant metamaterial absorbers for infrared spectral filtering: quasimodal analysis, design, fabrication, and characterization," J. Opt. Soc. Am. B **31**(6), 1339-1346 (2014).

19. J. Rosenkrantz de Lasson, P. T. Kristensen, J. Mørk, and N. Gregersen, "Semianalytical quasi-normal mode theory for the local density of states in coupled photonic crystal cavity-waveguide structures," Opt. Lett. **40**(24), 5790-5793 (2015).

20. F. Yang, H. Liu, H. Jia, and Y. Zhong, "Analytical description of quasi-normal mode in resonant plasmonic nano-cavities," J. Opt. **18**(3), 035003 (2016).

21. M. B. Doost, "Resonant-state-expansion Born approximation for waveguides with dispersion," Phys. Rev. A **93**(2), 023835 (2016).

22. M. B. Doost, W. Langbein, and E. A. Muljarov, "Resonant state expansion applied to three-dimensional open optical systems," Phys. Rev. A **90**(1), 013834 (2014).

23. G. W. Hanson and A. B. Yakovlev, *Operator Theory for Electromagnetics* (Springer, 2002).

24. J. Mäkitalo, M. Kauranen, and S. Suuriniemi, "Modes and resonances of plasmonic scatterers," Phys. Rev. B **89**(16), 165429 (2014).

25. R. Mennicken and M. Möller, *Non-Self-Adjoint Boundary Eigenvalue Problems*, Jan van Mill ed. (Elsevier, 2003).

26. R.-C. Ge and S. Hughes, "Quasinormal mode theory and modelling of electron energy loss spectroscopy for plasmonic nanostructures," J. Opt. **18**(5), 054002 (2016).

27. R.-C. Ge and S. Hughes, "Quantum dynamics of two quantum dots coupled through localized plasmons: An intuitive and accurate quantum optics approach using quasinormal modes," Phys. Rev. B **92**(20), 205420 (2015).

28. O. Lozan, M. Perrin, B. Ea-Kim, J.-M. Rampnoux, S. Dilhaire, and P. Lalanne, "Anomalous light absorption around subwavelength apertures in metal films," Phys. Rev. Lett. **112**(19), 193903 (2014).

29. P. Lalanne, J. P. Hugonin, H. T. Liu, and B. Wang, "A microscopic view of the electromagnetic properties of sub-λ metallic surfaces," Surf. Sci. Rep. **64**(10), 453-469 (2009).

30. A. E. Siegman, *Lasers* (University Science Books, 1986).

31. J. Waxenegger, A. Trügler, and U. Hohenester, "Plasmonics simulations with the MNPBEM toolbox: consideration of substrates and layer structures," Comput. Phys. Commun. **193**, 138-150 (2015).

32. R.-C. Ge, P. T. Kristensen, J. F. Young, and S. Hughes, "Quasinormal mode approach to modeling light-emission and propagation in nanoplasmonics," New J. Phys. **16**(11), 113048 (2014).

33. L. Novotny and B. Hecht, *Principles of Nano-Optics* (Cambridge 2006).

34. X. Zambrana-Puyalto and N. Bonod, "Purcell factor of spherical Mie resonators," Phys. Rev. B **91**(19), 195422 (2015).

35. R. Loudon, "Propagation of electromagnetic energy through an absorbing dielectric," J. Phys. A **3**(3), 233-245 (1970).

36. E. A. Muljarov and W. Langbein, "Resonant-state expansion of dispersive open optical systems: creating gold from sand," Phys. Rev. B **93**(7), 075417 (2016).

# 1. 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-6], sensors [7-9], devices for frequency conversion [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-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-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.

## 2. Formulation of the scattering problem

In the following, we consider monochromatic fieds, and use the $exp(+i\omega t)$ notation.

Two ways are commonly used to set the modal problem. This can be through a first order representation [13,14]:

$$ \nabla \wedge \tilde{\mathbf{H}}_m = i\tilde{\omega}_m \varepsilon(\tilde{\omega}_m, \mathbf{r})\tilde{\mathbf{E}}_m, \quad \text{and} \quad \nabla \wedge \tilde{\mathbf{E}}_m = -i\tilde{\omega}_m \mu(\tilde{\omega}_m, \mathbf{r})\tilde{\mathbf{H}}_m, \qquad (1) $$

or a second order representation [16,21,22]:

$$ [\nabla \wedge \nabla \wedge] \mathbf{E}_m(\mathbf{r}) = (\tilde{\omega}_m)^2 \mu_0 \varepsilon(\tilde{\omega}_m, \mathbf{r}) \mathbf{E}_m(\mathbf{r}), \qquad (2) $$

for non-magnetic materials. In Eqs. (1)-(2), $\varepsilon(\omega, \mathbf{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 $\varepsilon_b = n_b^2$ – 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 $\mathbf{E}_b(\omega, \mathbf{r})$, and we write the scattered field at frequency $\omega$, $\mathbf{E}_s(\omega, \mathbf{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]:

$$ \mathbf{E}_s(\omega, \mathbf{r}) = \sum_m \beta_m(\omega) \tilde{\mathbf{E}}_m(\mathbf{r}) + \mathbf{f}(\omega, \mathbf{r}). \qquad (3) $$

Note that, the QNM field of eigenvector $\tilde{k} = \tilde{\omega}n_b/c$ diverges as $exp(-i\tilde{k}r)/r$ 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 $\text{Im}(\tilde{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

$$ [\nabla \wedge \nabla \wedge] \mathbf{E}_s(\omega, \mathbf{r}) = (\omega)^2 \mu_0 \varepsilon(\omega, \mathbf{r}) \mathbf{E}_v(\omega, \mathbf{r}) - i\omega\mu_0 \mathbf{J}(\omega, \mathbf{r}), \qquad (4) $$

where $\mathbf{J}(\omega, \mathbf{r}) = i\omega\Delta\varepsilon(\omega, \mathbf{r})\mathbf{E}_b(\omega, \mathbf{r})$ is the current source for the scattered field [14], and $\Delta\varepsilon$ represents the permittivity change that models the nanoparticule(s).

Our goal is to find the expansion coefficients $\beta_m(\omega)$ that will give the best approximation to the scattered field, so that we will consider $\mathbf{f}(\omega, \mathbf{r}) = 0$. Then, once the $\beta_m(\omega)$ are obtained, the absorption and extinction cross sections can be respectively computed as [14]:

$$ \sigma_A = -\frac{\omega}{2S_0} \iiint_V \text{Im}[\varepsilon(\omega, \mathbf{r})] |\mathbf{E}_s(\omega, \mathbf{r}) + \mathbf{E}_b(\omega, \mathbf{r})|^2 d\mathbf{r}. \qquad (5) $$

and
$$ \sigma_E = -\frac{\omega}{2S_0} \iiint_V \mathrm{Im}(\Delta\varepsilon(\omega, \mathbf{r})(\mathbf{E}_b(\omega, \mathbf{r}) + \mathbf{E}_s(\omega, \mathbf{r})) \cdot \mathbf{E}_b^*(\omega, \mathbf{r})) \mathrm{d}\mathbf{r}. \quad (6) $$

Additionally, the Purcell factor of an electric dipole $\mathbf{J}_0$ located at $\mathbf{r}_0$ is obtained [13,14] as:

$$ P(\omega) = -\frac{6\pi\varepsilon_0 c^3}{\omega^2 n_b \|\mathbf{J}_0\|^2} \mathrm{Re}[\mathbf{J}_0 \cdot \mathbf{E}(\omega, \mathbf{r}_0)], \quad (7) $$

where $\mathbf{E}(\omega, \mathbf{r}) = \mathbf{E}_s(\omega, \mathbf{r}) + \mathbf{E}_{\mathrm{bDIPOLE}}(\omega, \mathbf{r})$ is the total field, and $\mathbf{E}_{\mathrm{bDIPOLE}}(\omega, \mathbf{r})$ represents the field emitted by the dipole in the absence of NP, inside the homogeneous background of index $n_b$.

# 3. Excitation of a single QNM

## 3.1 Normalization of the QNM

### 3.1.1 Field normalisation

In this section, we assume only one QNM is excited.

The QNM normalization can be 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 $\int_{\mathbb{R}^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]:

$$ \int_{\mathbb{R}^3} \tilde{\mathbf{E}}(\mathbf{r}) \frac{\partial[\omega\varepsilon]}{\partial\omega} \Big|_{(\tilde{\omega}, \mathbf{r})} \tilde{\mathbf{E}}(\mathbf{r}) - \tilde{\mathbf{H}}(\mathbf{r}) \frac{\partial[\omega\mu]}{\partial\omega} \Big|_{(\tilde{\omega}, \mathbf{r})} \tilde{\mathbf{H}}(\mathbf{r}) \mathrm{d}\mathbf{r} = 1, \quad (8) $$

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 $\tilde{\mathbf{H}}\tilde{\mathbf{H}} = -[\nabla \wedge \tilde{\mathbf{E}}] \cdot [\nabla \wedge \tilde{\mathbf{E}}] / (\mu\omega)^2$. Then, using the vectorial identity $\mathbf{B} \cdot [\nabla \wedge \mathbf{A}] - \mathbf{A} \cdot [\nabla \wedge \mathbf{B}] = \nabla(\mathbf{A} \wedge \mathbf{B})$, one obtains an equivalent normalization condition:

$$ \int_{\mathbb{R}^3} \tilde{\mathbf{E}}(\mathbf{r}) \frac{\partial[\omega^2 \varepsilon]}{\partial\omega} \Big|_{(\tilde{\omega}, \mathbf{r})} \tilde{\mathbf{E}}(\mathbf{r}) \mathrm{d}\mathbf{r} = \tilde{\omega} \quad (9) $$

Consequently, the QNM norm, as defined by [16] now reads:

$$ \langle \tilde{\mathbf{E}} | \tilde{\mathbf{E}} \rangle = \frac{1}{2\tilde{\omega}} \int_{\mathbb{R}^3} \tilde{\mathbf{E}}(\mathbf{r}) \frac{\partial[\omega^2 \varepsilon]}{\partial\omega} \Big|_{(\tilde{\omega}, \mathbf{r})} \tilde{\mathbf{E}}(\mathbf{r}) \mathrm{d}\mathbf{r} = \frac{1}{2} \quad (10) $$

In order to define $\tilde{\mathbf{E}}(\mathbf{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]:

$$ \mathbf{E}_1 \cdot (\nabla \wedge \nabla \wedge \mathbf{E}_2) - \mathbf{E}_2 \cdot (\nabla \wedge \nabla \wedge \mathbf{E}_1) = \nabla \cdot (\mathbf{E}_2 \wedge \nabla \wedge \mathbf{E}_1 - \mathbf{E}_1 \wedge \nabla \wedge \mathbf{E}_2), \quad (11) $$
Then, taking $E_1 = E_s(\omega, \mathbf{r})$ and $E_2 = \tilde{E}(\mathbf{r})$, and integrating on the whole space (including PML), one eventually gets, using Eq. (4):

$$ \beta(\omega) \int_{\mathbb{R}^3} \tilde{\mathbf{E}}(\mathbf{r}) [\tilde{\omega}^2 \varepsilon(\tilde{\omega}, \mathbf{r}) - \omega^2 \varepsilon(\omega, \mathbf{r})] \tilde{\mathbf{E}}(\mathbf{r}) \, d\mathbf{r} = \omega^2 \int_{\mathbb{R}^3} \Delta\varepsilon(\omega, \mathbf{r}) E_b(\omega, \mathbf{r}) \cdot \tilde{\mathbf{E}}(\mathbf{r}) \, d\mathbf{r}, \quad (12) $$

Lets now consider a frequency $\Omega$ close to the pole, $\tilde{\omega}$, so that a Taylor expansion of the permittivity is valid,

$$ \Omega^2 \varepsilon(\Omega, \mathbf{r}) = \tilde{\omega}^2 \varepsilon(\tilde{\omega}, \mathbf{r}) + (\Omega - \tilde{\omega}) \frac{\partial[\omega^2 \varepsilon]}{\partial \omega} \Bigg|_{(\tilde{\omega}, \mathbf{r})}, \quad (13) $$

and Eq. (12) becomes:

$$ \beta(\Omega) = \frac{\Omega^2 \int_{\mathbb{R}^3} \Delta\varepsilon(\Omega, \mathbf{r}) E_b(\Omega, \mathbf{r}) \cdot \tilde{\mathbf{E}}(\mathbf{r}) \, d\mathbf{r}}{\Omega - \tilde{\omega}}. \quad (14) $$

When $\Omega \to \tilde{\omega}$, the scattered field increases dramatically [14] and one can neglect all the contributions except the $\tilde{\omega}$ pole in the QNM expansion, to write $E_s(\Omega, \mathbf{r}) = \beta(\Omega) \tilde{\mathbf{E}}(\mathbf{r})$. Multiplying both sides of the latter equation by $\Delta\varepsilon(\Omega, \mathbf{r}) E_b(\Omega, \mathbf{r})$ and integrating, one gets:

$$ \int_V \Delta\varepsilon(\Omega, \mathbf{r}) E_b(\Omega, \mathbf{r}) \cdot E_s(\Omega, \mathbf{r}) \, d\mathbf{r} = \beta(\Omega) \int_V \Delta\varepsilon(\Omega, \mathbf{r}) E_b(\Omega, \mathbf{r}) \cdot \tilde{\mathbf{E}}(\mathbf{r}) \, d\mathbf{r}. \quad (15) $$

Inserting Eq. (14) in Eq. (15) and using that $\Omega \approx \tilde{\omega}$, one gets:

$$ [\beta(\Omega)]^2 = -\Omega \frac{\int_V \Delta\varepsilon(\Omega, \mathbf{r}) E_b(\Omega, \mathbf{r}) \cdot E_s(\Omega, \mathbf{r}) \, d\mathbf{r}}{\Omega - \tilde{\omega}}, \quad (16) $$

and eventually one obtains the QNM definition:

$$ \tilde{\mathbf{E}}(\mathbf{r}) = i \sqrt{\frac{\Omega - \tilde{\omega}}{\Omega}} \frac{\mathbf{E}_s(\Omega, \mathbf{r})}{\sqrt{\int_V \Delta\varepsilon(\Omega, \mathbf{r}) \mathbf{E}_b(\Omega, \mathbf{r}) \cdot \mathbf{E}_s(\Omega, \mathbf{r}) d\mathbf{r}}}, \quad (17) $$

where $\Omega$ 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 $\Omega$ converges towards $\tilde{\omega}$. Besides, using Eq. (17) to define $\tilde{\mathbf{E}}(\mathbf{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 qth azimuthal modes, whose eigen-fields vary in $\exp(iq\phi)$ 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):

$$V_{QNM} = \frac{1}{2\epsilon_0 n_b^2 [\tilde{\mathbf{E}}(\mathbf{r}_0) \cdot \mathbf{u}]^2}, \quad (18)$$

Where $\mathbf{r}_0$ is a specific point around the NP, and $\mathbf{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 meaningful comparison, all the data have been converted if necessary into the $\exp(i\omega t)$ convention we use for harmonic fields.

<table><caption>**Table 1. Comparing the mode volumes.**</caption><thead><tr><th>Nanosystem</th><th>λ̃ (nm) in exp(+iωt) convention</th><th>V<sub>QNM</sub> in exp(+iωt) convention</th></tr></thead><tbody><tr><td>nanorod <sup>A</sup> (this work)</td><td>918 - 47i</td><td>(5.3 + 0.32i) λ<sup>3</sup> / 10<sup>4</sup></td></tr><tr><td>nanorod <sup>A</sup> ([13])</td><td>920 - 47i</td><td>(5 + 0.4 i) λ<sup>3</sup> / 10<sup>4</sup></td></tr><tr><td rowspan="2">Dimer <sup>B</sup> (this work)</td><td>933 - 78i</td><td>(-3.3 + 7.0 i) λ<sup>3</sup> / 10<sup>4</sup></td></tr><tr><td>952 - 108i</td><td>(4.1-3.1 i) λ<sup>3</sup> / 10<sup>4</sup></td></tr><tr><td rowspan="2">Dimer <sup>B</sup> ([13])</td><td>933 - 77i</td><td>(-3 + 7 i) λ<sup>3</sup> / 10<sup>4</sup></td></tr><tr><td>958 - 108i</td><td>(4 - 3 i) λ<sup>3</sup> / 10<sup>4</sup></td></tr><tr><td>Dimer <sup>C</sup> (this work)</td><td>1023 - 70i</td><td>(35.6 + 0.52 i) r<sup>3</sup></td></tr><tr><td>Dimer <sup>C</sup> ([16])</td><td>1025 - 72i</td><td>(35.19 + 0.96 i) r<sup>3</sup></td></tr></tbody></table>

**Nanorod A:** 100nm length and 15nm radius gold nanorod in glass ($n_b = 1.5$). $\mathbf{r}_0$ is placed on axis, 10 nm above the rod, $\mathbf{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 $\lambda^3 / 10^4$ where $\lambda$ 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.

### 3.1.2 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 $\int_{\mathbb{R}^3} \tilde{\mathbf{E}}(\mathbf{r})\epsilon(\mathbf{r})\tilde{\mathbf{E}}(\mathbf{r})$ 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:
$$2 \int_{\mathbb{R}^3 \setminus V} \tilde{\mathbf{E}}(\mathbf{r}) \varepsilon_b(\mathbf{r}) \tilde{\mathbf{E}}(\mathbf{r}) \, \mathrm{d}\mathbf{r} = 1 - \frac{1}{\tilde{\omega}} \int_V \tilde{\mathbf{E}}(\mathbf{r}) \frac{\partial[\omega^2 \varepsilon]}{\partial \omega} \Bigg|_{(\tilde{\omega}, \mathbf{r})} \tilde{\mathbf{E}}(\mathbf{r}) \, \mathrm{d}\mathbf{r}. \quad (19)$$

where the symbol $\mathbb{R}^3 \setminus 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: $\varepsilon_b(\mathbf{r})$. The very important point in Eq. (19) is that it permits to connect an integral on $\mathbb{R}^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).

## 3.2 Derivation of expansion coefficients

Equation (19) indeed permits to express rigourously the $\mathbb{R}^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 coefficient for the second order problem is:

$$\beta_{\text{order}2}(\omega) = \frac{-\omega^2 \int_V \Delta\varepsilon(\omega, \mathbf{r}) \mathbf{E}_b(\omega, \mathbf{r}) \cdot \tilde{\mathbf{E}}(\mathbf{r}) \, \mathrm{d}\mathbf{r}}{(\omega^2 - \tilde{\omega}^2) \left( \frac{1}{2} + \int_V \tilde{\mathbf{E}}(\mathbf{r}) \left[ \frac{\omega^2 \varepsilon(\omega, \mathbf{r}) - \tilde{\omega}^2 \varepsilon(\tilde{\omega}, \mathbf{r})}{\omega^2 - \tilde{\omega}^2} - \frac{1}{2\tilde{\omega}} \frac{\partial[\omega^2 \varepsilon]}{\partial \omega} \Bigg|_{(\tilde{\omega}, \mathbf{r})} \right] \tilde{\mathbf{E}}(\mathbf{r}) \, \mathrm{d}\mathbf{r} \right)}, \quad (21)$$

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 $\varepsilon_b(\omega, \mathbf{r}) = g(\mathbf{r})h(\omega)$. Then, a term depending on $\frac{dh}{d\omega}$ 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:

$$\beta_{\text{order}1}(\omega) = -\omega \frac{\int_V \Delta\varepsilon(\omega, \mathbf{r}) \mathbf{E}_b(\omega, \mathbf{r}) \cdot \tilde{\mathbf{E}}(\mathbf{r}) \, \mathrm{d}\mathbf{r}}{(\omega - \tilde{\omega}) \left( 1 + \int_V \tilde{\mathbf{E}}(\mathbf{r}) \left[ \frac{\omega\varepsilon(\omega, \mathbf{r}) - \tilde{\omega}\varepsilon(\tilde{\omega}, \mathbf{r})}{\omega - \tilde{\omega}} - \frac{\partial[\omega\varepsilon]}{\partial\omega} \Bigg|_{(\tilde{\omega}, \mathbf{r})} \right] \tilde{\mathbf{E}}(\mathbf{r}) \, \mathrm{d}\mathbf{r} \right)}, \quad (22)$$

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 $\int_V \tilde{\mathbf{E}} \left[ \frac{\partial(\omega\varepsilon)}{\partial\omega} \Bigg|_{\omega} \right] \tilde{\mathbf{E}}$ 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 $\omega$ is close enough to $\tilde{\omega}$, both the first and second order expression reduce to the previously derived relations for the analytic expansion coefficient [14].
### 3.3 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 $ε_b = 2.25$. We took a Drude model for gold, assuming $ε(ω) = 1 - ϵ_p^2 / (ω^2 - iωγ)$, with $ω_p = 1.26 × 10^{16} s^{-1}$ and $γ = 1.41 × 10^{14} s^{-1}$ [14].

#### 3.3.1 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 $β(ω)$.

<figure><img src="image_2.png" /><figcaption>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.</figcaption></figure>

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 $λ_1 = 1.5 - 0.098i$ µm, $λ_2 = 0.8 - 0.023i$ µm and $λ_3 = 0.6 - 0.013i$ µm. 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 anti-symmetric 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.

<figure><img src="image_3.png" /><figcaption>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.</figcaption></figure>

### 3.3.2 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, $\int_V \Delta\epsilon(\omega, \mathbf{r}) E_{\text{bDIPOLE}}(\omega, \mathbf{r}) \tilde{E}(\mathbf{r}) d\mathbf{r}$, at each frequency, whereas the total field expansion simply necessitates to evaluate the scalar product between the dipole current $J_0$ and the QNM field at the dipole location, $\tilde{E}(r_0)$, which has a constant value. This later is therefore more convenient to use, even though 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:

$$ \alpha_{\text{order1}}(\omega) = \frac{i \mathbf{J}_0 \cdot \tilde{\mathbf{E}}(\mathbf{r}_0)}{(\omega - \tilde{\omega}) \left( 1 + \int_V \tilde{\mathbf{E}}(\mathbf{r}) \left[ \frac{\omega \epsilon(\omega, \mathbf{r}) - \tilde{\omega} \epsilon(\tilde{\omega}, \mathbf{r})}{\omega - \tilde{\omega}} - \frac{\tilde{\partial}[\omega \epsilon]}{\tilde{\partial}\omega} \right]_{(\tilde{\omega}, \mathbf{r})} \tilde{\mathbf{E}}(\mathbf{r}) d\mathbf{r} \right)}, \quad (23) $$
where the total field is well approximated by $E(\omega, \mathbf{r}) = \alpha(\omega)\tilde{E}(\mathbf{r})$.

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).

<figure><img src="image_4.png" /><figcaption>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 lines – either 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.</figcaption></figure>

# 4. Generalization to several QNMs

## 4.1 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:

$$[M_2^{pm}][\beta_2^m(\omega)] = -\omega^2 [\int \Delta\epsilon(\omega)E_b \cdot \tilde{E}_p], \quad (24)$$
where M is a non-diagonal matrix.

There, $\beta_m^m(\omega)$ is the vector of all expansion coefficients on $\tilde{E}_m$. 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 $\mathbb{R}$ and then splitting the integral in two domains, $\mathbb{R} \setminus V$ and $V$, to eventually obtain:

$$ \int_{\mathbb{R}^3 \setminus V} \varepsilon_b(\mathbf{r}) \tilde{E}_m \cdot \tilde{E}_p = - \int_V \tilde{E}_m \left[ \frac{\tilde{\omega}_m^2 \varepsilon(\tilde{\omega}_m) - \tilde{\omega}_p^2 \varepsilon(\tilde{\omega}_p)}{\tilde{\omega}_m^2 - \tilde{\omega}_p^2} \right] \tilde{E}_p, \quad (25) $$

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 $\mathbb{R}^3$ of two QNMs as an integral on a finite volume.

Finally, Eq. (19) permits to obtain the diagonal terms:

$$ M_2^{pp} = (\omega^2 - \tilde{\omega}_p^2) \left[ \frac{1}{2} + \int_V \tilde{E}_p \left[ \frac{\omega^2 \varepsilon(\omega) - \tilde{\omega}_p^2 \varepsilon(\tilde{\omega}_p)}{\omega^2 - \tilde{\omega}_p^2} - \frac{1}{2\tilde{\omega}} \frac{\partial(\omega^2 \varepsilon)}{\partial \omega} \bigg|_{\tilde{\omega}_p} \right] \tilde{E}_p \right], \quad (26) $$

whereas the off-diagonal are, from Eq. (25):

$$ M_2^{pm} = (\omega^2 - \tilde{\omega}_p^2) \int_V \tilde{E}_m \left[ \frac{\omega^2 \varepsilon(\omega) - \tilde{\omega}_p^2 \varepsilon(\tilde{\omega}_p)}{\omega^2 - \tilde{\omega}_p^2} - \frac{\tilde{\omega}_m^2 \varepsilon(\tilde{\omega}_m) - \tilde{\omega}_p^2 \varepsilon(\tilde{\omega}_p)}{\tilde{\omega}_m^2 - \tilde{\omega}_p^2} \right] \tilde{E}_p. \quad (27) $$

And for the first order scattering problem, we define in an analogous way:

$$ [M_1^{pm}] [\beta_1^m(\omega)] = -\omega \left[ \int \Delta\varepsilon(\omega) \mathbf{E}_b \cdot \tilde{\mathbf{E}}_p \right] \quad (28) $$

Using mode orthogonality at frequencies $\tilde{\omega}_m$ and $\tilde{\omega}_p$ – see Eq. (S4) of [13] –, one gets:

$$ \int_{\mathbb{R}^3 \setminus V} \tilde{E}_m \varepsilon_b \tilde{E}_p - \tilde{H}_m \mu \tilde{H}_p = - \int_V \tilde{E}_m \left[ \frac{\tilde{\omega}_m \varepsilon(\tilde{\omega}_m) - \tilde{\omega}_p \varepsilon(\tilde{\omega}_p)}{\tilde{\omega}_m - \tilde{\omega}_p} \right] \tilde{E}_p - \mu \tilde{H}_m \tilde{H}_p \quad (29) $$

And then, the diagonal elements are:

$$ M_1^{pp} = (\omega - \tilde{\omega}_p) \left[ 1 + \int_V \tilde{E}_p \left[ \frac{\omega \varepsilon(\omega) - \tilde{\omega}_p \varepsilon(\tilde{\omega}_p)}{\omega - \tilde{\omega}_p} - \frac{\partial(\omega \varepsilon)}{\partial \omega} \bigg|_{\tilde{\omega}_p} \right] \tilde{E}_p \right], \quad (30) $$

whereas the off-diagonal are:

$$ M_1^{pm} = (\omega - \tilde{\omega}_p) \int_V \tilde{E}_m \left[ \frac{\omega \varepsilon(\omega) - \tilde{\omega}_p \varepsilon(\tilde{\omega}_p)}{\omega - \tilde{\omega}_p} - \frac{\tilde{\omega}_m \varepsilon(\tilde{\omega}_m) - \tilde{\omega}_p \varepsilon(\tilde{\omega}_p)}{\tilde{\omega}_m - \tilde{\omega}_p} \right] \tilde{E}_p. \quad (31) $$

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 $\omega \to \tilde{\omega}_m$ or $\omega \to \tilde{\omega}_p$, the off-
diagonal term $M^{pm}$ 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 non-
homogeneous 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 $\epsilon_b(\omega, \mathbf{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 permittivity, 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 $\int_V \tilde{E}_m \tilde{E}_p = 0$, for $m \neq p$, both modes are orthogonal (see Table 2).

## 4.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.

<figure><img src="image_5.png" /><figcaption>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₁ = 20 nm, L₁ = 80 nm; D₂ = 85 nm, L₂ = 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 Table 1.</figcaption></figure>
<table><caption>Table 2. Coupling strongly two QNMs by a nanogap.</caption><thead><tr><th>g (nm)</th><th>δ̃<sub>1</sub> × 10<sup>-15</sup> (s<sup>-1</sup>)</th><th>δ̃<sub>2</sub> × 10<sup>-15</sup> (s<sup>-1</sup>)</th><th>∫<sub>V</sub> Ẽ<sub>1</sub> Ẽ<sub>1</sub> × 10<sup>-9</sup></th><th>∫<sub>V</sub> Ẽ<sub>1</sub> Ẽ<sub>2</sub> × 10<sup>-9</sup></th><th>∫<sub>V</sub> Ẽ<sub>2</sub> Ẽ<sub>2</sub> × 10<sup>-9</sup></th></tr></thead><tbody><tr><td>15</td><td>1.786 + 0.180i</td><td>2.127 + 0.197i</td><td>0.762 + 0.384i</td><td>0.436-0.118i</td><td>0.834-0.108i</td></tr><tr><td>45</td><td>1.953 + 0.220i</td><td>2.003 + 0.169i</td><td>0.136 + 1.428i</td><td>1.282 + 0.619i</td><td>1.481-1.114i</td></tr><tr><td>70</td><td>1.978 + 0.104i</td><td>1.982 + 0.292i</td><td>1.332-0.110i</td><td>0.195 + 0.365i</td><td>0.295 + 0.419i</td></tr></tbody></table>

The dimer size and geometry are described Fig. 4. g is the gap size

To do so, we compute the excitation coefficient from Eq. (24) and Eq. (28). The results, see Fig. 4(b), show that the coefficients $\beta(\omega)$ 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].

<figure><img src="image_6.png" /><figcaption>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.</figcaption></figure>

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 $\mathbf{f}(\omega, \mathbf{r})$ do not depend on $\omega$, and denote it $\mathbf{f}_{\omega}(\mathbf{r})$, what turns Eq. (3) into:

$$ \mathbf{E}_{\omega}(\omega, \mathbf{r}) = \sum_{m} \beta_{m}(\omega) \tilde{\mathbf{E}}_{m}(\mathbf{r}) + \mathbf{f}_{\omega}(\mathbf{r}). \qquad (32) $$

Then, by making a single additional computation, one computes $\mathbf{f}_{\omega}(\mathbf{r})$ as the difference between the exact scattered field at a given wavelength (we choose the local minimum of extinction, $\lambda=0.9296$ µm, in between the two extinction peaks), and the QNM expansion at this wavelength, using the second order expansion Eq. (24). Afterwards, $\mathbf{f}_{\omega}(\mathbf{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 $\mathbf{f}_{\omega}(\mathbf{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 $\mathbf{f}=0$), the second order model Eq. (24) gives the best results.

## 5. 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.

## Funding

French National Research Agency (ANR), ANR-13-BS08-0003 / PhotoLighT State of Aquitaine, SMOLED project (number 2014 - 1R6306).

## Acknowledgment

I acknowledge Bernadette Bergeret for the help provided in gathering some references, and Frédéric Gruy, for discussions on [25].