Gaussian regularization for resonant states: open and dispersive optical systems

TL;DR

This paper introduces Gaussian regularization techniques to analytically compute inner products of resonant states in open, dispersive optical systems, enabling accurate spectral analysis aligned with Mie theory results.

Contribution

It presents a novel Gaussian regularization approach for calculating resonant state inner products in open systems with dispersive media, bridging the gap with exact Mie theory.

Findings

Gaussian regularization yields analytical inner product results.

Method accurately reproduces Mie theory calculations.

Applicable to systems with electric and magnetic dispersion.

Abstract

Resonant States (RS), also known as Quasi-Normal Modes (QNMs), are eigenstates that arise in spectral expansions of linear response functions of open systems. Manipulation of these spatially `divergent' oscillating functions requires a departure from the usual definitions of inner product, normalization and orthogonality typically encountered in the studies of closed systems. We show that once RS fields are expanded on a multipole basis, Gaussian regularization methods provide \emph{analytical} results for crucial RS inner product integrals \added{in the problematic region exterior to the scattering system}. Our demonstrations are carried out in the context of light scattering by spatially bounded objects composed of both electrically and magnetically dispersive media, with demonstrative analytic calculations being shown to \emph{completely} retrieve the results of exact Mie theory.

Tables

Table 1. Table 1: RS eigenvalues, z α subscript 𝑧 𝛼 z_{\alpha} , and associated residues, ℛ α = i / 𝒩 α 2 subscript ℛ 𝛼 𝑖 subscript superscript 𝒩 2 𝛼 \mathcal{R}_{\alpha}=i/\mathcal{N}^{2}_{\alpha} , for a non-dispersive dielectric medium with dielectric contrast of ε = 16 𝜀 16 \varepsilon=16 and no magnetic contrast, μ = 1 𝜇 1 \mu=1 .

$\alpha (q,n,\left[\ell \right])$	$z_{\alpha }$	${ℛ}_{\alpha }$
$(e,1,\left[1\right])$	$1.0395-i0.500935$	$-0.236682+i0.231492$
$(e,1,\left[2\right])$	$1.05273-i0.0723549$	$0.0659905-i0.0579972$
$(e,1,\left[3\right])$	$1.92043-i0.082005$	$0.0748408-i0.0282738$
$(e,1,\left[4\right])$	$2.7227-i0.073007$	$0.00279437-i0.0683107$
$(e,2,\left[0\right])$	$-i1.6797303$	$i0.146892$
$(e,2,\left[1\right])$	$1.377484-i0.0118433$	$0.00184613-i0.0118059$
$(e,2,\left[2\right])$	$2.071446-i0.667649$	$-0.305381+i0.277002$
$(h,1,\left[0\right])$	$-i1.250038$	$i0.136765$
$(h,1,\left[1\right])$	$0.7537823-i0.0240302$	$-0.00601759-i0.0229898$
$(h,1,\left[2\right])$	$1.5414631-i0.0459254$	$-0.0394075-i0.0195948$
$(h,2,\left[1\right])$	$0.870513-i1.75259$	$-0.0521306+i0.140046$
$(h,2,\left[2\right])$	$1.0957165-i0.00684025$	$-0.000482964-i0.00681678$

Table 2. Table 2: Illustrative electric dipole RSs, α ( e , n = 1 , ℓ ) \alpha(e,n=1,\ell) , for plasmonic particles: for the dispersion models of Eq.​​   51b and Eq.​​   52 for gold and Eq.​​   51a for silver. The complex wavelengths and dielectric constants are given in addition to the Mie coefficient residue, ℛ n , α subscript ℛ 𝑛 𝛼 \mathcal{R}_{n,\alpha} , associated with RS normalization is given in the last column.

$R$ (nm)	$\ell$	Disp.	${\lambda }_{\alpha }$ (nm)	${\epsilon }_{\alpha }$	${ℛ}_{\alpha }$
100	1	(51b)	$606.976+i239.112$	$-13.7606-i12.5419$	$-0.217942+i0.034889$
100	1	(52)	$592.227+i210.097$	$-5.90132-i12.9326$	$-0.309295-i0.00587748$
80	1	(51b)	$505.163+i174.433$	$-9.48069-i7.58828$	$-0.2111-i0.00385708$
80	1	(52)	$547.815+i80.062$	$-3.64747-i3.52209$	$0.0167796-i0.10005$
100	1	(51a)	$600.211+i231.333$	$-11.1356-i14.0631$	$-0.236629+i0.0266621$
$\cdot$	2	$\cdot$	$290.678+i33.6325$	$0.704812-i0.966352$	$0.178962+i0.124692$
$\cdot$	3	$\cdot$	$217.845+i20.8224$	$2.58005-i0.449824$	$0.0139492-i0.235957$
80	1	$\cdot$	$500.306+i156.443$	$-6.78066-i7.89519$	$-0.235268-i0.0509908$
$\cdot$	2	$\cdot$	$281.965+i50.9047$	$1.03036-i1.44187$	$0.12746+i0.24896$
$\cdot$	3	$\cdot$	$192.72+i20.4675$	$3.11019-i0.394101$	$0.10372-i0.210644$

Equations

$$|\Psi ⟩\equiv |\begin{array}{c}\mathbf{E}\\ \mathbf{H}\end{array}⟩,|J⟩\equiv |\begin{array}{c}\mathbf{j}\\ 0\end{array}⟩,$$

$$\epsilon (\mathbf{r},\omega )\equiv {\epsilon }_s(\mathbf{r},\omega )/{\epsilon }_b,\mu (\mathbf{r},\omega )\equiv {\mu }_s(\mathbf{r},\omega )/{\mu }_b,$$

$$\omega \Gamma (\omega )|\Psi ⟩=\mathbb{L}|\Psi ⟩+\frac{1}{i}|J⟩,$$

$$\Gamma (\omega )\equiv \left[\begin{array}{cc}\epsilon (\mathbf{r},\omega )& 0\\ 0& -\mu (\mathbf{r},\omega )\end{array}\right]\mathrm{and}\mathbb{L}\equiv i{c}_b\left[\begin{array}{cc}0& \nabla \times \\ \nabla \times & 0\end{array}\right].$$

$${\omega }_{\alpha }\Gamma ({\omega }_{\alpha })|{\Psi }_{\alpha }⟩=\mathbb{L}|{\Psi }_{\alpha }⟩,$$

$${\omega }_{\alpha }\equiv {\omega }_{\alpha }^{\prime }+i{\omega }_{\alpha }^{\prime \prime },{\omega }_{\alpha }^{\prime \prime }<0,$$

$$-{\omega }_{\alpha }^{\ast }\Gamma (-{\omega }_{\alpha }^{\ast })|{\Psi }_{\alpha }^{\ast }⟩=\mathbb{L}|{\Psi }_{\alpha }^{\ast }⟩,$$

$$\mathbb{G}(\omega )\to i\sum _{\alpha }\frac{1}{⟨{\Psi }_{\alpha }^{(arb.)}|{\left[\omega \Gamma \right]}_{\alpha }^{\prime }|{\Psi }_{\alpha }^{(arb.)}⟩}\frac{|{\Psi }_{\alpha }^{(arb.)}⟩⟨{\Gamma }_{\alpha }{\Psi }_{\alpha }^{(arb.)}|}{\omega -{\omega }_{\alpha }}+{\mathbb{G}}_{n.r.}(\omega ),$$

$$⟨\mathbf{r}|{\Psi }_{\alpha }⟩\equiv {\Psi }_{\alpha }(\mathbf{r})=\sum _{q=0}^1\sum _{n=1}^{\infty }\sum _{m=0}^n{c}_{q,n,m}^{(\alpha )}{\Phi }_{q,n,m}^{(+)}({\omega }_{\alpha },\mathbf{r})r>R_{\mathrm{out}},$$

$$\begin{array}{rl}{\Phi }_{h,n,m}^{(+)}(\omega ,\mathbf{r})& \equiv {\Phi }_{0,n,m}^{(+)}(\omega ,\mathbf{r})\equiv k^{3/2}{i}^n\left(\begin{array}{c}{\mathbf{M}}_{h,n,m}^{(+)}(k\mathbf{r})\\ -i{\mathbf{N}}_{h,n,m}^{(+)}(k\mathbf{r})\end{array}\right)\\ {\Phi }_{e,n,m}^{(+)}(\omega ,\mathbf{r})& \equiv {\Phi }_{1,n,m}^{(+)}(\omega ,\mathbf{r})\equiv k^{3/2}{i}^{n-1}\left(\begin{array}{c}{\mathbf{N}}_{e,n,m}^{(+)}(k\mathbf{r})\\ -i{\mathbf{M}}_{e,n,m}^{(+)}(k\mathbf{r})\end{array}\right),\end{array}$$

$${\Psi }_{\alpha }(\mathbf{r})=\sum _{q=0}^1\sum _{n=1}^{\infty }\sum _{m=0}^n{d}_{q,n,m}^{(\alpha )}{\Phi }_{q,n,m}^{(1)}({\omega }_{\alpha },\mathbf{r})r<R_{\mathrm{in}},$$

$$\epsilon (\omega )\equiv \frac{{\epsilon }_s(\omega )}{{\epsilon }_b},\mu (\omega )\equiv \frac{{\mu }_s(\omega )}{{\mu }_b},\rho (\omega )\equiv \sqrt{\frac{{\epsilon }_s(\omega ){\mu }_s(\omega )}{{\epsilon }_b{\mu }_b}},$$

$$\begin{array}{rl}{\Phi }_{h,n,m}^{(1)}\left(\omega ,\mathbf{r}\right)& \equiv k^{3/2}{i}^n\left(\begin{array}{c}{\mathbf{M}}_{h,n,m}^{(1)}(k\mathbf{r}\rho (\omega ))\\ -i\sqrt{\frac{\epsilon (\omega )}{\mu (\omega )}}{\mathbf{N}}_{n,m}^{(1)}(k\mathbf{r}\rho (\omega ))\end{array}\right)\\ {\Phi }_{e,n,m}^{(1)}(\omega ,\mathbf{r})& \equiv k^{3/2}{i}^{n-1}\left(\begin{array}{c}{\mathbf{N}}_{n,m}^{(1)}(k\mathbf{r}\rho (\omega ))\\ -i\sqrt{\frac{\epsilon (\omega )}{\mu (\omega )}}{\mathbf{M}}_{n,m}^{(1)}(k\mathbf{r}\rho (\omega ))\end{array}\right),\end{array}$$

$$\Psi_{\alpha}(\boldsymbol{r})=\left\{\begin{array}[c]{c}d_{\alpha}\Phi_{q,n,m}^{(1)}(\omega_{\alpha},\boldsymbol{r})\;\ ;\qquad r<R\\ \;c_{\alpha}\Phi_{q,n,m}^{(+)}(\omega_{\alpha},\boldsymbol{r})\ \ ;\qquad r>R\end{array}\right.\;,$$

$${\gamma }_{\alpha (e,n,\ell )}$$

$${\gamma }_{\alpha (h,n,\ell )}$$

$${\epsilon }_{\alpha }\equiv \frac{{\epsilon }_s({\omega }_{\alpha })}{{\epsilon }_b},{\mu }_{\alpha }\equiv \frac{{\mu }_s({\omega }_{\alpha })}{{\mu }_b},{\rho }_{\alpha }\equiv \sqrt{\frac{{\epsilon }_s({\omega }_{\alpha }){\mu }_s({\omega }_{\alpha })}{{\epsilon }_b{\mu }_b}},$$

$$\Psi_{\alpha}(\boldsymbol{r})=\frac{1}{\mathcal{N}_{\alpha}}\left(\begin{array}[c]{c}\gamma_{\alpha}\Phi_{q,n,m}^{(1)}(\omega_{\alpha},\boldsymbol{r})\;;\qquad r\leq R\\ \ \ \qquad\Phi_{q,n,m}^{(+)}(\omega_{\alpha},\boldsymbol{r})\;;\quad r\geq R\end{array}\right)\;,$$

$$⟨{\Psi }_{\beta }|{\Psi }_{\alpha }⟩\equiv \underset{{\mathcal{S}}_{\mathrm{o}}\to \infty }{\lim }{\int }_{\mathcal{V}}d\mathbf{r}{\left[{\Psi }_{\beta }(\mathbf{r})\right]}^t.{\Psi }_{\alpha }(\mathbf{r})\equiv {\int }_{V_{\infty }}d\mathbf{r}\left\{{\mathbf{E}}_{\beta }(\mathbf{r})\cdot {\mathbf{E}}_{\alpha }(\mathbf{r})+{\mathbf{H}}_{\beta }(\mathbf{r})\cdot {\mathbf{H}}_{\alpha }(\mathbf{r})\right\}.$$

$$⟨{\Psi }_{\beta }|\mathbb{L}{\Psi }_{\alpha }⟩=⟨\mathbb{L}{\Psi }_{\beta }|{\Psi }_{\alpha }⟩.$$

$$\frac{1}{c_b}\left[⟨{\Psi }_{\beta }|\mathbb{L}{\Psi }_{\alpha }⟩-⟨\mathbb{L}{\Psi }_{\beta }|{\Psi }_{\alpha }⟩\right]=0=⟨{\Psi }_{\beta }|\left[k_{\alpha }\Gamma ({\omega }_{\alpha })-k_{\beta }\Gamma ({\omega }_{\beta })\right]|{\Psi }_{\alpha }⟩,$$

$$\begin{array}{rl}\frac{1}{i{c}_b}\left[{\Psi }_{\beta }^t.\mathbb{L}{\Psi }_{\alpha }-{\left(\mathbb{L}{\Psi }_{\beta }\right)}^t.{\Psi }_{\alpha }\right]& =\left[{\mathbf{E}}_{\beta },{\mathbf{H}}_{\beta }\right]\left[\begin{array}{cc}0& \nabla \times \\ \nabla \times & 0\end{array}\right]\left[\begin{array}{c}{\mathbf{E}}_{\alpha }\\ {\mathbf{H}}_{\alpha }\end{array}\right]-\left[{\mathbf{E}}_{\alpha },{\mathbf{H}}_{\alpha }\right]\left[\begin{array}{cc}0& \nabla \times \\ \nabla \times & 0\end{array}\right]\left[\begin{array}{c}{\mathbf{E}}_{\beta }\\ {\mathbf{H}}_{\beta }\end{array}\right]\\ & ={\mathbf{H}}_{\beta }\cdot \nabla \times {\mathbf{E}}_{\alpha }-{\mathbf{E}}_{\alpha }\cdot \nabla \times {\mathbf{H}}_{\beta }-\left[{\mathbf{H}}_{\alpha }\cdot \nabla \times {\mathbf{E}}_{\beta }-{\mathbf{E}}_{\beta }\cdot \nabla \times {\mathbf{H}}_{\alpha }\right]\\ & =\nabla \cdot \left({\mathbf{E}}_{\alpha }\times {\mathbf{H}}_{\beta }-{\mathbf{E}}_{\beta }\times {\mathbf{H}}_{\alpha }\right).\end{array}$$

$$\begin{array}{rl}\frac{1}{c_b}{\int }_{{\mathcal{V}}_a}d\mathbf{r}\left[{\Psi }_{\beta }^t.\mathbb{L}{\Psi }_{\alpha }-{\left(\mathbb{L}{\Psi }_{\beta }\right)}^t.{\Psi }_{\alpha }\right]& =i◯\int {\int }_{{\mathcal{S}}_2}\left[{\mathbf{E}}_{\alpha }(k_{\alpha }\mathbf{r})\times {\mathbf{H}}_{\beta }(k_{\beta }\mathbf{r})-{\mathbf{E}}_{\beta }(k_{\beta }\mathbf{r})\times {\mathbf{H}}_{\alpha }(k_{\alpha }\mathbf{r})\right]\cdot d{\mathcal{S}}_2\\ & -i◯\int {\int }_{\mathcal{S}}\left[{\mathbf{E}}_{\alpha }(k_{\alpha }\mathbf{r})\times {\mathbf{H}}_{\beta }(k_{\beta }\mathbf{r})-{\mathbf{E}}_{\beta }(k_{\beta }\mathbf{r})\times {\mathbf{H}}_{\alpha }(k_{\alpha }\mathbf{r})\right]\cdot d\mathcal{S}.\end{array}$$

$$\underset{S_2\to \infty }{\lim }◯\int {\int }_{{\mathcal{S}}_2}\left[{\mathbf{E}}_{\alpha }(k_{\alpha }\mathbf{r})\times {\mathbf{H}}_{\beta }(k_{\beta }\mathbf{r})-{\mathbf{E}}_{\beta }(k_{\beta }\mathbf{r})\times {\mathbf{H}}_{\alpha }(k_{\alpha }\mathbf{r})\right]\cdot d{\mathcal{S}}_2=0.$$

$$\begin{array}{rl}{\int }_{\mathcal{V}}& d\mathbf{r}{\Psi }_{\beta }^t(\mathbf{r}).\left[k_{\alpha }\Gamma ({\omega }_{\alpha })-k_{\beta }\Gamma ({\omega }_{\beta })\right].{\Psi }_{\alpha }(\mathbf{r})\\ & +i◯\int {\int }_{\mathcal{S}}\left[{\mathbf{E}}_{\alpha }(k_{\alpha }\mathbf{r})\times {\mathbf{H}}_{\beta }(k_{\beta }\mathbf{r})-{\mathbf{E}}_{\beta }(k_{\beta }\mathbf{r})\times {\mathbf{H}}_{\alpha }(k_{\alpha }\mathbf{r})\right]\cdot d\mathcal{S}=0.\end{array}$$

$$\begin{array}{rl}0& =⟨{\Psi }_{\beta }|\Gamma |{\Psi }_{\alpha }⟩=⟨{\Psi }_{\beta }|\Gamma {\Psi }_{\alpha }⟩=⟨\Gamma {\Psi }_{\beta }|{\Psi }_{\alpha }⟩\\ & ={\int }_{V_{\infty }}d\mathbf{r}\left[{\mathbf{E}}_{\beta }(\mathbf{r}),{\mathbf{H}}_{\beta }(\mathbf{r})\right].\left[\begin{array}{cc}\epsilon (\mathbf{r})& 0\\ 0& -\mu (\mathbf{r})\end{array}\right].\left[\begin{array}{c}{\mathbf{E}}_{\alpha }(\mathbf{r})\\ {\mathbf{H}}_{\alpha }(\mathbf{r})\end{array}\right]\\ & ={\int }_{V_{\infty }}d\mathbf{r}\left\{\epsilon (\mathbf{r}){\mathbf{E}}_{\beta }(\mathbf{r})\cdot {\mathbf{E}}_{\alpha }(\mathbf{r})-\mu (\mathbf{r}){\mathbf{H}}_{\beta }(\mathbf{r})\cdot {\mathbf{H}}_{\alpha }(\mathbf{r})\right\},\end{array}$$

$${\int }_{R_{\mathrm{out}}}^{\infty }d\mathbf{r}{\left[{\Psi }_{\beta }^t(\mathbf{r})\right]}^t.\Gamma .{\Psi }_{\alpha }(\mathbf{r})=\sum _{q,n,m}{c}_{\beta (q,n,m)}{c}_{\alpha (q,n,m)}{\int }_{R_{\mathrm{out}}}^{\infty }dr{\Phi }_{q,n,m}^{(+)}(k_{\beta },\mathbf{r}).\Gamma .{\Phi }_{q,n,m}^{(+)}(k_{\alpha },\mathbf{r}),$$

$${\int }_0^{R_{\mathrm{in}}}d\mathbf{r}{\Psi }_{\beta }^t(\mathbf{r}).\Gamma .{\Psi }_{\alpha }(\mathbf{r})=\sum _{q,n,m}{d}_{\beta (q,n,m)}{d}_{\alpha (q,n,m)}{\int }_0^{R_{\mathrm{in}}}d\mathbf{r}{\left[{\Phi }_{q,n,m}^{(1)}({\rho }_{\beta }{k}_{\beta },\mathbf{r})\right]}^t.\Gamma .{\Phi }_{q,n,m}^{(1)}({\rho }_{\alpha }{k}_{\alpha },\mathbf{r}),$$

$${\int }_0^{R}d\mathbf{r}\mu {\mathbf{H}}_{e,n,m}(k_{\beta }\mathbf{r})\cdot {\mathbf{H}}_{e,n,m}(k_{\alpha }\mathbf{r})=-\frac{z_{\alpha }^{3/2}{z}_{\beta }^{3/2}{\gamma }_{\alpha (e,n)}{\gamma }_{\beta (e,n)}}{{\mathcal{N}}_{\alpha }{\mathcal{N}}_{\beta }}\epsilon \left\{{\int }_0^1{\tilde{r}}^2{j}_n(\rho z_{\alpha }\tilde{r})j_n(\rho z_{\beta }\tilde{r})d\tilde{r}\right\}$$

$$=\frac{z_{\alpha }^{3/2}{z}_{\beta }^{3/2}{\gamma }_{\alpha (e,n)}{\gamma }_{\beta (e,n)}}{{\mathcal{N}}_{\alpha }{\mathcal{N}}_{\beta }}\epsilon \left\{\frac{z_{\beta }{j}_n(\rho z_{\alpha })j_n^{\prime }(\rho z_{\beta })-z_{\alpha }{j}_n^{\prime }(\rho z_{\alpha })j_n(\rho z_{\beta })}{\rho \left(z_{\alpha }^2-z_{\beta }^2\right)}\right\}$$

$$=\frac{z_{\alpha }^{3/2}{z}_{\beta }^{3/2}}{{\mathcal{N}}_{\alpha }{\mathcal{N}}_{\beta }}\frac{z_{\alpha }{h}_n^{\prime }(z_{\alpha })h_n(z_{\beta })-z_{\beta }{h}_n^{\prime }(z_{\beta })h_n(z_{\alpha })}{z_{\alpha }^2-z_{\beta }^2},$$

$${\int }_R^{\infty }d\mathbf{r}{\mathbf{H}}_{e,n,m}(k_{\beta }\mathbf{r})\cdot {\mathbf{H}}_{e,n,m}(k_{\alpha }\mathbf{r})⟶-\frac{z_{\alpha }^{3/2}{z}_{\beta }^{3/2}}{{\mathcal{N}}_{\alpha }{\mathcal{N}}_{\beta }}\underset{\eta \to 0}{\lim }{\int }_1^{\infty }{\tilde{r}}^2{h}_n(z_{\alpha }\tilde{r})h_n(z_{\beta }\tilde{r})e^{-\eta {\tilde{r}}^2}dr$$

Full text

1 Introduction

The eigenstates of open systems, known as Resonant States (RSs), can be viewed as a means of simultaneously describing both the spatial and temporal behavior of the oscillatory ‘eigenmodes’ of an open system. As such, RSs have been a fruitful, albeit problematic, concept in many branches of physics for over a century, [1] where they have been given a wide variety of names like: quasi-normal modes, transient modes, and the famous leaky modes of wave-guides. The RS concept was widely developed in quantum mechanics [2] after its pioneering success by Gamow in describing alpha decay [3], and later by Siegert for characterizing nuclear reactions [4]. Throughout the 20th century, quantum RS descriptions continued to draw the attention of prominent physicists like Wigner [5], Peierls [6], and Zel’dovich [7, 8], the latter of whom described the mathematical difficulties associated with RSs as follows: “An exponentially decaying state that describes, for example, the phenomenon of a decay, is characterized by a complex value of the energy, the imaginary part of the energy giving the decay probability. The wave function of this state increases exponentially in absolute value at large distances, and therefore the usual methods of normalization, of perturbation theory, and of expansion in terms of eigenfunctions do not apply to this state”. For us, the pivotal words in Zel’dovich’s remark are ‘usual methods’. Indeed, this work shows that ‘non-traditional’ normalization methods allow the RSs to quantitatively express open system spectral response functions.

The utility of RSs is not limited to the quantum realm, and in the 1990s Leung et al.[9, 10] studied various mathematical aspects of the RS states in ‘classical field’ applications, renaming them “Quasi-Normal-Modes” (QNMs), a terminology often adopted in photonics and colliding black holes.[11] It thus appears that RSs are a nearly universal tool for describing open systems in arbitrary spatial dimensions, and we restrict attention here to 3D electromagnetic systems for the purpose of discussion.

Combining the assumption that energy outside a system is ‘radiated’ into a background media governed by a Helmholtz type equation satisfying outgoing boundary conditions, leads to frequency domain eigenstates characterized by complex eigenfrequencies, $\omega_{\alpha}=\omega_{\alpha}^{\prime}+i\omega_{\alpha}^{\prime\prime}$, with causality constraining $\omega_{\alpha}^{\prime\prime}$ to be negative (given that we adopted an $\exp(-i\omega t)$ inverse temporal Fourier transform/time harmonic convention). In terms of wavenumber, $k_{\alpha}=\omega_{\alpha}/c$, the associated far-field radial dependence is proportional to $\exp\left[ik_{\alpha}(r-ct)\right]/k_{\alpha}r$, giving RS excitations a physically desirable exponential temporal decay, but necessarily accompanied by exponentially increasing spatial oscillations at large distances, (sometimes referred to as the ‘exponential catastrophe’ [12]). An example of the generality of the RS approach is found in electromagnetic scattering problems [13, 14, 15], where one encounters both vector and scalar wave RS problems. Mathematical descriptions of vector RSs are more complex than those for scalar fields, and this was a major motivation for this work (as well as refs. [16, 17, 18], where the mathematical groundwork and physical motivation are developed).

For a long time, the mathematical difficulties posed by the spatial divergence of RSs led to a widely held belief they were principally of phenomenological interest, but it has become increasingly appreciated in recent years that RSs can be used as the basis for quantitative calculations. Notably, causality considerations regularize the exponentially diverging behavior when one transforms from frequency domain descriptions back into the time domain.[1, 19, 20, 21, 22] Determining the correct normalization of the RSs in the frequency domain long remained a matter of debate with a number of recent papers in photonics advocating various RS normalization schemes, many of which involve a regularization of the RS inner product integrals (see refs. [23, 24, 25, 26, 27, 28] and articles cited in two recent reviews [29, 30]). A notable feature of the RS paradigm is that normalization adjusts the eigenstate phases in addition to their amplitudes.

For the propose of discussion, normalization schemes can be classified into three principal categories. A first technique relies on the use of perfect matched layers (PML) to regularize the normalization integral [23, 27, 29], which can be interpreted as a deformation of spatial integration paths into the complex plane. This technique is strongly linked to the complex scaling approach developed in the context of quantum scattering theory [2, 31]. Two other methods make use of a surface integral to complement the volume integral  [24, 25, 26, 30]. The definitions for the surface integral in these two approaches differ however and they therefore lead to two different definitions for the norm. The first definition [24, 30] has been inspired by an earlier work of Lei et al. [32] which makes use of the Silver-Muller radiation condition that describes the asymptotic behavior of outgoing waves to evaluate the surface term. The volume integral then needs to be calculated over a volume that is sufficiently large for the surface integral to be located in the far-field region. Another definition of the norm integral makes use of analytic continuation of the RS field along with some properties of vector analysis to express the surface term.[25, 26]

Here we illustrate yet another regularization scheme, along the lines first proposed by Zel’dovich [7], to regularize the inner product by introducing a Gaussian ‘killing’ function, $\exp(-\eta r^{2})$ into RS inner product integrals and then taking the limit $\eta\rightarrow 0$ after integration (a technique which nowadays is considered as belonging to the general distribution techniques of regularization by ‘good functions’ [33], and which agrees with the old technique of fractional calculus dating back to Euler and Leibnitz, and with that of Mellin transforms [34]). One advantage of this technique is that it rigorously leads to analytic expressions for RS type integrals in the exterior regions where the field divergences occur, (shown in a 1992 paper by one of the present authors and two colleagues citeMDS). This technique was recently extended to 3D electromagnetic RSs where it was referred to somewhat picturesquely (or melodically) as “Killing Mie Softly” with the derivations involving distribution theory mathematics by two of the present authors [16]. Besides its value in replacing numerical integrals with analytic formulas, we will see that this approach also sheds new light on the other regularization schemes, and helps justify their physical justification and equivalence.

This article is organized as follows: section 2 reviews our RS formulation of electromagnetic scattering problems followed by a presentation of RS multipolar representations. In section 2.3, we derive orthogonality relations for the RSs of scatterers possessing electric and/or magnetic dispersion. We then show that the analytical formulas for Gaussian regularized multipole integrals given in Ref.[16] and Appendix C, verify the general orthogonality formulas in the particular case of spherical scatterers. Next, we derive resonant state normalization formulas for fully dispersive scattering materials, followed by a Gaussian multipole application in the case of electrically and magnetically dispersive spheres. These analytic formulas open the door to RS studies of meta-materials and illustrate electromagnetic duality in a manner that was not achieved previously. The physical significance of these results and their relationship with other formulations and results in the literature are also discussed in this section, with additional comparisons with the literature being detailed in Appendix A.

The application of RS spectral expansions to Mie theory, which we have developed and exploited for a number of years now[17, 18, 22, 35, 36, 37], is reviewed in section 3. Advantages of this approach is that it enables RS expansions to remain accurate even far from the resonant frequencies by including the crucial non-resonant terms in the RS expansion. Section 4 is dedicated to high precision numerical results and applications for dispersive scatterers. These are intended to serve as benchmarks for those using our methods or elaborating purely numerical models for nanophotonic and metamaterial problems. The conclusion discusses some of the avenues for future investigations and applications. The appendices A-G, elaborate on the mathematical definitions and derivations presented in the main text.

2 Resonant states

When formulating Resonant States (RSs) in electromagnetism, it is convenient to express the electromagnetic field as a six component ‘ket’ state, $\left|\Psi\right\rangle$, and appropriately dimensioned source currents, $\left|J\right\rangle$, as multi-component fields,

[TABLE]

where $\boldsymbol{E}$ and $\boldsymbol{H}$ are respectively the vector electric and magnetic fields. We found it convenient to adopt, ‘field theoretic’ units of m*-3/2*, but Gaussian units or impedance adjusted SI units are acceptable alternatives since these also assign the same units to both electric and magnetic fields (although we found that m*-3/2* field dimensions are particularly well adapted to Green’s functions, energy density, LDOS [38] applications). Six component electromagnetic field formulations are nothing new, and notably $\left\langle\boldsymbol{r}|\Psi\right\rangle$ in our notation is quite similar (but not quite identical) to the 6-component field, $\overrightarrow{\mathbb{F}}(\boldsymbol{r})$, that appeared in an RS study during the writing of an initial version of this manuscript.[39]

The scattering particles are taken to be characterized by spatially local materials with sharp boundaries and temporally dispersive constitutive parameters, $\varepsilon(\boldsymbol{r},\omega)$, and/or $\mu(\boldsymbol{r},\omega)$. The scatterers are assumed to be immersed in an isotropic non-dispersive background medium, described by real-valued constitutive parameters $\varepsilon_{b}$ or $\mu_{b}$, with the background media wave velocity given by $c_{b}=c_{\rm v}\sqrt{\varepsilon_{b}\mu_{b}}$, where $c_{\rm v}$ is the vacuum speed of light. Henceforth, the constitutive parameters of the scatterers will be consistently normalized with respect to the background medium properties,

[TABLE]

which notably simplifies various derivations, and renders almost all further formulas independent of unit conventions.

Using the above definitions and normalizations, the frequency domain Maxwell equations can be written as a single convenient equation,

[TABLE]

Resonant states are defined as source-free eigensolutions of Eq.​​ 3b,

[TABLE]

obeying outgoing boundary conditions with complex eigenfrequencies,

[TABLE]

with $\omega_{\alpha}^{\prime\prime}<0$ a necessary condition for exponential temporal decay.

The frequency dependence of the medium metric, $\Gamma(\omega)$, is required to satisfy Kramers-Kronig relations and $\Gamma^{\ast}(\omega)=\Gamma(-\omega^{\ast})$, so that fields in the time domain are both causal and real valued. Consequently, even though negative frequency RSs must be included in our analysis, they are not independent of the positive frequency solutions. One remarks that Eq.​​ 4 is a self-consistent eigenvalue equation in terms of the complex wavenumber (frequency), and consequently the determination of RS eigenvalues must almost always be carried out numerically.

Taking the complex conjugate of Eq.​​ 5 for a given eigenstate $\omega_{\alpha}$, and remarking that $\mathbb{L}^{\ast}=-\mathbb{L}$, one obtains,

[TABLE]

which is of the same form as Eq.​​ 4, so that for any RS eigenvalue, $\omega_{\alpha}$, there is an associated RS with eigenvalue, $-\omega_{\alpha}^{\ast}$, which also satisfies outgoing boundary conditions.

Since the RS frequencies are discrete, their index, $\alpha$, can be assigned integer values, but when symmetries permit additional quantum numbers, following the discussion in section 2.1 below, we will designate RS indices by, $\alpha(q,...,\ell)$ where $q,...$, denotes symmetry based quantum numbers, while the $\ell$ number adopts both positive and negative integer indices, $\ell=\pm\left((0),1,2,\ldots\right)$ ; the presence of an $\ell=0$ index depends on the mode symmetry properties. From the discussion of the previous paragraph, we assign the $\ell$ indices such that, $\omega_{\alpha(q,...;-\ell)}=-\omega^{\ast}_{\alpha(q,...,\ell)}$, which requires all $\omega_{\alpha(q,...;0)}$ to have purely imaginary values. In the cases considered here, there will generally be at most one RS eigenvalue on the negative imaginary axis for a given set of quantum numbers, so we can denote such states $\omega_{\alpha(q,...,\ell=0)}$ when they occur.

A Green function response operator, $\mathbb{G}$, by definition, produces the electromagnetic field equations of Eq.​​ 3b when acting on arbitrary source currents, $\left|J\right\rangle$. Adopting a non-conventional ‘bra’ and ‘ket’ notation discussed below in section 2.3, and using a first order Taylor expansion of $\Gamma(\omega)$ about $\Gamma(\omega_{\alpha})$, the above RS formalism tells us that the spectral expansion of the Green function operator in the frequency domain takes the form,

[TABLE]

which employs the shorthand notations, $\Gamma_{\alpha}\equiv\Gamma(\omega_{\alpha})$ and $\left[\omega\Gamma\right]_{\alpha}^{\prime}\equiv\frac{d}{d\omega}\left[\omega\Gamma(\omega)\right]_{\omega=\omega_{\alpha}}$. The operator, $\mathbb{G}_{\rm n.r.}(\omega)$, in Eq.​​ 7 is taken to include contributions like static poles [40] and other non-resonant state contributions associated with the non-uniqueness of the Green operator and boundary conditions.

The state, $|\Psi_{\alpha}^{\rm(arb.)}\rangle$, in Eq.​​ 7 stands for a state which is ‘arbitrarily’ normalized (in the sense that neither Eq.​​ 4 nor Eq.​​ 7 determines how RSs should be normalized). The Green’s function can however be made to ‘determine’ the RS normalization if we assign a dimensionless predetermined number to $\langle\Psi_{\alpha}^{\rm(arb.)}|\left[\omega\Gamma\right]_{\alpha}^{\prime}|\Psi_{\alpha}^{\rm(arb.)}\rangle$ (a value of ‘$1$’ appearing to be the ‘natural’ choice since this leads to a particularly simple expression for $\mathbb{G}$ in Eq.​​ 7). We will adopt a different normalization choice, in section 2.4, which facilitates the expressions of linear response theory which possibly give the RS a more physically intuitive interpretation.

In the following sections, we first show how RSs can be developed on a multipole basis, and how RS product integrals can be evaluated analytically via Gaussian regularization. This methodology is then used to validate a number of RS properties in this section, finally arriving at analytic RS orthogonality and normalization formulas in the case of spherical scatterers.

2.1 Multipolar expansions of resonant states

In any homogeneous region of the scattering system, any three dimensional RS can be expanded in terms of multipole wave functions evaluated at the RS frequency, $\omega_{\alpha}$. Specifically, let us consider the generalizable example of a single scattering particle as shown in Fig.​​ 1. In the homogeneous background media lying outside a sphere of radius $R_{\rm out}$, a RS field, $\Psi_{\alpha}(\boldsymbol{r})$, can be expanded in terms of multipolar basis functions, $\Phi_{q,n,m}^{(+)}(\omega,\boldsymbol{r})$:

[TABLE]

where $c_{q,n,m}^{(\alpha)}$ are the RS expansion coefficients.

The six component functions, $\Phi_{q,n,m}^{(+)}(\omega,\boldsymbol{r})$, can be characterized according to being either of magnetic ($h$) field types (denoted by a $q=0$ index), or of electric ($e$) field types (denoted by the index $q=1$). The other indices are the multipole (angular momentum) numbers, $n\geq 1$, and the azimuthal (projection) numbers $m<n$. Using multipole analysis, the $\Phi_{q,n,m}^{(+)}$ functions can be defined as:

[TABLE]

where $k\equiv\omega/c_{b}$ is the exterior medium wavenumber, while $\boldsymbol{M}^{(+)}_{q,n,m}(k\boldsymbol{r})$ and $\boldsymbol{N}^{(+)}_{q,n,m}(k\boldsymbol{r})$, are the dimensionless outgoing Vector Partial Waves (VPWs), (defined in Appendix B with additional details and references).

If one chooses the origin of the coordinate system to lie inside a scattering particle of interest, one can determine a homogeneous spherical region of radius $R_{\rm in}$ containing only the scattering medium. Inside this region, the RS fields can be expressed,

[TABLE]

where the real-valued RS expansion coefficients are designated $d_{q,n,m}^{(\alpha)}$ and the $\Phi_{q,n,m}^{(1)}$ are formulated in terms of the ‘regular’ (i.e. singularity free) VPWs, and the normalized (possibly dispersive) constitutive parameters of the scattering media:

[TABLE]

so that the regular, medium dependent media wave functions, $\Phi_{q,n,m}^{(1)}$ can be defined as:

[TABLE]

where $\boldsymbol{M}_{n,m}^{(1)}$ and $\boldsymbol{N}_{n,m}^{(1)}$ are regular (i.e. singularity free) VPWs.

One way to solve the RS eigenstate values, $\omega_{\alpha}$, would be to numerically solve the propagation of each $\Phi^{(1)}_{q,n,m}(\omega,\boldsymbol{r})$ function through the inhomogeneous region, $R_{\rm in}<r<R_{\rm out}$, and then determine the values of $\omega_{\alpha}$ which allow sets of coefficients $c_{q,n,m}^{(\alpha)}$ and $d_{q,n,m}^{(\alpha)}$ to be determined which satisfy both electric and magnetic boundary values at $R_{\rm in}$ and $R_{\rm out}$. An example of this procedure for spherical particles is given in section 2.2 below, but for particles of arbitrary shape, alternative schemes, employing other types of basis functions are generally preferred (cf. articles cited in [29]).

The generality of the multipole technique stems from the fact that regardless of the means by which RS eigensolutions have been obtained, their solutions in homogeneous regions can be reexpressed in terms of the spherical wave basis described in this section (using for example the techniques described in ref. [41]). As mentioned earlier, solving an eigenvalue problems like Eq.​​ 4, only determines the coefficients $c_{q,n,m}^{(\alpha)}$ and $d_{q,n,m}^{(\alpha)}$ up to an overall ‘normalization’ factor, henceforth designated $\mathcal{N}_{\alpha}$, which must be determined by additional criteria which will be specified in section 2.4.

2.2 Resonant states of spherical particles

For spherically symmetric particles, there is no mixing of the multipole indices and for each distinct set of multipole numbers ($q$, $n$, $m$), there exists an infinite discrete set of resonant states (degenerate with respect to the $m$ number). As discussed after Eq.​​ 6 above, an additional ‘quantum’ index, $\ell$, enumerates the different RS eigenfrequencies, $\omega_{\alpha}$ within a given multipole set so that $\alpha(q,n,m,\ell)$ designates one and only one RS. Due to the symmetries, the RS multipole developments of Eqs.​ 8 and 10 simplify here to:

[TABLE]

with the RS index, $\alpha$, uniquely identified by the full set of quantum numbers, i.e. $\alpha(q,n,m,\ell)$.

For spherically symmetric cases, the overall normalization factor, $\mathcal{N}_{\alpha}$, can be associated with the external field coefficient as $c_{\alpha}=1/\mathcal{N}_{\alpha}$. Given the expressions in Eq.​​ 9 and Eq.​​ 12 for the $\Phi^{(+)}$ and $\Phi^{(1)}$ basis functions, the continuity of the transverse electric and transverse magnetic fields [42] respectfully take the form of analytic linear relationships between internal field and external coefficients, which we henceforth write as, $\gamma_{\alpha}\equiv d_{\alpha}/c_{\alpha}$:

[TABLE]

where $z_{\alpha}\equiv k_{\alpha}R$, is the complex size parameter, $j_{n}(z)$ the spherical Bessel functions, and $h_{n}(z)$ the outgoing spherical Hankel functions while $\psi_{n}(z)\equiv zj_{n}(z)$ and $\xi_{n}\equiv zh_{n}(z)$ are their respective Ricatti-Bessel function counterparts (cf. Appendix B). Since the respective pairs of continuity conditions for $\gamma_{\alpha}$ in Eq.​​ 14 are different from one another at arbitrary frequencies, the RS frequencies are defined as being the set of discrete frequencies for which both conditions on the respective $\gamma_{\alpha}$ are satisfied.

The notations $\varepsilon_{\alpha}$ and $\mu_{\alpha}$, and $\rho_{\alpha}$, introduced in Eq.​​ 14 are shorthands for designating the constitutive parameters evaluated at the RS frequency, i.e.,

[TABLE]

which are normalized with respect to the exterior medium in accordance with the notation introduced in Eq.​​ 3b. To summarize, the exact expressions for the RSs of spherical particles can henceforth be written,

[TABLE]

The value of the normalization factor, $\mathcal{N}_{\alpha}$, will be determined analytically in section 2.4 for the spherical scatterer case. Alternative compact expressions for the spherical particle RS expression of Eq.​​ 16 and comparisons with expressions in the literature are given in Appendix A.

One can readily verify that satisfying the equalities in Eq.​​ 14 is identical to the condition for the occurrence of poles in the electric and magnetic Mie coefficients respectively (as can be seen by examining section D.2 of Appendix D). We will show in section 3 that the determination of the values of $z_{\alpha}$ is almost all that we need in order to perfectly reconstruct Mie response theory, once the normalization factors, $\mathcal{N}_{\alpha}$, have been determined as functions of $\omega_{\alpha}$ in Eq.​​ 37.

2.3 Inner products and RS orthogonalization

In conventional quantum mechanics, inner products are defined using complex conjugated ‘bra’ states, but in a lossy system this would transform loss to gain which is something we need to avoid[38, 43]. We thus define ‘bra’ states as, $\left\langle\Psi_{\alpha}\right|\equiv\left\langle\boldsymbol{E}_{\alpha},\boldsymbol{H}_{\alpha}\right|$, without complex conjugation of the fields. The inner product of any two electromagnetic states, $\Psi_{\alpha}$ and $\Psi_{\beta}$ is thus defined as the integral over a volume, $\mathcal{V}$, inside a surface $\mathcal{S}_{\rm o}$ that is sent to infinity:

[TABLE]

Given the fundamental exponential divergence of the RSs in the far-field discussed in the introduction, it could appear that the RS inner product integral of Eq.​​ 17 is ‘ill-defined’, but the Gaussian regularization method assigns unambiguous, finite, values to RS products like those of Eq.​​ 17 in a manner which allows the differential operator, $\mathbb{L}$ of Eq.​​ 3b, to remain symmetric when acting on regularized RS products, i.e.,

[TABLE]

Applying the RS equations of motion in Eq.​​ 3b to this relation provides an orthogonality relation for the product states, valid even in the case of temporally dispersive media,

[TABLE]

where $\Gamma(\omega)$ is the frequency dependent medium metric.

Although the inner product orthogonality relation of Eq.​​ 19 is defined as an integration over all space, it can alternatively be reexpressed as an integration over a finite volume supplemented by an appropriate surface integral. This is done by remarking that vector identities allow the inner product integrands of the left hand side of Eq.​​ 19 to be rewritten:

[TABLE]

Next, let us consider an arbitrary ‘annular’ closed volume, ${\mathcal{V}_{a}}$, with an inner surface, $\mathcal{S}$, exterior to the scattering system and an outer surface, $\mathcal{S}_{2}$, exterior to $\mathcal{S}$ as illustrated in Fig.​​ 2.

Integrating over the annular volume integral and applying the divergence theorem to Eq.​​ 20 yields:

[TABLE]

Letting the outer surface tend to infinity, the result can only be consistent with Eq.​​ 19 provided that,

[TABLE]

The vanishing of the surface integral, $\mathcal{S}_{2}$, at infinity arises from distribution theory and is discussed at length in ref.[16] and reviewed in Appendix C, but it can be intuitively understood as a consequence of the fact that the Gaussian regularization program suppresses the fields at distances far from the origin thanks to an $e^{-\eta r^{2}}$ factor.

Eq.​​ 22 allows the orthogonality relation of Eq.​​ 19 to alternatively be expressed as a volume integral over any finite closed region, $\mathcal{V}$, containing the system, supplemented by an integral on the surface, $\mathcal{S}$, of $\mathcal{V}$:

[TABLE]

For lossless particles, the medium metric, $\Gamma$, is frequency independent and real valued, so that Eq.​​ 19 reads $(k_{\alpha}-k_{\beta})\left\langle\Psi_{\beta}\right|\Gamma\left|\Psi_{\alpha}\right\rangle=0$, which for $k_{\alpha}\neq k_{\beta}$ leads to a frequency independent eigenstate orthogonality condition for lossless systems:

[TABLE]

which can of course also be expressed in the finite volume-surface form of Eq.​​ 23.

It is worth remarking at this stage that the presence of the medium metric, $\Gamma$, gives a Lagrangian aspect to RS field products (as opposed to the more common Hamiltonian type products), which will become even more apparent when we consider RS normalization in section 2.4.

2.3.1 Regularizing multipole expansion products

Once the RS fields are developed according to Eq.​​ 8, then the analytic expressions for RS product integrals, given in Appendix C, allow one to rapidly evaluate all the product integrals required for orthogonalization and normalization in the problematic far-field region. Concretely, for an arbitrary pair of resonant states, $\Psi_{\alpha}$ and $\Psi_{\beta}$, developed according to Eq.​​ 8, the volume integral of RS product in the region $r>R_{\rm out}$ (cf. Fig.​​ 1) is:

[TABLE]

where the orthogonality of the VPWs over angular integration results in there being only a single sum over basis indices rather than the double sum that would generally occur using other basis functions. Analytic expressions for all the integrals on the right hand side of Eq.​​ 25 are given in Eq.​​ 88 and Eq.​​ 89 of Appendix C for $\alpha\neq\beta$ and $\alpha=\beta$ respectively.

Inside a homogeneous region, $r<R_{\rm in}$, lying completely inside a scattering object (cf. Fig.​​ 1), one can again develope the RS field on a multipole basis according to Eq.​​ 10, and the RS volume integrals are:

[TABLE]

with ordinary analytic expressions for the integrals on the right hand side of Eq.​​ 26 being given in Eq.​​ 90 and Eq.​​ 91 of Appendix C. Numerical integration is still required for volume product integrals of RS fields in heterogeneous regions, like the region $R_{\rm in}<r<R_{\rm out}$ in Fig.​​ 1, but for spherical scatterers, one can take $R_{\rm in}=R_{\rm out}=R$, and the inner product integrals for RS orthogonalization and normalization of become entirely analytic as we demonstrate below.

2.3.2 Orthogonalization : Example with a spherical scatterer

Let us consider the case of a lossless sphere with relative permittivity, $\varepsilon=16$ and null permeability contrast, $\mu=1$. A few low frequency electric mode RSs are graphically represented in Fig.​​ 6 and given to high numerical precision in Table 1. Since the sphere is lossless, the frequency independent RS orthogonality relation of Eq.​​ 24 applies. This orthogonality relation requires that either the electric field and magnetic product integrals exactly cancel one another, or that they are both null. However, we know that for a system with material losses, knowledge of the electric field completely determines the magnetic field and vice versa, which leads us to the conclusion that RS orthogonality for dispersionless scatterers corresponds to the magnetic and electric field product integrals being independently zero. Rather than mathematically proving this statement, we next invoke the formulas of Appendix C to show that this is indeed true for a spherical scatterer.

Let us consider a magnetic field product of two distinct electric type RSs ($q=0$) of a dispersionless spherical scatterer. Rescaling the radial variable by the sphere radius, $\tilde{r}\equiv r/R$, a volume integral of the (complex valued) magnetic fields inside the sphere can then be evaluated as,

[TABLE]

where Eq.​​ 27a is obtained after analytical integration of the angular variables, while the term in brackets of (27b) is a direct application of the bracketed integral in Eq.​​ 27a as given by Eq.(85a) of [16] (cf. also Watson ref. [44], and Eqs.​ 85 and 90c). The final result of Eq.​​ 27c exploits the expressions of Eq.​​ 14a for $\gamma_{\alpha(e,n)}$, and employs algebraic manipulations that are only true when $z_{\alpha}$ and $z_{\beta}$ are distinct RS size parameter eigenvalues (i.e. $z_{\alpha}\neq z_{\beta}$).

Thanks to the Gaussian regularization program, we find an analytic expression for integral in the region outside the sphere (where $\mu(r)=1$ due to medium normalization),

[TABLE]

where we used Eq.(103a) in ref. [16] to evaluate the integral in the first line of this equation (reproduced in Eq.​​ 89b). The ‘Killing’ regularization, described in detail in ref. [16] and Appendix C, tames the integrand’s exponentially divergent amplitude with a $e^{-\eta r^{2}}$ factor in order to yield the finite analytic expression of Eq.​​ 28b in the $\eta\rightarrow 0$ limit. We ignored irrelevant overall sign factors in the above demonstration because it is now clear that the integrals in Eqs.​ 27c and 28b exactly cancel, so that indeed there is a magnetic field orthogonality of RS states in spheres for dispersion free media to yield,

[TABLE]

where we dropped both the multiple indices, $n$ and the field type index ‘$e$’ $(q=0)$ in Eq.​​ 29 since as shown in Appendix C that analogous orthogonality relations holds for magnetic, (‘$h$’ $(q=1)$), type RS modes.

In order to appreciate the complexity of what is accomplished in Eq.​​ 29, the real and imaginary parts of a RS magnetic product integrand of the first two electric dipole RSs at $z_{\alpha(e,n=1,\ell=1)}\simeq 1.0395-i0.5009$ and $z_{\beta(e,n=1,\ell=2)}\simeq 1.0527-i0.072355$ of a high index dielectric sphere ($\varepsilon=16$) are plotted in Fig.​​ 3 (cf. Table 1 and Fig.​​ 6). Specifically, the real and imaginary parts of the magnetic field product integrands given in Eq.​​ 27a ($r<R$) and Eq.​​ 28a ($r>R$) are plotted first for $r/R$ in Fig.​​ 3(a), up to $1.2$ and then extended in Fig.​​ 3(b) up to values of $r/R=25$, where a Gaussian killing factor of $\eta=0.025$ is adopted for values of $r>R$. High accuracy numerical calculations confirm that the integrals of both the real and imaginary parts of the total magnetic field integral do indeed tend towards zero when a sufficiently small value of $\eta$ in the region $r/R>1$ is adopted.

Analogous result holds for RS orthogonality in terms of an electric field despite presence of discontinuities in the electric field component normal to the scatterer surface. Analytical formulas for electric field orthogonality are given explicitly in Appendix E and one finds,

[TABLE]

when $\alpha\neq\beta$, a result which is illustrated graphically by plotting the electric field integrands of Eq.​​ 97a and Eq.​​ 97d for $r/R$ up to $1.2$ and Fig.​​ 3(d) for $r/R$ up to $25$. As shown in Appendix E, the mathematics for the electric field integration was somewhat more complicated than it was for the magnetic field integrals of Eq.​​ 28b since the required integrals were not given in standard references like that of Watson [44], but they are provided in reference [16].

An advantage of Gaussian regularization is that it can replace numerical integrations by analytical expressions (once the multipole development of the RS fields outside the system are determined). Although Gaussian regularization could also be used as numerical method, this may not be the most efficient choice for the following reason: a numerical evaluation of a Gaussian regularized RS integral should adopt a small value of $\eta$ , so as to be close to the $\eta\to 0$ limit, but small values of $\eta$ lead to large amplitude far-field oscillations (cf. Fig.​​ 4), which renders numerical integrations time consuming. For numerical evaluations, it is often preferable to generalize the radial coordinate to complex values and deform radial integration path so that it goes to infinity at a finite angle from the real axis in the complex plane. This procedure is commonly employed in temporal analysis of applied mathematics and lies at the heart of the PML method.[45]

We can readily understand the numerical advantage of the complex contour method from Fig.​​ 4, which plots the real part of Eq.​​ 28, i.e. ${\rm Re}\left[r^{2}h_{n}(z_{\alpha}r)h_{n}(z_{\beta}r)e^{-\eta r^{2}}\right]$ (with $\eta=0.025$) in the complex plane. The Gaussian regularized contour on the real axis is plotted as a thick blue line, while the red dashed curve indicates the path of a contour rotated by an angle of $\pi/4$ into the complex plane. The integrals along the two contours must be the same, but the red dashed PML type contour is clearly more manageable numerically and indeed a Gaussian regularization is superfluous for the rotated path which could be evaluated directly without regularization (i.e. for $\eta=0$). Although it is visually clear from Fig.​​ 3 and Fig.​​ 4 that RS product integrands oscillate around zero outside the scatterer, they do not quite average to zero since the RS orthogonality relation of Eq.​​ 29 requires the exterior product integral to be equal and opposite to the field product integral inside the sphere, (this interior integral is generally not null, as seen in this example by inspection of Fig.​​ 3(a) and Fig.​​ 3(b)).

The utility of analytic Gaussian regularization formulas is further highlighted by remarking that if one also evaluates the electric field product integrals (cf. Appendix E in Appendix E), then after some algebra one finds, when $k_{\alpha}\neq k_{\beta}$, that:

[TABLE]

where the first term is a volume integral of fields inside a sphere of radius, $R$, denoted $\mathcal{V}_{R}$, and the second term is a quantity evaluated at the sphere’s surface. It is worth remarking that although $R$ can be ‘naturally’ taken to be the particle radius, Eq.​​ 31 remains valid for any sphere surrounding the scatterer.

An examination of Eq.​​ 31 shows that the volume integral and the second, ‘surface’, contribution are none other than special cases (spherical volume and surface) of the more general volume and surface integrals derived in Eq.​​ 23. Although the explicit agreement between the analytic Gaussian regularization method with Eq.​​ 23 was only shown here for the case of a dispersionless sphere, analogous manipulations confirm this agreement for the case of dispersive spheres, at the cost of some additional algebra.

2.4 Resonant state normalization

A variety of methodologies have been adopted to determine RS normalization [24, 25, 26, 27, 46], which has given rise to ‘different’ normalization conditions being proposed for photonic problems, but Kristenssen et al.[24] and others have argued that varying normalization formulas should be at least inter-consistent. We agree that final response functions of an RS theory should be independent of regularization schemes, but as we saw that if we adopted Eq.​​ 7 as Green function response, the ‘correctness’ of the RS normalization factors can be considered a matter of convention.

We found it advantageous to adopt a field theory based route to normalization and impose the RS norms in Eq.​​ 7 be proportional to their associated RS frequencies, $\omega_{\alpha}$, which given the adimensional requirements scalar product leads to:

[TABLE]

where we recall that $z_{\alpha}=\omega_{\alpha}R/c$. Since the frequency derivatives of the constitutive factors in $d\left[\omega\varepsilon(\boldsymbol{r},\omega)\right]/d\omega$ and $d\left[\omega\mu(\boldsymbol{r},\omega)\right]/d\omega$ can be interpreted as accounting for the energy associated with material degrees of freedom [47], one can interpret the Lagrangian energy of the RS as $\mathcal{E}_{{\alpha}}^{\rm L}=\hbar c\left\langle\Psi_{\alpha}\right|\frac{d}{d\omega}(\omega\Gamma)\left|\Psi_{\alpha}\right\rangle/R=\hbar\omega_{\alpha}$. The Lagrangian interpretation of Eq.​​ 32, is more transparent in dispersionless media where it simplifies to:

[TABLE]

Despite the divergence of RS field amplitudes, the infinite volume integrals in Eq.​​ 32 can be regularized by either Gaussian regularization or by a PML type rotation of the integration path into the complex plane as illustrated in Fig.​​ 4. Yet another way to determine the normalization factors is said to avoid the RS divergences by integrating only over a finite region of space, analogous to the formula found for RS orthogonalization in section 2.3. Indeed, it suffices to divide both sides of Eq.​​ 23 by $k_{\alpha}-k_{\beta}$, and then take the limit $k_{\alpha}\to k_{\beta}$ to obtain a finite volume expression of $\left\langle\Psi_{\alpha}\right|\left.\frac{d\left(k\Gamma\right)}{dk}\right|_{k_{\alpha}}\left|\Psi_{\alpha}\right\rangle$ being equal to a constant (which we assign as $z_{\alpha}$), so that a finite volume normalization expression reads:

[TABLE]

where $\mathcal{V}$ is a finite volume surrounding the system and $\mathcal{S}$ its surface. Although not necessary, we used $k=\omega/c_{b}$, to formulate Eq.​​ 34 in terms of exterior medium wavenumber, $k$, (as opposed to angular frequency, $\omega$). Several works in recent years have obtained formulas quite similar to Eq.​​ 34[25, 26, 46], but the formula with the most explicit agreement is Eq.(20) of Ref.[39], derived in the context a six-by-six component Green’s function (except for a differing RS convention - obtained by replacing $z_{\alpha}$ on the right hand side of Eq.​​ 34 by $1$).

The frequency derivative (or $k$ derivative) of fields in the surface integral of Eq.​​ 34 could be deemed undesirable for numerical analysis, but we can take advantage of the scale invariance in the region outside the system to replace $d/dk$ by $\frac{\boldsymbol{r}\cdot\boldsymbol{\nabla}}{k}$, as done in Ref.[39, 46]. Those favoring a surface integral approach to RS normalization often choose to formulate normalization in terms of volume integrals involving only the electric field (in contrast to the ‘Lagrangian’ field formulation of this work). One can transform from one formulation to the other by invoking a Poynting vector analysis which yields:

[TABLE]

which in the case of null magnetic permeability contrast, i.e. $\mu\rightarrow 1$ allows the normalization constraint to be written in terms of an electric field volume integral:

[TABLE]

where we used the replacement $\frac{d}{d\omega}\left[\varepsilon(\omega)\right]_{\omega_{\alpha}}\rightarrow 2\frac{d}{dk^{2}}\left[\varepsilon(k^{2})\right]_{k^{2}_{\alpha}}$. The surface terms in Eq.​​ 36 are radial derivatives, and the first appearance of such a normalization in photonics seems to be have been derived in 2010 for dispersionless media [48].

Even though the normalization formulas of Eqs.​ 32, 34 and 36 may appear distinct from one another, the Gaussian regularization perspective shows their equivalence, under the proper conditions, since they are all determined by the same underlying condition. This equivalence can be demonstrated directly for a spherical surface, $\mathcal{S}$, of radius $R$ by using Hankel function recursion relations to show that the surface integral on the right hand side of Eq.​​ 34 of an arbitrary spherical surface of radius $R$ is analytically identical to the formula of Eq.​​ 102 for a Gaussian regularized volume integration of the electromagnetic field carried out in the region outside this sphere. Such analytic equivalencies can also be extended to arbitrary scatterer geometries by invoking the multipole decomposition of Eq.​​ 8.

2.4.1 Analytic formulas for fully dispersive spherical scatterers

Analytical RS normalization expressions are derived in Appendix F using a Gaussian regularization of the full volume integration of Eq.​​ 32 for an electrically and magnetically dispersive sphere adopting the expressions for the RSs in section 2.1 and using the formulas in Appendix C and Ref.[16], we find:

[TABLE]

For dispersionless media, light scattering becomes scale invariant and Eq.​​ 37 simplifies to expressions which only depend on the size parameter of the mode, $z_{\alpha}$, and the relative constitutive parameters, $\varepsilon$ and $\mu$ (now real-valued) as was previously derived in Ref.[17]:

[TABLE]

Comparisons of the above normalization formulas for spherical scatterers with other analytic formulations are discussed in Appendix A.

It is insightful to plot the radial dependence of the normalization integrals of Eq.​​ 33 (after angular integration) for a lossless dielectric of $\varepsilon=16$ as shown in Fig.​​ 5. The convergence factor was simply set to $\eta=0$ in these graphs since the exponential divergence only becomes relevant for $r/R$ significantly larger than the $r/R<4$ range in the RSs plots for the $z_{2}$, $z_{3}$, and $z_{4}$ RSs. For such low loss RSs, a numerical ‘verification’ of Eq.​​ 33 or Eq.​​ 32 is simple because the exponential divergence only becomes non-negligible after passing through a long region where the integrand is nearly zero. Therefore, truncations of the low loss RS integrals will approximately verify Eq.​​ 33 as one can guess by visual inspection of Figs.5(b)-(d).

The situation is quite different however for poorly confined RSs with large imaginary parts, like $z_{1}\simeq 1.04-i0.5$. As one can see in 5(a), the exponential divergence sets in immediately outside the sphere and a regularization of the RS normalization is essential if we look to numerically reproduce the analytic normalization result of Eq.​​ 33. Similar behavior is observed for the RSs of lossy scatterers studied in section 4.3 below.

Unlike the physics of closed systems, where the phase of the modes can often be treated as arbitrary, the RS ‘normalization’ factors, $\mathcal{N}_{\alpha}$, must correctly adjust the phase of the RSs which is a crucial information for reconstructing response functions. This point will be reinforced in section 3, where we will see that the RS phase determines the complex residues in spectral developments of response functions.

The ultimate justification for setting the right hand sides of Eqs.​ 32, 34 and 36 to $z_{\alpha}$, is that this choice generates particularly simple spectral response functions presented in section 3 below. A more widespread choice is obtained by defining a rescaled RS wavefunction $\widetilde{\Psi}_{\alpha}(\boldsymbol{r})\equiv\Psi_{\alpha}(\boldsymbol{r})/z_{\alpha}^{1/2}$ for which the RS normalization reads:

[TABLE]

which is reminiscent of normalization in closed systems like conventional quantum mechanics.

2.4.2 Mode volumes

Recent derivations of the Purcell factor for RSs define a mode volume, $V_{\alpha}$.[23, 24, 26, 37] These derivations generally invoke a Green’s function formalism and the local density of states (LDOS), which generally leads to a ‘complex’ valued mode volume[49]. A typical definition for $V_{\alpha}$ derived along the LDOS Green function approach is:

[TABLE]

where $\boldsymbol{r}_{c}$ is a ‘chosen’ position (often the field maxima) while $\widetilde{\boldsymbol{\rm f}}_{\alpha}(\boldsymbol{r})$ is one of the commonly employed photonic notations[24] for the RS electric field, and the double ‘$\langle\langle...\rangle\rangle$’ is an inner product notation ref.[24]. Following the usual LDOS route to mode volume in the framework of this work leads to,

[TABLE]

which up to possible proportionality factors agrees with more recent mode volume definitions that continue to be debated even now[50].

Mode volume assignments along the lines of Eqs.​ 40 and LABEL:{lagvol} tell us that $|V_{\alpha}|$ is proportional $1/\left|\mathcal{N}_{\alpha}\right|^{2}$, so we generally expect small values of $\left|\mathcal{N}_{\alpha}\right|$ to correspond to large mode volumes. One indeed sees from Table 2 that poorly confined modes like that $z_{1}$ (cf. Fig.​​ 5(a)) have smaller $\left|\mathcal{N}_{\alpha}\right|$, than the more tightly confined modes illustrated in the other graphics of Fig.​​ 5.

Inspection of Eq.​​ 38 leads to the physically intuitive conclusion that weak constitutive contrasts lead to large mode volumes (i.e. small values of $\left|\mathcal{N}_{\alpha}\right|$). Nevertheless, the $\mathcal{N}_{\alpha}$ factors in Eq.​​ 38 do not vanish in the limit of zero contrast since the vanishing $(\varepsilon-1)$ and $(\mu-1$) factors, are multiplied by factors expressed in terms of spherical Hankel functions evaluated at the RS frequencies. The fact that the $\mathcal{N}_{\alpha}$ normalization factors remain finite in the limit of null contrast is in fact a necessary ingredient for retrieving a vanishing $T$-matrix response, (cf. section 3), in the zero contrast limit. A detailed analysis of this intriguing result will be addressed elsewhere.

3 Resonant state expansions of Mie theory response functions

One can obtain physical response functions directly from the Green’s function, but we have argued in recent years the advantages of employing elegant and powerful scattering formalisms involving $S$ and $T$-matrices and reaction matrices (cf. [17, 18, 22, 35, 36, 37, 51] for more details and applications). An essential idea is seen by remarking that the $T$-matrix operator, $\mathbb{T}$, is related to the Green’s function as:[52, 53, 54, 55, 56]

[TABLE]

where $\mathbb{G}_{0}$ is the Green function of the homogeneous background. An advantage of this formalism is that the $\mathbb{T}$-operator describes electromagnetic response in terms of fields (as opposed to source current for Green’s functions). This occurs because the $\mathbb{G}_{0}$ operators on the left and right side of $\mathbb{T}$ in Eq.​​ 42 respectively generate outgoing and incident multipole fields, (the normalization in Eq.​​ 32 accounts for differences in the fields generated by $\mathbb{G}_{0}$ operators and the conventional multipole fields).[52, 43, 55] Once projected onto the multipole basis, the $\mathbb{T}$-‘operator’ takes the form of a ‘matrix’ on the multipole basis functions (which is precisely the justification of the terminology $T$-matrix).

Given the observations made around Eq.​​ 7 and in the previous paragraph, we conclude that the RS residues of the scattering $T$-matrix are $\mathcal{R}_{\alpha}\equiv i/\mathcal{N}^{2}_{\alpha}$, as will be confirmed below for spherical scatterers, since the $T$-matrices in this case are diagonal on the multipole basis (with matrix elements given analytically by Mie theory for complex frequencies). This agreement is remarkable in that $\mathcal{N}_{\alpha}$ was obtained through RS regularization, quite outside the context of traditional Mie theory. Although it is clear from Eq.​​ 42 that $T$-matrices are just a reformulation of the Green function of Eq.​​ 7, we have found in our prior works that the $T$-matrix formulation has the advantage of facilitating the derivation of the non-resonant contributions to response functions [17, 18, 22, 35, 36, 37, 51] (notably through considerations of causality and energy conservation [21]). Relevant results of this approach, and some new derivations in the context of RS Mie theory, are reviewed below.

3.1 RS spectral expansions

In the interest of completeness, essential elements of Mie theory and linear response theory are recalled in Appendix D. Since Mie theory predates the theory of matrix response functions, it remains common in the literature to keep the traditional notations [42] for the diagonal ‘$T$-matrix’ elements; $T_{h,n}\equiv-b_{n}$ and $T_{e,n}\equiv-a_{n}$, and $\Omega_{h,n}=c_{n}$ and $\Omega_{e,n}=d_{n}$ for the internal field coefficients.

The meromorphic expansion of the $T_{n}$ functions in terms of frequency, $\omega$, have been determined to be:[17, 18, 22, 51]

[TABLE]

The residues, $r_{\alpha}$, of the RSs in the RS expansion of the $T$-matrix in Eq.​​ 43a are,

[TABLE]

where the $e^{2iz_{\alpha}}$ phase factor doesn’t affect the $T$-matrix residue Eq.​​ 43a since it is canceled by the $e^{-2iz}$ multiplicative factor so the RS residues are indeed $\mathcal{R}_{\alpha}=i/\mathcal{N}^{2}$ as expected.

The internal field coefficients of Mie thery, $\Omega_{q,n}$, are quite similar $T$-matrix expressions with the exception of their response being expressed entirely in terms of RS states, (i.e. there are no non-resonant contributions):

[TABLE]

As one could deduce from Eq.​​ 7 and the resonant state expression in Eq.​​ 16, the residues of internal RS coefficients, are related to the $T$-matrix residues in Eq.​​ 45, as $r_{\alpha}^{(\Omega)}=\gamma_{\alpha}r_{\alpha}$. We will see below that this relation can also be deduced directly from Mie theory.

3.2 Spectral residues in Mie theory

Given the expressions for $\Omega_{n}$, and $T_{n}$ in Eqs.​ 43b and 45 in terms of the $N_{q,n}$ and $D_{q,n}$ functions defined in Eqs.​ 94 and 95 of section D.2, the residues in Eqs.​ 45 and 43b can alternatively be calculated from Mie theory, as:

[TABLE]

One should note that for non-dispersive media, the scale invariance of the underlying electromagnetic equations imposes that the residues, $r_{\alpha}$ and $r_{\alpha}^{(\Omega)}$, only depend on particle radius, $R$, as a function of, $z_{\alpha}=\omega_{\alpha}R$, which results in the spectral response functions of Eq.​​ 45 and Eq.​​ 43 becoming solutions for all particle sizes and frequencies in this case.

We remark from the second equality in Eq.​​ 46b that the relation between the $T$-matrix residues, $r_{n,\alpha}$, and the internal field residues, $r_{n,\alpha}^{(\Omega)}$, is purely analytic. This relation can be rewritten in a more transparent form by invoking the respective RS conditions of Eq.​​ 14 such that the expressions of $N_{e,n}$ and $N_{h,n}$ can be rewritten:

[TABLE]

for the magnetic modes. In both Eqs.​ 47a and 47b we invoked the Wronskian relation,

[TABLE]

and recalled the definitions of the $\gamma_{\alpha}$ functions defined in Eq.​​ 14.

Mie theory therefore analytically validates the relation

[TABLE]

that we deduced after Eq.​​ 45 using arguments based on RS spectral expansions of the response function. Although the derivation is a bit lengthy here, one can also directly derive the $T$-matrix residue, $\mathcal{R}_{\alpha}=i/\mathcal{N}^{2}$, directly from Mie theory. We found this fact quite striking since Mie theory makes no use whatsoever of the RS regularization methodology. We interpret this as striking evidence for the mathematical validity of the RS regularization program.

4 Numerical implementations

This section provides precise calculations of RS eigenvalues, normalization, and response function reconstructions in the case of spherical scatterers. As discussed in section 2, generalizations to non-spherical shapes are quite analogous, but they will not be considered here since they are not amenable to such high precision analytics.

Even for spherical geometries, solving for the RS frequencies must be carried out numerically, but thanks to the field continuity conditions of Eq.​​ 14, the RS frequencies are determined as solutions to analytic expressions that can be expressed in terms of the reduced logarithmic derivatives of the Ricatti-Bessel functions:

[TABLE]

where we recall from Eq.​​ 15 that $\varepsilon_{\alpha}$, $\mu_{\alpha}$, and $\rho_{\alpha}$ are in general frequency dependent constitutive parameters evaluated at the RS frequency, $\omega_{\alpha}$, such that the first terms in Eq.​​ 50 have a non-trivial frequency dependence.

4.1 Resonant states and normalizations of high refractive index dielectric spheres

There is currently considerable interest in non-dispersive scatterers since high index dielectrics are being considered in nano-optics as an alternative to plasmonics [57]. Even for such energy conserving systems, one still needs to solve Eq.​​ 50 numerically, but the difficulty is greatly reduced and the problem becomes scale invariant, which means that one solves the problem for all particle sizes and/or frequencies simultaneously. Furthermore, the RS normalization coefficients are calculated with the simple expressions found in Eq.​​ 38.

Numerical results for a few low-order modes of a dielectric with no magnetic permeability contrast, $\mu=1$ and high, $\varepsilon=16$, permittivity contrast are shown in Table 1 (further analysis and discussion of such modes may be also found in [18]). There are an infinite number of modes for each value of $n$, and they follow an asymptotic pattern resembling that of zeros of Bessel functions. Beginning with electric modes, table 1 gives the values of $z_{\alpha}$ and the normalization factors for the first two dipole and quadrupole modes. An interesting feature of the RSs in this case is that for $n$ even (odd), there is one electric (magnetic) mode with purely imaginary $z$, labeled with $\ell=0$. One can also remark from table 1 that for imaginary RS eigenvalues, the interaction ‘strengths’ (residues), $\mathcal{R}_{\alpha}$, have purely imaginary values.

The high number of significant figures given in Table 1 for the RS eigenvalue size parameters, $z_{\alpha}$ and normalizations, $\mathcal{R}_{\alpha}=i/\mathcal{N}^{2}_{\alpha}$ were adopted because these values were indeed calculated up to such high order accuracy which underscores the results of this work to serve as benchmarks for more numerically based approaches.

4.2 Numerical verification of RS expansions of Mie theory.

The steps to implement the results of this paper up to the construction of Mie theory response functions can be summarized as follows :

For each required multipole order $n$, and mode type $q=(0,1)$, (i.e. $(h)$ and $(e)$ respectively), the only non-trivial step is to find a sufficient number of the resonant state size parameters, $z_{\alpha}$, by numerically solving the transcendental Eq.​​ 14 (or equivalently Eq.​​ 50). The sufficient number of RSs depends on the particle radius $R$ (or frequency of interest) and the dispersion relations of the constitutive parameters, but can be as low as 1 for some applications in dispersive materials. When only a few RSs are required, graphical solutions with numerical refinements may suffice (see section 4.3 for examples, but one should exercise caution with such techniques since important RSs can occur far from the real $z$ axis as can be seen in Fig.​​ 6(b) and Fig.​​ 6(e)). 2. 2.

Calculate the square of the complex-valued RS norms, $\mathcal{N}_{\alpha(q,n,\ell)}^{2}$, using Eq.​​ 37. 3. 3.

Insert the values of $\mathcal{N}_{\alpha(q,n,\ell)}^{2}$ found in step 2 into Eq.​​ 44 to determine the residues, $r_{\alpha(q,n,\ell)}$, for the meromorphic expression of Eq.​​ 43b for the scattered field Mie coefficients, $a_{n}=-T_{e,n}$ and $b_{n}=-T_{h,n}$. Cross section contributions, $\sigma_{n}$, can then be determined by the standard formulas recalled in Eq.​​ 96 of section D.2. 4. 4.

If one needs to determine the fields inside the sphere, use the values of $r_{\alpha}$ found in step 3 and the $\gamma_{\alpha(q,n,\ell)}$ coefficients of Eq.​​ 14 to calculate the residues, $r_{\alpha(q,n,\ell)}^{\left(\Omega\right)}$, of Eq.​​ 49 that appear in the meromorphic expression of Eq.​​ 45, for the internal field Mie coefficients $c_{n}=\Omega_{h,n}$ and $d_{n}=\Omega_{e,n}$ (defined in Eq.​​ 93b).

The results displayed in Table 1 of section 4.1 were obtained by applying steps 1 and 2 of the above steps. An illustration of carrying out step 3 above, with the same material parameters ($\varepsilon=16$, $\mu=1$) used in Fig.​​ 5, is illustrated graphically for both electric and magnetic mode contributions and orbital quantum numbers, $n=1$ and $n=4$. The first few RSs for each mode are given as blue dots in the lower half plane with ${\rm Im}(z)$ on the vertical axis, while the associated contributions to the cross section efficiencies, $Q\equiv\sigma/(\pi R^{2})$, are plotted on the positive vertical scale. In Fig.​​ 6(c) and Fig.​​ 6(f), we verified that this procedure reproduces Mie theory up to arbitrary computational accuracy for all multipole orders (the curves in Fig.​​ 5 were generated both from RS calculations and Mie theory, but are identical since they agreed here up to six figure accuracy).

Similar calculations can reproduce Mie theory for dispersive media, but in that case, the problem is no longer scale invariant and the position of the RS eigenfrequencies will depend on the particle size and will be seen in section 4.3 below.

4.3 Resonant states and normalizations of ‘Drude’ model conductors

The calculation of the complex wavenumbers of RS’s depends on an extrapolation of the experimental data on complex permittivities (or permeabilities in cases involving magnetic materials) from the real axis of frequencies or wavelengths. Such an extrapolation depends critically on both the accuracy of the experimental data and on the accuracy of the analytic model used for the extrapolation.

Data on the complex refractive index of gold as a function of frequency are available in Volume 1 of he comprehensive collection due to Palik [58]. The simplest widely used analytic expression for this data is that of the Drude model, which may be written in terms of the angular frequency, $\omega$, for gold [23] and silver [59] as:

[TABLE]

The comparison with experimental data given in Fig.​​ 11 of Appendix G shows that the Drude model works relatively well for silver at wavelengths longer than around $0.30\mu{\rm m}$, while the comparison of the Drude model for gold in Fig.​​ 12 shows that the Drude model for gold only works well for wavelengths longer than $0.60\mu{\rm m}$.

To achieve a better fit for gold valid down to shorter wavelengths, extra Lorentz resonant terms need to be added to the Drude model of Eq.​​ 51b. For example, Sikdar and Kornyshev [60] give the parameters for a model for gold taking into account inter-band transitions via two additional resonances:

[TABLE]

Here the parameters are given in eV: $\epsilon_{\infty}=5.9752$, $\hslash\omega_{pD}=8.8667{\rm eV}$, $\hslash\Gamma_{D}=0.03799{\rm eV}$, $s_{1}=1.76$, $\hslash\omega_{p1,L}=3.6{\rm eV}$, $\hslash\Gamma_{1,L}=1.3{\rm eV}$, $s_{2}=0.952$, $\hslash\omega_{p2,L}=2.8{\rm eV}$ and $\hslash\Gamma_{2,L}=0.737{\rm eV}$. The conversion of $\omega$ to electron volts is achieved by dividing the frequency results by $e/\hbar=1.51927\times 10^{15}$. As can be seen from Fig.​​ 12, this model works well down to around $0.40\mu{\rm m}$. An interesting feature of the dispersion model of Eq.​​ 52 is that it supports bulk plasmons. [61] These occur for complex wavelengths or frequencies at which $\varepsilon(\omega)=0$. For the dispersion model of Eq.​​ 52 they are at: $\lambda_{\rm Lg}=0.257778+i0.0206709\mu{\rm m}$, $\lambda_{\rm Lg}=0.395618+i0.0536046\mu{\rm m}$ , $\lambda_{\rm Lg}=0.502518+i0.048786\mu{\rm m}$. Bulk plasmons are longitudinal waves, which do not couple to transverse electromagnetic waves and cannot be excited by or scattered into them. This translates as the fact that one obtains a normalization residue factor, $\mathcal{R}_{n,\alpha}$, of zero for these longitudinal ‘resonant states’.

Fig.​​ 7 shows the modulus of the dispersion equation of Eq.​​ 50a for $n=1$ and for a gold sphere of radius $100{\rm nm}$, both for the Drude model of Eq.​​ 51b and for the more elaborate model of Eq.​​ 52. The former has its minimum at $\lambda_{\alpha}=0.606976+i0.239112\mu{\rm m}$, in good agreement with the value reported in Sauvan et al. [23]. The latter has a slightly different value: $\lambda_{\alpha}=0.592227+i0.210097\mu{\rm m}$. From Table 2 the values of $\varepsilon_{\alpha}$ at the resonance and consequently the normalization factor disagree more significantly for the more elaborate model.

Fig.​​ 8 shows the inverse modulus of the $n=1$ dispersion relation of Eq.​​ 50 for a gold sphere of radius $80{\rm nm}$. The plot is more complicated for the more elaborate dispersion model of Eq.​​ 52 since the decrease in radius moves the wavelength range of interest into that for which the zeros, poles and branch cuts of this dispersion model become evident (see the data listed above for the wavelengths corresponding to bulk plasmons). The complex wavelength of the RS for the Drude model is $\lambda_{\alpha}=0.505163+i0.174433\mu{\rm m}$. With the model of Eq.​​ 52, this becomes $\lambda_{\alpha}=0.547815+i0.080062\mu{\rm m}$. Note the significant decrease in the imaginary part around this wavelength, as seen in Fig.​​ 12 shows, since the Drude model is only a crude approximation to the experimental data.

The aforementioned electric dipole RSs for the different dispersion models are given in Table 2 for RSs in silver and gold spheres of radius 100 and 80nm. In order to calculate the mode normalization factors, numerical differentiation, as described in Eq.​​ 46b, can be used to find the necessary dispersion derivatives. The data given in Table 2 for the resonances of silver spheres of radius 100 and 80nm shows complex wavelengths and normalization factors which are broadly similar to the corresponding values for gold.

5 Conclusions and perspectives

We used Gaussian regularization to demonstrate that resonant states can be assigned finite normalization and orthogonality properties that are closely analogous to definitions adopted for the eigenstates of closed systems. A notable difference however is that the RS normalization factors, $\mathcal{N}_{\alpha}$, are complex valued, (in contrast to their real-valued counterparts in closed systems). Although, “irrelevant” phase factors are a familiar feature of the time evolution of eigenstates in dispersionless closed systems, phase plays a crucial physical role in an open system’s response to an external excitation (and consequently can no longer be treated as arbitrary). The RS normalization factor must therefore ensure the physically correct RS phase, in addition to its amplitude. Recognizing the importance of phase factors from the outset helps one to understand the fact that when RS physics is imposed on closed system concepts, like the Purcell factor, the complex normalization factor inevitably leads to the introduction of apparently incongruous notions like complex valued mode volumes.[50, 49, 24]

For systems of arbitrary geometry, the determination of the RS eigenvalues and the regularization of their inner products are often carried out numerically. This work combined a multipolar expansion of arbitrary RSs with Gaussian regularization in order to analytically evaluate the inner product integrals in the region lying outside the system of interest. Our analytical analysis can both serve as benchmarks for purely numerical treatments, as shown in section 4, and also increase the speed of numerical approaches by replacing certain numerical integrals with analytical formulas. The results of this work for RS normalization have been shown to be consistent with other formulas in the literature for non-spherical scatterers and we have explicitly generalized previous analytical results of the literature for dispersive spherical scatterers with both permittivity and permeability contrasts.

The analytic formulas obtained herein enabled us to completely reformulate the ‘Mie theory’ response functions of dispersive spherical objects in terms of resonant state expansions, as summarized in section 3, and illustrated in Fig.​​ 5. A notable by-product of the response function formulation was that it allowed an analytical determination of the controversial ‘non-resonant’ contributions that have repeatedly appeared as snags in RS reconstructions of response functions. Given the vast popularity and practical utility of Mie theory in applications, there can be little doubt that the ability to reformulate Mie theory in terms of the conceptually and mathematically simple RS expansions will find future, and likely unforeseen, implementations.

6 Acknowledgements

Research conducted within the context of the International Associated Laboratory for Photonics between France and Australia (LIA ALPhA). This work has been carried out thanks to the support of the A*MIDEX project (no. ANR-11-IDEX-0001-02) funded by the Investissements d’Avenir French Government program, managed by the French National Research Agency (ANR). We would like to thank Jean-Paul Hugonin for calling our attention to the numerical equivalence between the Gaussian normalization and the PML deformation into the complex plane as demonstrated in Figure 2. We would also thank the referees whose remarks helped improve this manuscript, most notably concerning notation and more thorough the comparisons with results in literature.

Appendix A Alternative normalizations and literature comparisons

Given the wide variety of constantly evolving notations and conventions in the QNM/RS literature, it is not realistic to carry out exhaustive comparisons, so in this appendix we show that our normalization formulas agree with other analytic formulations in the special case of purely dielectric scatterers that have previously been given in the literature, notably those by Muljarov et al. and in particular in Refs.[26, 46].

Although we adopted a six component electromagnetic field formulation in contrast to the more common formulations entirely in terms of the 3-component electric field, such differences should not affect the values of RS normalization. One remarks that Ref.[46] labeled RS indices ‘$n$’ as opposed to our generalized index $\alpha$, and they used an index, $l$, for multipole number as opposed to our ‘Mie theory’ choice of index $n$. Comparisons between our formulas and those given in Refs.[46, 26] are facilitated by introducing the definitions of Ref.[46]. These authors used radial dependent constitutive parameters of a spherical scatterer are defined as:

[TABLE]

and defined a function $R_{n}$:

[TABLE]

together with

[TABLE]

where the various spherical Bessel functions, $j_{n},\psi_{n},h_{n}$, and $\xi_{n}$ are defined here in Appendix B.

Invoking our definitions of Eq.​​ 14 together with the RS relations for spherical scatterers in Eq.​​ 50, our six component RSs Eq.​​ 16 can finally be compactly reexpressed:

[TABLE]

where the final set of parentheses are identical to the parenthesis of the central expressions; thus establishing the connection between the normalization coefficients $A^{\rm TE}_{\alpha}$ and $A^{\rm TM}_{\alpha}$ of Ref.[46], with our normalization factors, $\mathcal{N}_{\alpha}$. Note that we normalized the coefficients of Ref.[46] $\overline{A}_{\alpha}\equiv A_{\alpha}\sqrt{n(n+1)}$, in order to compensate for the $1/\sqrt{n(n+1)}$ normalization factor of spherical harmonics that was included in our vector spherical harmonics $\boldsymbol{X}_{q,n,m}$, $\boldsymbol{Y}_{q,n,m}$, and $\boldsymbol{Z}_{q,n,m}$, given in Appendix B (also note that $\sqrt{n(n+1)}$ reads $\sqrt{l(l+1)}$ in Ref.[46]).

The electric fields in the parenthesis of Eq.​​ 65 can be directly identified with the RS electric fields in Eq.(26) and Eq.(27), of Ref.[46] (again up to the $1/\sqrt{n(n+1)}$ factor). In order to make our RS fields of Eq.​​ 16 match the RS field expressions of Ref.[46], we had to divide our multipolar fields $\Phi_{\alpha}$ by $h_{n}(z_{\alpha})$, which we then had to compensate for via a modified normalization factor, $\widetilde{\mathcal{N}}_{\alpha}$:

[TABLE]

which results in particularly compact and symmetric expressions:

[TABLE]

where $\varphi_{n}^{(+)}(z)$ is the reduced logarithmic derivative of the outgoing Riccati Hankel function, $\xi_{n}(z)$:

[TABLE]

Finally, we remark that we had to multiply the $A_{\alpha}$ coefficients of Ref.[46] by an additional $z_{\alpha}$ factor in order to account for the fact that we normalized the RS to a $z_{\alpha}$ in Eq.​​ 32 as opposed to a normalization of the RSs to unity. Finally, reading off the relation between the $\overline{A}_{\alpha}$ and $\widetilde{\mathcal{N}}_{\alpha},$ in Eq.​​ 65 we find our expressions of $\widetilde{\mathcal{N}}_{\alpha}$ factors provided we first set $\mu_{\alpha}(\omega,r)=1$ and $\varepsilon_{\alpha}\rightarrow\varepsilon=\rho^{2}$ with $\varepsilon$ and $\rho$ real valued:

[TABLE]

The agreement with the $A_{\alpha}^{\mathrm{TE}}$ factor with our formulas of Ref.[46] is immediate by remaking that $\varepsilon\rightarrow n_{R}^{2}$ where $n_{R}$ is the notation of Ref.[46] for relative index contrast. Our expression in Eq.​​ 69b for $\overline{A}_{\alpha}^{\mathrm{TM}}$ also agrees with that of Ref.[46] even though this is not immediately clear by inspection. In order reveal this agreement, we first recall that the magnetic mode RS condition of Eq.​​ 50b can be expressed:

[TABLE]

which allows us to write,

[TABLE]

and if we express this in a manner of Ref.[46], Eq.​​ 71 becomes,

[TABLE]

which in their notation reads,

[TABLE]

so our expressions finally agree exactly provided we take $\mu(\omega,r)=1$ and set $n_{R}\rightarrow\sqrt{\varepsilon}=\rho$.

In the above, we only compared the formulas for dispersionless media normalization of Ref.[46], but proceeding along the same lines permits the retrieval of the comparison with later formulas derived for dispersive media [26].

Appendix B Vector partial waves and vector spherical harmonics

The regular Vector Partial Waves (VPWs) are defined as:

[TABLE]

while for the outgoing $(+)$ and incoming ($-$) waves,

[TABLE]

In Eqs.​ 74 and 75, we used the Riccati-Bessel functions $\psi_{n}$ and $\xi_{\pm,n}$ that are functions of a complex variable $z$, defined by,

[TABLE]

and the prime is the derivative with respect to the argument, i.e.,

[TABLE]

The angle dependent vector spherical harmonics, $\boldsymbol{X}$, $\boldsymbol{Y}$, and $\boldsymbol{Z}$ of the ‘electric’ mode type ($q=1$) are defined for $m=0,1,...,n:$

[TABLE]

while the magnetic VSHs are defined for $m=1,...,n$,

[TABLE]

We remark that in the classic text of ref. [42], the magnetic and electric modes were referred to as ‘odd’ and ‘even’ modes respectively, but we found the terminology magnetic and electric to have more physical relevance. We also point a recent convention of associating the magnetic (‘odd’) modes with negative $m$ and the electric (‘even’) modes with non-negative $m$ values. Although such a notation has the advantage of eliminating the need of a separate index to distinguish electric and magnetic type modes, we found it useful to keep the separate index here for reasons of clarity in explanations and notation.

The normalized associated Legendre functions are defined,

[TABLE]

with a normalization factor,

[TABLE]

and we adopted a common definition of (positive $m$) associated Legendre functions as:

[TABLE]

Appendix C Killing Mie softly

This section demonstrates how an integral of special functions with diverging amplitudes can yield finite results by first taming the relevant integrals by a Gaussian factor, $\exp(-\eta x^{2})$ and then taking of the limit $\eta\rightarrow 0$. The analytic arguments presented here and in [16] may also be of help to those developing purely numerical methods for nanophotonics. As described in [16], if the integrand oscillates about a non-zero mean value, this needs to be treated separately, as it may generate a delta function contribution (such contributions cancel out in the cases treated here). The oscillations about a mean of zero will in general give a zero contribution for the integration variable tending to infinity, as in the spirit of generalized function theory.

We now give an example of the results in the references [16, 62]. The integral chosen is that over the product of the spherical Bessel functions $j_{n}(Kx)$ and $y_{n}(kx)$, where $K$ and $k$ are (possibly complex) wavenumbers: see equation LABEL:snr18.

[TABLE]

The analytic result for this integral from [62] is given by the rightmost expression in (LABEL:snr18). Here $h_{-1,b}$ denotes an associated Bessel function [63], and ${\cal H}(b,k,K,\eta)$ is the following finite-range integral:

[TABLE]

Using expansions given by Luke [63] and evaluating the finite integral in (84) by direct numerical integration, the expression (LABEL:snr18) may be verified for arbitrary choices of the parameters $K$, $k$ and $\eta$.

The asymptotic treatment given in [62] which takes the limit as $\eta\rightarrow 0$ is lengthy and complicated. However, it yields a simple result, which is in keeping with the well-known result from Watson [44]:

[TABLE]

where ${\cal C}_{\mu}$ and ${\cal D}_{\mu}$ are cylinder functions of integer or real order $\mu$. For the integral (LABEL:snr18) we need to take ${\cal C}_{\mu}=J_{n+1/2}$ and ${\cal D}_{\mu}=Y_{n+1/2}$. The expression (85) is used (with a minus sign) to give the contribution from the lower limit to the integral (LABEL:snr18). In the case chosen, the asymptotics show that the contribution from the upper limit (infinity) is precisely zero. (For other cases, like the integrals ${\cal I}_{jj}$ and ${\cal I}_{yy}$, there is a delta function contribution which comes from the upper limit; this arises if the asymptotic expansion of the integrand in expressions like (LABEL:snr18) with $\eta=0$ as $x\rightarrow\infty$ contain a constant term in addition to terms oscillating around zero.) The result is in this case:

[TABLE]

We now give an example of the effectiveness of the Gaussian “killing” function technique for the integral (LABEL:snr18): see Fig.​​ 9. Even with a Gaussian with $\eta$ only equal to 0.01, the divergent integrand is replaced by one which can be integrated accurately. The numerical integration of the Gaussian form with $\eta=0.01$ gives $0.0164787-0.0138487i$, for integration with upper limit 80 or beyond. The analytic value for the integral (86) is $0.0163332-0.0135188i$. For $\eta=0.005$, the numerical integral gives $0.0164062-0.0136812i$, slightly closer to the exact answer, while for $\eta=0.001$, the numerical integration in Mathematica fails.

As a second example, we give in Fig.​​ 10 plots of the integrand in the following normalization integral:

[TABLE]

As $\eta$ gets smaller, the oscillating real and imaginary parts of the integrand become for larger $x$ concentrated between Gaussian envelopes. The envelopes peak round $x=-\Im(k)/\eta$, and have a $1/e$ full width of $2/\sqrt{\eta}$. The peak modulus of the envelope is $\exp\{[\Im(k)]^{2}/\eta\}$, and the oscillations within the envelope go as $\exp[2i\Re(k)x]$. The oscillatory behavior means that the value of the integral over this peaked region alternates between positive and negative values, and has on average the value zero. To carry out the integral directly for small $\eta$ becomes increasingly difficult, particularly if $\Re(k)$ is small in magnitude.

The peak value of the envelope increases rapidly as $\eta$ decreases, while there are more and more cancelling positive and negative contributions to the integral. In order to obtain an accurate numerical estimate for the integral, more and more decimal places have to be used in evaluations of the elements of the integrand, if the oscillations are to cancel the diverging peak value.

This means, for calculations with a fixed accuracy (say 10-12 decimal places) there will be a smallest value of $\eta$ for which the numerical estimate is good. The calculations will get more and more difficult as the ratio $|\Im(k)/\Re(k)|$ increases. These remarks emphasize the value of the analytic studies [16] which have shown there is a simple explicit and general answer for normalization integrals such as (87).

In region outside the sphere, the specific integrals of the more general exposition of Ref. [16] that are needed for this work are:

[TABLE]

while for $\beta=\alpha$ one finds:

[TABLE]

where $z_{\alpha}\equiv k_{\alpha}R$, with $h_{n}(z)$ the outgoing spherical Hankel functions, and $\xi_{n}(z)\equiv zh_{n}(z)$ the outgoing Ricatti-Hankel functions. We recall that the above multipole integrals are the only one for a finite sized system at these analytically regularized results are the only

The integrals of the regular partial waves inside a dispersionless homogeneous sphere of arbitrary radius $R$ inside the particle can also be carried out analytically without the need for regularization:

[TABLE]

while when $\beta=\alpha$, one finds:

[TABLE]

These formulas are particularly useful for spherical scatterers but could also prove useful in non-spherical particles and multiple scattering situations.

Appendix D Mie theory

D.1 Multipolar decomposition of fields

Mie theory relies on developing fields in different regions in terms of vector partial waves (VPWs), which are homogeneous media electromagnetic wave solutions in an angular momentum eigenstate basis. The excitation field, $\boldsymbol{E}_{\mathrm{e}}$, is the field created by source currents in the absence of the scatterer and it can be developed on the complete basis of ‘regular’ multipolar fields denoted by a superscripted $(1)$,

[TABLE]

The total field inside a spherical particle can also be developed in terms of regular VPWs provided that position vector in the VPWs is weighted by the field frequency dependent wavenumber inside the sphere, $k\rho_{\omega}$ using the notation of Eq.​​ 11,

[TABLE]

Mie theory also involves the ‘scattered field’ which has its traditional definition as the incident field subtracted from the total field in presence of the scatterer, which can be developed in terms of outgoing VPWs,

[TABLE]

which can be viewed as replacing the spherical Bessel functions in the definitions of the regular VPWs by outgoing Hankel functions as shown in Eq.​​ 75.

D.2 Mie coefficients and cross sections

A consequence of linear response of any scatterer is that if one expresses the scattered and internal fields coefficients, $f$ and $s$, as infinite dimensional column matrices, then the scattering and internal field coefficients can be expressed in terms of the column matrix of excitation field coefficients via infinite matrices, $T$ and $\Omega$, such that, $f=T.e$ and $s=\Omega.e$. The fact that scattering by spheres does not change the angular momentum, both $T$ and $\Omega$ become diagonal matrices for spherically symmetric scatterers and traditional ‘Mie’ theory provides algebraic expressions for these diagonal elements, henceforth denoted $T_{q,n}$ and $\Omega_{q,n}$. For the different multipole orders, one has,

[TABLE]

Mie theory also provides expressions for the lesser known $c_{n}$ and $d_{n}$ coefficients between the internal field and the incident field decomposition which can be directly expressed in terms of denominator functions,

[TABLE]

The ‘numerator’ and ‘denominator’ functions in section D.2 are functions of the angular frequency, $\omega$ and particle radius, $R$:

[TABLE]

and,

[TABLE]

where we have again used $z=kR=\omega R/c$ and used the condensed notation of Eq.​​ 11 for $\varepsilon_{\omega}$, $\mu_{\omega}$, and $\rho_{\omega}$, for the (possibly frequency dependent) constitutive parameters: An inspection of Eqs.​ 94 and 95 readily shows for non-dispersive media, the, $N_{q,n}$ and $D_{q,n}$ reduce to functions of size parameter, $z=kR$, only.

Cross sections like those plotted in Fig.​​ 6 are obtained from the Mie coefficients of Eq.​​ 93a via well-known relations:

[TABLE]

Appendix E Resonant state electric field orthogonalization

In the main text, we showed the orthogonality of electric mode RS product integrals for a spherical scatterer when integrating the magnetic field product of the RSs products over all space when $\alpha\neq\beta$. Here we show that RS orthogonality also holds for products of the electric fields of the electric type modes (n.b. the same integrals also arise in magnetic field integration of magnetic type modes). The volume inner product integral inside the sphere for two RSs $\alpha$ and $\beta$ is zero unless they share the same $q,n,m$ numbers, and in the case where they do share the same $q,n,m$, their product integral inside the sphere is,

[TABLE]

One finishes the electric field inner product by carrying out the volume integral in the region outside the scatterer (i.e. for $r>R$),

[TABLE]

where we used Eq.(111a) in ref. [16] to evaluate the integral in the first line of this equation. Simply adding Eqs.​ 97c and 97e then immediately leads to the orthogonality of the electric type RSs in lossless media when integrating their electric fields over all space,

[TABLE]

All the calculations for the magnetic modes when $\alpha\neq\beta$ are completely analogous to those of this section and those of section 2.3.1 in the main text up to overall sign factors (which we did not explicit here for sake of clarity).

Appendix F Resonant state normalization for Mie theory: dispersive media

F.1 Magnetic mode normalization

The normalization for magnetic (TE) mode fields is found through,

[TABLE]

For clarity, the overall magnetic multipole mode sign factor, $(-1)^{n}$, resulting from the $i^{n}$ phase factor in multipole fields (cf. Eqs.​ 9 and 12) is omitted in this section since it has no effect on normalization and is explicitly included in the response functions expansions of Appendix D.

Denoting by, $\varepsilon(\omega)$, the permittivity of the sphere normalized with respect to exterior medium, the contribution to the electric field integral from the volume inside the sphere is,

[TABLE]

Denoting by $\mu(\omega)$, the magnetic permeability of the sphere divided by permeability of the external medium, the magnetic field integral of Eq.​​ 99 is,

[TABLE]

where the $\Xi^{(+)}_{h,n}$ function was again defined in Eq.​​ 37c.

The exterior integrals are again finite despite their divergent kernels following our regularization approach,

[TABLE]

and the sum of the above two integrals yields:

[TABLE]

Finally, putting together the results of Eqs.​ 100a, 100b and 102 into Eq.​​ 99, one finds the normalization factor of Eq.​​ 37a for magnetic type resonant states for spherical scatterers with full temporal dispersion.

F.2 Electric mode normalization

The normalization condition for electric (TM) modes is strictly analogous to those of the magnetic (TE) modes above with:

[TABLE]

where the overall electric multipole mode sign factor, $(-1)^{n-1}$, will be henceforth suppressed for simplicity like we did above for the magnetic modes above.

The $E$-field integration inside the sphere is,

[TABLE]

The $H$-field integral inside the sphere is,

[TABLE]

where $\Xi_{e,n}^{(-)}$ is defined in Eq.​​ 37d.

The integrals outside the sphere are identical to those of Eq.​​ 101 with the roles of electric and magnetic fields reversed, so we obtain the same result as Eq.​​ 102 for the field integrals in a region exterior to a sphere of radius $R$. Inserting the results of Eqs.​ 104a, 104b and 102 into Eq.​​ 103, one finds the normalization factor of Eq.​​ 37b for electric type resonant states for spherical scatterers with full temporal dispersion.

Appendix G Drude models for fitting dispersion relations

Drude models and its extensions can be used to fit experimental measurements like those of Johnson and Christy for silver and gold as shown respectively in Fig.​​ 11, and Fig.​​ 12.