Spectral expansions of open and dispersive optical systems: Gaussian regularization and convergence

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, ω_α . Specifically, let us consider the generalizable example of a single scattering particle as shown in figure 1. In the homogeneous background media lying outside a sphere of radius R_out, an RS field, Ψ_α ( r ), can be expanded in terms of multipolar basis functions, ${{\Phi}}_{q,n,m}^{(+)}(\omega ,\boldsymbol{r})$:

$$\begin{equation}\langle \boldsymbol{r}\vert {{\Psi}}_{\alpha }\rangle \equiv {{\Psi}}_{\alpha }(\boldsymbol{r})=\sum\limits _{q=0}^{1}\sum\limits _{n=1}^{\infty }\sum\limits _{m=0}^{n}{c}_{q,n,m}^{(\alpha )}{{\Phi}}_{q,n,m}^{(+)}({\omega }_{\alpha },\boldsymbol{r})\quad r{ >}{R}_{\text{out}},\end{equation} \tag{ 8 }$$

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

Zoom In Zoom Out Reset image size

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 ⩾ 1, and the azimuthal (projection) numbers m < n. Using multipole analysis, the ${{\Phi}}_{q,n,m}^{(+)}$ functions can be defined as:

$$\begin{equation}\begin{aligned}{{\Phi}}_{h,n,m}^{(+)}(\omega ,\boldsymbol{r})& \equiv {{\Phi}}_{0,n,m}^{(+)}(\omega ,\boldsymbol{r})\equiv {k}^{3/2}{\text{i}}^{n}\left(\begin{matrix}\hfill {\boldsymbol{M}}_{h,n,m}^{(+)}(k\boldsymbol{r})\hfill \\ \hfill -\text{i}{\boldsymbol{N}}_{h,n,m}^{(+)}(k\boldsymbol{r})\hfill \end{matrix}\right)\quad \\ {{\Phi}}_{e,n,m}^{(+)}(\omega ,\boldsymbol{r})& \equiv {{\Phi}}_{1,n,m}^{(+)}(\omega ,\boldsymbol{r})\equiv {k}^{3/2}{\text{i}}^{n-1}\left(\begin{matrix}\hfill {\boldsymbol{N}}_{e,n,m}^{(+)}(k\boldsymbol{r})\hfill \\ \hfill -\text{i}{\boldsymbol{M}}_{e,n,m}^{(+)}(k\boldsymbol{r})\hfill \end{matrix}\right),\end{aligned}\end{equation} \tag{ 9 }$$

where k ≡ ω/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_in containing only the scattering medium. Inside this region, the RS fields can be expressed,

$$\begin{equation}{{\Psi}}_{\alpha }(\boldsymbol{r})=\sum\limits _{q=0}^{1}\sum\limits _{n=1}^{\infty }\sum\limits _{m=0}^{n}{d}_{q,n,m}^{(\alpha )}{{\Phi}}_{q,n,m}^{(1)}({\omega }_{\alpha },\boldsymbol{r})\qquad r{< }{R}_{\text{in}},\end{equation} \tag{ 10 }$$

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:

$$\begin{equation}\varepsilon (\omega )\equiv \frac{{\varepsilon }_{\text{s}}(\omega )}{{\varepsilon }_{\text{b}}},\qquad \mu (\omega )\equiv \frac{{\mu }_{\text{s}}(\omega )}{{\mu }_{\text{b}}},\qquad \rho (\omega )\equiv \sqrt{\frac{{\varepsilon }_{\text{s}}(\omega ){\mu }_{\text{s}}(\omega )}{{\varepsilon }_{\text{b}}{\mu }_{\text{b}}}},\end{equation} \tag{ 11 }$$

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

$$\begin{equation}\begin{aligned}{{\Phi}}_{h,n,m}^{(1)}\left(\omega ,\boldsymbol{r}\right)& \equiv {k}^{3/2}{\text{i}}^{n}\left(\begin{matrix}\hfill {\boldsymbol{M}}_{h,n,m}^{(1)}(k\boldsymbol{r}\rho (\omega ))\hfill \\ \hfill -\text{i}\sqrt{\frac{\varepsilon (\omega )}{\mu (\omega )}}{\boldsymbol{N}}_{n,m}^{(1)}(k\boldsymbol{r}\rho (\omega ))\hfill \end{matrix}\right)\\ {{\Phi}}_{e,n,m}^{(1)}(\omega ,\boldsymbol{r})& \equiv {k}^{3/2}{\text{i}}^{n-1}\left(\begin{matrix}\hfill {\boldsymbol{N}}_{n,m}^{(1)}(k\boldsymbol{r}\rho (\omega ))\hfill \\ \hfill -\text{i}\sqrt{\frac{\varepsilon (\omega )}{\mu (\omega )}}{\boldsymbol{M}}_{n,m}^{(1)}(k\boldsymbol{r}\rho (\omega ))\hfill \end{matrix}\right),\end{aligned}\end{equation} \tag{ 12 }$$

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, ω_α , would be to numerically solve the propagation of each ${{\Phi}}_{q,n,m}^{(1)}(\omega ,\boldsymbol{r})$ function through the inhomogeneous region, R_in < r < R_out, and then determine the values of ω_α 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_in and R_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, 30]).

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 reference [45]). As mentioned earlier, solving an eigenvalue problems like equation (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.

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 RSs (degenerate with respect to the m number). As discussed after equation (6) above, an additional 'quantum' index, ℓ, enumerates the different RS eigenfrequencies, ω_α within a given multipole set so that α(q, n, m, ℓ) designates one and only one RS. Due to the symmetries, the RS multipole developments of equations (8) and (10) simplify here to:

$$\begin{equation}{{\Psi}}_{\alpha }(\boldsymbol{r})=\begin{cases}_{\alpha }{{\Phi}}_{q,n,m}^{(1)}({\omega }_{\alpha },\boldsymbol{r});\quad r{< }R\quad \\ {c}_{\alpha }{{\Phi}}_{q,n,m}^{(+)}({\omega }_{\alpha },\boldsymbol{r});\quad r{ >}R\quad \end{cases},\end{equation} \tag{ 13 }$$

with the RS index, α, uniquely identified by the full set of quantum numbers, i.e. α(q, n, m, ℓ).

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 equations (9) and (12) for the Φ⁽⁺⁾ and Φ⁽¹⁾ basis functions, the continuity of the transverse electric and transverse magnetic fields [46] respectfully take the form of analytic linear relationships between internal field and external coefficients, which we henceforth write as, γ_α ≡ d_α /c_α :

$$\begin{equation}{\gamma }_{\alpha (e,n,\ell )}=\frac{{\rho }_{\alpha }}{{\varepsilon }_{\alpha }}\frac{{h}_{n}({z}_{\alpha })}{{j}_{n}({\rho }_{\alpha }{z}_{\alpha })}={\rho }_{\alpha }\frac{{\xi }_{n}^{\prime }({z}_{\alpha })}{{\psi }_{n}^{\prime }({\rho }_{\alpha }{z}_{\alpha })}\end{equation} \tag{ 14a }$$

$$\begin{equation}{\gamma }_{\alpha (h,n,\ell )}=\frac{{h}_{n}({z}_{\alpha })}{{j}_{n}({\rho }_{\alpha }{z}_{\alpha })}={\mu }_{\alpha }\frac{{\xi }_{n}^{\prime }({z}_{\alpha })}{{\psi }_{n}^{\prime }({\rho }_{\alpha }{z}_{\alpha })},\end{equation} \tag{ 14b }$$

where z_α ≡ k_α R, is the complex size parameter, j_n (z) the spherical Bessel functions, and h_n (z) the outgoing spherical Hankel functions, while ψ_n (z) ≡ zj_n (z) and ξ_n ≡ zh_n (z) are their respective Ricatti–Bessel function counterparts (cf appendix B). Since the respective pairs of continuity conditions for γ_α in equation (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 γ_α are satisfied.

The notations ɛ_α and μ_α , and ρ_α , introduced in equation (14) are shorthands for designating the constitutive parameters evaluated at the RS frequency, i.e.

$$\begin{equation}{\varepsilon }_{\alpha }\equiv \frac{{\varepsilon }_{\text{s}}({\omega }_{\alpha })}{{\varepsilon }_{\text{b}}},\qquad {\mu }_{\alpha }\equiv \frac{{\mu }_{\text{s}}({\omega }_{\alpha })}{{\mu }_{\text{b}}},\qquad {\rho }_{\alpha }\equiv \sqrt{\frac{{\varepsilon }_{\text{s}}({\omega }_{\alpha }){\mu }_{\text{s}}({\omega }_{\alpha })}{{\varepsilon }_{\text{b}}{\mu }_{\text{b}}}},\end{equation} \tag{ 15 }$$

which are normalized with respect to the exterior medium in accordance with the notation introduced in equation (2). To summarize, the exact expressions for the RSs of spherical particles can henceforth be written,

$$\begin{equation}{{\Psi}}_{\alpha }(\boldsymbol{r})=\frac{1}{{\mathcal{N}}_{\alpha }}\left(\begin{matrix}\hfill {\gamma }_{\alpha }{{\Phi}}_{q,n,m}^{(1)}({\omega }_{\alpha },\boldsymbol{r})\quad r{\leqslant}R\hfill \\ \hfill {{\Phi}}_{q,n,m}^{(+)}({\omega }_{\alpha },\boldsymbol{r})\quad r{\geqslant}R\hfill \end{matrix}\right),\end{equation} \tag{ 16 }$$

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 equation (16) and comparisons with expressions in the literature are given in appendix A.

One can readily verify that satisfying the equalities in equation (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 equation (85) of appendix D). We will show in section 3 that the determination of the values of z_α 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 ω_α in equation (37).

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 [40, 47]. We thus define 'bra' states as, $\left\langle {{\Psi}}_{\alpha }\right\vert \equiv \left\langle {\boldsymbol{E}}_{\alpha },{\boldsymbol{H}}_{\alpha }\right\vert$, without complex conjugation of the fields. The inner product of any two electromagnetic states, Ψ_α and Ψ_β is thus defined as the integral over a volume, $\mathcal{V}$, inside a surface ${\mathcal{S}}_{\mathrm{o}}$ that is sent to infinity:

$$\begin{align}\left\langle {{\Psi}}_{\beta }\vert {{\Psi}}_{\alpha }\right\rangle & \equiv \underset{\mathcal{S}\to \infty }{\mathrm{lim}}{\int }_{\mathcal{V}}\mathrm{d}\boldsymbol{r}{\left[{{\Psi}}_{\beta }(\boldsymbol{r})\right]}^{t}.{{\Psi}}_{\alpha }(\boldsymbol{r})\equiv {\int }_{{V}_{\infty }}\mathrm{d}\boldsymbol{r}\left\{{\boldsymbol{E}}_{\beta }(\boldsymbol{r})\cdot {\boldsymbol{E}}_{\alpha }(\boldsymbol{r})+{\boldsymbol{H}}_{\beta }(\boldsymbol{r})\cdot {\boldsymbol{H}}_{\alpha }(\boldsymbol{r})\right\}.\end{align} \tag{ 17 }$$

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 equation (17) is 'ill-defined', but the Gaussian regularization method assigns unambiguous, finite, values to RS products like those of equation (17) in a manner which allows the differential operator, $\mathbb{L}$ of equation (2), to remain symmetric when acting on regularized RS products, i.e.

$$\begin{equation}\left\langle {{\Psi}}_{\beta }\vert \mathbb{L}{{\Psi}}_{\alpha }\right\rangle =\left\langle \mathbb{L}{{\Psi}}_{\beta }\vert {{\Psi}}_{\alpha }\right\rangle .\end{equation} \tag{ 18 }$$

Applying the RS equations of motion in equation (2) to this relation provides an orthogonality relation for the product states, valid even in the case of temporally dispersive media,

$$\begin{equation}\frac{1}{{c}_{\text{b}}}\left[\left\langle {{\Psi}}_{\beta }\vert \mathbb{L}{{\Psi}}_{\alpha }\right\rangle -\left\langle \mathbb{L}{{\Psi}}_{\beta }\vert {{\Psi}}_{\alpha }\right\rangle \right]=0=\left\langle {{\Psi}}_{\beta }\right\vert \left[{k}_{\alpha }{\Gamma}({\omega }_{\alpha })-{k}_{\beta }{\Gamma}({\omega }_{\beta })\right]\left\vert {{\Psi}}_{\alpha }\right\rangle ,\end{equation} \tag{ 19 }$$

where Γ(ω) is the frequency dependent medium metric.

Although the inner product orthogonality relation of equation (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 equation (19) to be rewritten:

$$\begin{align}\frac{1}{\text{i}{c}_{\text{b}}}\left[{{\Psi}}_{\beta }^{t}.\mathbb{L}{{\Psi}}_{\alpha }-{\left(\mathbb{L}{{\Psi}}_{\beta }\right)}^{t}.{{\Psi}}_{\alpha }\right]& =\left[{\boldsymbol{E}}_{\beta },{\boldsymbol{H}}_{\beta }\right]\left[\begin{matrix}\hfill 0\hfill & \hfill \mathbf{\nabla }{\times}\hfill \\ \hfill \mathbf{\nabla }{\times}\hfill & \hfill 0\hfill \end{matrix}\right]\left[\begin{matrix}\hfill {\boldsymbol{E}}_{\alpha }\hfill \\ \hfill {\boldsymbol{H}}_{\alpha }\hfill \end{matrix}\right]\\ & \quad -\left[{\boldsymbol{E}}_{\alpha },{\boldsymbol{H}}_{\alpha }\right]\left[\begin{matrix}\hfill 0\hfill & \hfill \mathbf{\nabla }{\times}\hfill \\ \hfill \mathbf{\nabla }{\times}\hfill & \hfill 0\hfill \end{matrix}\right]\left[\begin{matrix}\hfill {\boldsymbol{E}}_{\beta }\hfill \\ \hfill {\boldsymbol{H}}_{\beta }\hfill \end{matrix}\right]\\ & ={\boldsymbol{H}}_{\beta }\cdot \mathbf{\nabla }{\times}{\boldsymbol{E}}_{\alpha }-{\boldsymbol{E}}_{\alpha }\cdot \mathbf{\nabla }{\times}{\boldsymbol{H}}_{\beta }\\ & \quad -\left[{\boldsymbol{H}}_{\alpha }\cdot \mathbf{\nabla }{\times}{\boldsymbol{E}}_{\beta }-{\boldsymbol{E}}_{\beta }\cdot \mathbf{\nabla }{\times}{\boldsymbol{H}}_{\alpha }\right]\\ & =\mathbf{\nabla }\cdot \left({\boldsymbol{E}}_{\alpha }{\times}{\boldsymbol{H}}_{\beta }-{\boldsymbol{E}}_{\beta }{\times}{\boldsymbol{H}}_{\alpha }\right).\end{align} \tag{ 20 }$$

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 figure 2.

Zoom In Zoom Out Reset image size

Integrating over the annular volume integral and applying the divergence theorem to equation (20) yields:

$$\begin{align}\frac{1}{{c}_{\text{b}}}{\int }_{{\mathcal{V}}_{a}}\mathrm{d}\boldsymbol{r}\left[{{\Psi}}_{\beta }^{t}.\mathbb{L}{{\Psi}}_{\alpha }-{\left(\mathbb{L}{{\Psi}}_{\beta }\right)}^{t}.{{\Psi}}_{\alpha }\right]& =\text{i}{\unicode{9675}}\hspace{-11.99995pt}\int \hspace{-8.00003pt}{\int }_{{\mathcal{S}}_{2}}\left[{\boldsymbol{E}}_{\alpha }({k}_{\alpha }\boldsymbol{r}){\times}{\boldsymbol{H}}_{\beta }({k}_{\beta }\boldsymbol{r})-{\boldsymbol{E}}_{\beta }({k}_{\beta }\boldsymbol{r})\right.\\ & \left.\quad {\times}{\boldsymbol{H}}_{\alpha }({k}_{\alpha }\boldsymbol{r})\right]\cdot {\mathrm{d}\boldsymbol{\mathcal{S}}}_{2}-\text{i}{\unicode{9675}}\hspace{-11.99995pt}\int \hspace{-8.00003pt}{\int }_{\mathcal{S}}\left[{\boldsymbol{E}}_{\alpha }({k}_{\alpha }\boldsymbol{r})\right.\\ & \left.\quad {\times}{\boldsymbol{H}}_{\beta }({k}_{\beta }\boldsymbol{r})-{\boldsymbol{E}}_{\beta }({k}_{\beta }\boldsymbol{r}){\times}{\boldsymbol{H}}_{\alpha }({k}_{\alpha }\boldsymbol{r})\right]\cdot \mathrm{d}\boldsymbol{\mathcal{S}}.\end{align} \tag{ 21 }$$

Letting the outer surface tend to infinity, the result can only be consistent with equation (19) provided that,

$$\begin{equation}\underset{{\mathcal{S}}_{2}\to \infty }{\mathrm{lim}}\enspace {\unicode{9675}}\hspace{-11.99995pt}\int \hspace{-8.00003pt}{\int }_{{\mathcal{S}}_{2}}\left[{\boldsymbol{E}}_{\alpha }({k}_{\alpha }\boldsymbol{r}){\times}{\boldsymbol{H}}_{\beta }({k}_{\beta }\boldsymbol{r})-{\boldsymbol{E}}_{\beta }({k}_{\beta }\boldsymbol{r}){\times}{\boldsymbol{H}}_{\alpha }({k}_{\alpha }\boldsymbol{r})\right]\cdot {\mathrm{d}\boldsymbol{\mathcal{S}}}_{2}=0.\end{equation} \tag{ 22 }$$

The vanishing of the surface integral, ${\mathcal{S}}_{2}$, at infinity arises from distribution theory and is discussed at length in reference [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 ${\text{e}}^{-\eta {r}^{2}}$ factor.

Equation (22) allows the orthogonality relation of equation (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}$:

$$\begin{align}{\int }_{\mathcal{V}}& \mathrm{d}\boldsymbol{r}\enspace {{\Psi}}_{\beta }^{t}(\boldsymbol{r}).\left[{k}_{\alpha }{\Gamma}({\omega }_{\alpha })-{k}_{\beta }{\Gamma}({\omega }_{\beta })\right].{{\Psi}}_{\alpha }(\boldsymbol{r})\\ & +\enspace \text{i}{\unicode{9675}}\hspace{-11.99995pt}\int \hspace{-8.00003pt}{\int }_{\mathcal{S}}\left[{\boldsymbol{E}}_{\alpha }({k}_{\alpha }\boldsymbol{r}){\times}{\boldsymbol{H}}_{\beta }({k}_{\beta }\boldsymbol{r})-{\boldsymbol{E}}_{\beta }({k}_{\beta }\boldsymbol{r}){\times}{\boldsymbol{H}}_{\alpha }({k}_{\alpha }\boldsymbol{r})\right]\cdot \mathrm{d}\boldsymbol{\mathcal{S}}=0.\end{align} \tag{ 23 }$$

For lossless particles, the medium metric, Γ, is frequency independent and real valued, so that equation (19) reads $({k}_{\alpha }-{k}_{\beta })\left\langle {{\Psi}}_{\beta }\right\vert {\Gamma}\left\vert {{\Psi}}_{\alpha }\right\rangle =0$, which for k_α ≠ k_β leads to a frequency independent eigenstate orthogonality condition for lossless systems:

$$\begin{equation}\begin{aligned}0& =\left\langle {{\Psi}}_{\beta }\right\vert {\Gamma}\left\vert {{\Psi}}_{\alpha }\right\rangle =\left\langle {{\Psi}}_{\beta }\vert {\Gamma}{{\Psi}}_{\alpha }\right\rangle =\left\langle {\Gamma}{{\Psi}}_{\beta }\vert {{\Psi}}_{\alpha }\right\rangle \\ & ={\int }_{{\mathcal{V}}_{\infty }}\mathrm{d}\boldsymbol{r}\left[{\boldsymbol{E}}_{\beta }(\boldsymbol{r})\enspace ,{\boldsymbol{H}}_{\beta }(\boldsymbol{r})\right].\left[\begin{matrix}\hfill \varepsilon (\boldsymbol{r})\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill -\mu (\boldsymbol{r})\hfill \end{matrix}\right].\left[\begin{matrix}\hfill {\boldsymbol{E}}_{\alpha }(\boldsymbol{r})\hfill \\ \hfill {\boldsymbol{H}}_{\alpha }(\boldsymbol{r})\hfill \end{matrix}\right]\\ & ={\int }_{{\mathcal{V}}_{\infty }}\mathrm{d}\boldsymbol{r}\left\{\varepsilon (\boldsymbol{r}){\boldsymbol{E}}_{\beta }(\boldsymbol{r})\cdot {\boldsymbol{E}}_{\alpha }(\boldsymbol{r})-\mu (\boldsymbol{r}){\boldsymbol{H}}_{\beta }(\boldsymbol{r})\cdot {\boldsymbol{H}}_{\alpha }(\boldsymbol{r})\right\},\end{aligned}\end{equation} \tag{ 24 }$$

which can of course also be expressed in the finite volume-surface form of equation (23).

It is worth remarking at this stage that the presence of the medium metric, Γ, 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 equation (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 RSs, Ψ_α and Ψ_β , developed according to equation (8), the volume integral of RS product in the region r > R_out (cf figure 1) is:

$$\begin{align}{\int }_{{R}_{\text{out}}}^{\infty }\mathrm{d}\boldsymbol{r}{\left[{{\Psi}}_{\beta }^{t}(\boldsymbol{r})\right]}^{t}.{\Gamma}.{{\Psi}}_{\alpha }(\boldsymbol{r})& =\sum\limits _{q,n,m}{c}_{\beta (q,n,m)}{c}_{\alpha (q,n,m)}{\int }_{{R}_{\text{out}}}^{\infty }\mathrm{d}r\enspace {{\Phi}}_{q,n,m}^{(+)}({k}_{\beta },\boldsymbol{r}).{\Gamma}.{{\Phi}}_{q,n,m}^{(+)}({k}_{\alpha },\boldsymbol{r}),\end{align} \tag{ 25 }$$

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 equation (25) are given in equations (81a) and (81d) of appendix C for α ≠ β and α = β respectively.

Inside a homogeneous region, r < R_in, lying completely inside a scattering object (cf figure 1), one can again develop the RS field on a multipole basis according to equation (10), and the RS volume integrals are:

$$\begin{align}{\int }_{0}^{{R}_{\text{in}}}\mathrm{d}\boldsymbol{r}\enspace {{\Psi}}_{\beta }^{t}(\boldsymbol{r}).{\Gamma}.{{\Psi}}_{\alpha }(\boldsymbol{r})& =\sum\limits _{q,n,m}{d}_{\beta (q,n,m)}{d}_{\alpha (q,n,m)}{\int }_{0}^{{R}_{\text{in}}}\mathrm{d}\boldsymbol{r}{\left[{{\Phi}}_{q,n,m}^{(1)}({\rho }_{\beta }{k}_{\beta },\boldsymbol{r})\right]}^{t}.{\Gamma}.{{\Phi}}_{q,n,m}^{(1)}({\rho }_{\alpha }{k}_{\alpha },\boldsymbol{r}),\end{align} \tag{ 26 }$$

with ordinary analytic expressions for the integrals on the right-hand side of equation (26) being given in equations (82) and (83) of appendix C. Numerical integration is still required for volume product integrals of RS fields in heterogeneous regions, like the region R_in < r < R_out in figure 1, but for spherical scatterers, one can take R_in = R_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, ɛ = 16 and null permeability contrast, μ = 1. A few low frequency electric mode RSs are graphically represented in figure 6 and given to high numerical precision in table 1. Since the sphere is lossless, the frequency independent RS orthogonality relation of equation (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 without 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.

Table 1. RS eigenvalues, z_α , and associated residues, ${\mathcal{R}}_{\alpha }=i/{\mathcal{N}}_{\alpha }^{2}$, for a non-dispersive dielectric medium with dielectric contrast of ɛ = 16 and no magnetic contrast, μ = 1.

α(q, n, [ℓ])	zα	${\mathcal{R}}_{\alpha }\equiv i/{\mathcal{N}}_{\alpha }^{2}$
(e, 1, [1])	1.0395 − i0.500 935	−0.236 682 + i0.231 492
(e, 1, [2])	1.052 73 − i0.072 3549	0.065 9905 − i0.057 9972
(e, 1, [3])	1.920 43 − i0.082 005	0.074 8408 − i0.028 2738
(e, 1, [4])	2.7227 − i0.073 007	0.002 794 37 − i0.068 3107
(e, 2, [0])	−i1.679 7303	i0.146 892
(e, 2, [1])	1.377 484 − i0.011 8433	0.001 846 13 − i0.011 8059
(e, 2, [2])	2.071 446 − i0.667 649	−0.305 381 + i0.277 002
(h, 1, [0])	−i1.250 038	i0.136 765
(h, 1, [1])	0.753 7823 − i0.024 0302	−0.006 017 59 − i0.022 9898
(h, 1, [2])	1.541 4631 − i0.045 9254	−0.039 4075 − i0.019 5948
(h, 2, [1])	0.870 513 − i1.752 59	−0.0521 306 + i0.140 046
(h, 2, [2])	1.095 7165 − i0.006 840 25	−0.000 482 964 − i0.006 816 78

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, $~{r}\equiv r/R$, a volume integral of the (complex valued) magnetic fields inside the sphere can then be evaluated as,

$$\begin{align}{\int }_{0}^{R}\mathrm{d}\boldsymbol{r}\enspace \mu {\boldsymbol{H}}_{e,n,m}({k}_{\beta }\boldsymbol{r})\cdot {\boldsymbol{H}}_{e,n,m}({k}_{\alpha }\boldsymbol{r})& =-\frac{{z}_{\alpha }^{3/2}{z}_{\beta }^{3/2}{\gamma }_{\alpha (e,n)}{\gamma }_{\beta (e,n)}}{{\mathcal{N}}_{\alpha }{\mathcal{N}}_{\beta }}\varepsilon \left\{{\int }_{0}^{1}{~{r}}^{2}{j}_{n}(\rho {z}_{\alpha }~{r}){j}_{n}(\rho {z}_{\beta }~{r})\mathrm{d}~{r}\right\}\end{align} \tag{ 27a }$$

$$\begin{equation}=\frac{{z}_{\alpha }^{3/2}{z}_{\beta }^{3/2}{\gamma }_{\alpha (e,n)}{\gamma }_{\beta (e,n)}}{{\mathcal{N}}_{\alpha }{\mathcal{N}}_{\beta }}\varepsilon \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\}\end{equation} \tag{ 27b }$$

$$\begin{equation}=\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}},\end{equation} \tag{ 27c }$$

where equation (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 equation (27a) as given by equation (85a) of [16] (cf also Watson reference [48], and equations (78) and (82c)). The final result of equation (27c) exploits the expressions of equation (14a) for γ_α(e,n), and employs algebraic manipulations that are only true when z_α and z_β are distinct RS size parameter eigenvalues (i.e. z_α ≠ z_β ).

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

$$\begin{align}{\int }_{R}^{\infty }\mathrm{d}\boldsymbol{r}\enspace {\boldsymbol{H}}_{e,n,m}({k}_{\beta }\boldsymbol{r})\cdot {\boldsymbol{H}}_{e,n,m}({k}_{\alpha }\boldsymbol{r})\to & -\frac{{z}_{\alpha }^{3/2}{z}_{\beta }^{3/2}}{{\mathcal{N}}_{\alpha }{\mathcal{N}}_{\beta }}\enspace \underset{\eta \to 0}{\mathrm{lim}}{\int }_{1}^{\infty }{~{r}}^{2}{h}_{n}({z}_{\alpha }~{r}){h}_{n}({z}_{\beta }~{r}){\text{e}}^{-\eta {~{r}}^{2}}\enspace \mathrm{d}r\end{align} \tag{ 28a }$$

$$\begin{equation}=-\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}({z}_{\alpha }){h}_{n}^{\prime }({z}_{\beta })}{{z}_{\alpha }^{2}-{z}_{\beta }^{2}},\end{equation} \tag{ 28b }$$

where we used equation (103a) in reference [16] to evaluate the integral in the first line of this equation (reproduced in equation (81e)). The 'killing' regularization, described in detail in reference [16] and appendix C, tames the integrand's exponentially divergent amplitude with a ${\text{e}}^{-\eta {r}^{2}}$ factor in order to yield the finite analytic expression of equation (28b) in the η → 0 limit. We ignored irrelevant overall sign factors in the above demonstration because it is now clear that the integrals in equations (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,

$$\begin{equation}{\int }_{{\mathcal{V}}_{\infty }}\mathrm{d}\boldsymbol{r}\enspace \mu (\boldsymbol{r}){\boldsymbol{H}}_{q,n,m}({k}_{\beta }\boldsymbol{r})\cdot {\boldsymbol{H}}_{q,n,m}({k}_{\alpha }\boldsymbol{r})=0,\end{equation} \tag{ 29 }$$

where we stopped refering only to the field type index 'e' (q = 0) in equation (29) since as shown in appendix C, analogous orthogonality relations hold for magnetic, ('h' (q = 1)), type RS modes.

In order to appreciate the complexity of what is accomplished in equation (29), the real and imaginary parts of a solid angle averaged RS magnetic product integrand of the first two electric dipole RSs at z_{α(e,n=1,ℓ=1)} ≃ 1.0395 − i0.5009 and z_{β(e,n=1,ℓ=2)} ≃ 1.0527 − i0.072 355 of a high index dielectric sphere (ɛ = 16) are plotted in figure 3 (cf table 1 and figure 6). Specifically, the real and imaginary parts of the angle averaged magnetic field product integrands given in equation (27a) (r < R) and equation (28a) (r > R) are plotted first for r/R in figure 3(a), up to 1.2 and then extended in figure 3(b) up to values of r/R = 25, where a Gaussian killing factor of η = 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 η in the region r/R > 1 is adopted.

Zoom In Zoom Out Reset image size

Analogous result hold 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,

$$\begin{equation}{\int }_{{\mathcal{V}}_{\infty }}\mathrm{d}\boldsymbol{r}\enspace \varepsilon (\boldsymbol{r}){\boldsymbol{E}}_{e,n,m}({k}_{\beta }\boldsymbol{r})\cdot {\boldsymbol{E}}_{e,n,m}({k}_{\alpha }\boldsymbol{r})=0,\end{equation} \tag{ 30 }$$

when α ≠ β, a result which is illustrated graphically by plotting the electric field integrands of equations (89a) and (89d) for r/R up to 1.2 and figure 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 equation (28b) since the required integrals were not given in standard references like that of Watson [48], 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 η, so as to be close to the η → 0 limit, but small values of η lead to large amplitude far-field oscillations (cf figure 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 [49].

Zoom In Zoom Out Reset image size

We can readily understand the numerical advantage of the complex contour method from figure 4, which plots the real part of equation (28), i.e. $\mathrm{R}\mathrm{e}\left[{r}^{2}{h}_{n}({z}_{\alpha }r){h}_{n}({z}_{\beta }r){\text{e}}^{-\eta {r}^{2}}\right]$ (with η = 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 π/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 a Gaussian regularization is superfluous for the rotated path which could be evaluated directly without regularization (i.e. for η = 0). Although it is visually clear from figures 3 and 4 that RS product integrands oscillate around zero outside the scatterer, they do not quite average to zero since the RS orthogonality relation of equation (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 figures 3(a) and (b)).

The utility of analytic Gaussian regularization formulas is further highlighted by remarking that if one also evaluates the electric field product integrals (cf equation (89) in appendix E), then after some algebra one finds, when k_α ≠ k_β , that:

$$\begin{align}& \left({k}_{\alpha }-{k}_{\beta }\right){\int }_{{\mathcal{V}}_{R}}\mathrm{d}\boldsymbol{r}\left[\varepsilon (\boldsymbol{r}){\boldsymbol{E}}_{q,n,m}({k}_{\beta }\boldsymbol{r})\cdot {\boldsymbol{E}}_{q,n,m}({k}_{\alpha }\boldsymbol{r})-\mu (\boldsymbol{r}){\boldsymbol{H}}_{q,n,m}({k}_{\beta }\boldsymbol{r})\cdot {\boldsymbol{H}}_{q,n,m}({k}_{\alpha }\boldsymbol{r})\right]\\ & \quad +{r}^{2}{(-1)}^{(n-q)}\frac{{k}_{\alpha }^{3/2}{k}_{\beta }^{3/2}}{{\mathcal{N}}_{\alpha }{\mathcal{N}}_{\beta }}\left[{h}_{n}(R{k}_{\alpha })\frac{{\xi }_{n}^{\prime }(R{k}_{\beta })}{R{k}_{\beta }}-{h}_{n}(R{k}_{\beta })\frac{{\xi }_{n}^{\prime }(R{k}_{\alpha })}{R{k}_{\alpha }}\right]=0,\end{align} \tag{ 31 }$$

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, equation (31) remains valid for any sphere surrounding the scatterer.

An examination of equation (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 equation (23). Although the explicit agreement between the analytic Gaussian regularization method with equation (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.

A variety of methodologies have been adopted to determine RS normalization [24–27, 50], which has given rise to 'different' normalization conditions being proposed for photonic problems, but Kristensen 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 we saw when we adopted equation (7) for the 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 equation (7) be proportional to their associated RS frequencies, ω_α , which since adimensional scalar products were desired, led us to adopt:

$$\begin{align}\left\langle {{\Psi}}_{\alpha }\right\vert {\left[\frac{\mathrm{d}}{\mathrm{d}\omega }(\omega {\Gamma})\right]}_{{\omega }_{\alpha }}\left\vert {{\Psi}}_{\alpha }\right\rangle & ={\int }_{{\mathcal{V}}_{\infty }}\mathrm{d}\mathbf{r}\left\{{\left.\frac{\mathrm{d}\left[\omega \varepsilon (\boldsymbol{r},\omega )\right]}{\mathrm{d}\omega }\right\vert }_{{\omega }_{\alpha }}{\boldsymbol{E}}_{\alpha }^{2}(\boldsymbol{r})-{\left.\frac{\mathrm{d}\left[\omega \mu (\boldsymbol{r},\omega )\right]}{\mathrm{d}\omega }\right\vert }_{{\omega }_{\alpha }}{\boldsymbol{H}}_{\alpha }^{2}(\boldsymbol{r})\right\}={z}_{\alpha },\end{align} \tag{ 32 }$$

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

$$\begin{equation}{\mathcal{E}}_{\alpha }^{\text{L}}\propto \left\langle {{\Psi}}_{\alpha }\right\vert {\Gamma}\left\vert {{\Psi}}_{\alpha }\right\rangle ={\int }_{{\mathcal{V}}_{\infty }}\mathrm{d}\mathbf{r}\left\{\varepsilon (\boldsymbol{r}){\boldsymbol{E}}_{\alpha }^{2}(\boldsymbol{r})-\mu (\boldsymbol{r}){\boldsymbol{H}}_{\alpha }^{2}(\boldsymbol{r})\right\}\;\;\text{no}\;\text{disp.}{=}\;\;{z}_{\alpha }.\end{equation} \tag{ 33 }$$

Despite the divergence of RS field amplitudes, the infinite volume integrals in equation (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 figure 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 equation (23) by k_α − k_β , and then take the limit k_α → k_β to obtain a finite volume expression of $\left\langle {{\Psi}}_{\alpha }\right\vert {\left.\frac{\mathrm{d}\left(k{\Gamma}\right)}{\mathrm{d}k}\right\vert }_{{k}_{\alpha }}\left\vert {{\Psi}}_{\alpha }\right\rangle$ being equal to a constant (which we assign as z_α ), so that a finite volume normalization expression reads:

$$\begin{align}\left\langle {{\Psi}}_{\alpha }\right\vert & {\left.\frac{\mathrm{d}\left(k{\Gamma}\right)}{\mathrm{d}k}\right\vert }_{{k}_{\alpha }}\left\vert {{\Psi}}_{\alpha }\right\rangle \to {\int }_{\mathcal{V}}\mathrm{d}\boldsymbol{r}\enspace {{\Psi}}_{\alpha }^{t}(\boldsymbol{r}).\frac{\mathrm{d}}{\mathrm{d}k}{\left[k{\Gamma}(\boldsymbol{r},k)\right]}_{{k}_{\alpha }}.{{\Psi}}_{\alpha }(\boldsymbol{r})+\text{i}{\unicode{9675}}\hspace{-11.99995pt}\int \hspace{-8.00003pt}{\int }_{\mathcal{S}}\left\{{\left[\frac{\mathrm{d}}{\mathrm{d}k}{\boldsymbol{E}}_{\alpha }(k\boldsymbol{r})\right]}_{{k}_{\alpha }}{\times}{\boldsymbol{H}}_{\alpha }\right.\\ & \left.\qquad -{\boldsymbol{E}}_{\alpha }{\times}{\left[\frac{\mathrm{d}}{\mathrm{d}k}{\boldsymbol{H}}_{\alpha }(k\boldsymbol{r})\right]}_{{k}_{\alpha }}\right\}\cdot \mathrm{d}\boldsymbol{\mathcal{S}}={z}_{\alpha },\end{align} \tag{ 34 }$$

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

The frequency derivative (or k derivative) of fields in the surface integral of equation (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 \mathbf{\nabla }}{k}$, as done in reference [41, 50]. 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:

$$\begin{equation}{\int }_{\mathcal{V}}\mathrm{d}\boldsymbol{r}\enspace {\boldsymbol{E}}_{\alpha }\cdot {\varepsilon }_{\alpha }(\boldsymbol{r})\cdot {\boldsymbol{E}}_{\alpha }+{\int }_{\mathcal{V}}\mathrm{d}\boldsymbol{r}\enspace {\boldsymbol{H}}_{\alpha }\cdot {\mu }_{\alpha }(\boldsymbol{r})\cdot {\boldsymbol{H}}_{\alpha }=\frac{{c}_{\text{b}}}{\text{i}{\omega }_{\alpha }}{\unicode{9675}}\hspace{-11.99995pt}\int \hspace{-8.00003pt}{\int }_{\mathcal{S}}\left({\boldsymbol{E}}_{\alpha }{\times}{\boldsymbol{H}}_{\alpha }\right)\cdot \mathrm{d}\boldsymbol{\mathcal{S}},\end{equation} \tag{ 35 }$$

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

$$\begin{align}& 2{\int }_{\mathcal{V}}\mathrm{d}\boldsymbol{r}\enspace {\boldsymbol{E}}_{\alpha }\cdot \frac{\mathrm{d}}{\mathrm{d}{k}^{2}}{\left[{k}^{2}\varepsilon ({k}^{2})\right]}_{{k}_{\alpha }^{2}}\cdot {\boldsymbol{E}}_{\alpha }+\frac{\text{i}}{{k}_{\alpha }}{\unicode{9675}}\hspace{-11.99995pt}\int \hspace{-8.00003pt}{\int }_{\mathcal{S}}\left({\boldsymbol{E}}_{\alpha }{\times}{\boldsymbol{H}}_{\alpha }\right)\cdot \mathrm{d}\boldsymbol{\mathcal{S}}\\ & \quad +\frac{\text{i}}{{k}_{\alpha }}{\unicode{9675}}\hspace{-11.99995pt}\int \hspace{-8.00003pt}{\int }_{\mathcal{S}}\left\{\left[\boldsymbol{r}\cdot \mathbf{\nabla }{\boldsymbol{E}}_{\alpha }\right]{\times}{\boldsymbol{H}}_{\alpha }-{\boldsymbol{E}}_{\alpha }{\times}\left[\boldsymbol{r}\cdot \mathbf{\nabla }{\boldsymbol{H}}_{\alpha }\right]\right\}\cdot \mathrm{d}\boldsymbol{\mathcal{S}}={z}_{\alpha },\end{align} \tag{ 36 }$$

where we used the replacement $\frac{\mathrm{d}}{\mathrm{d}\omega }{\left[\varepsilon (\omega )\right]}_{{\omega }_{\alpha }}\to 2\frac{\mathrm{d}}{\mathrm{d}{k}^{2}}{\left[\varepsilon ({k}^{2})\right]}_{{k}_{\alpha }^{2}}$. The surface terms in equation (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 [52].

Even though the normalization formulas of equations (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 equation (34) of an arbitrary spherical surface of radius R is analytically identical to the formula of equation (94) 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 equation (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 equation (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 reference [16], we find that the normalization factors introduced in equation (16) have the following analytic expressions:

$$\begin{align}{\mathcal{N}}_{\alpha \enspace \enspace (h,n,\ell )}^{2}& =({\varepsilon }_{\alpha }-1){\xi }_{n}^{2}({z}_{\alpha })+({\mu }_{\alpha }-1)\left\{{\left[{\xi }_{n}^{\prime }({z}_{\alpha })\right]}^{2}+\frac{n(n+1){h}_{n}^{2}({z}_{\alpha })}{{\mu }_{\alpha }}\right\}\\ & \qquad +\frac{{\omega }_{\alpha }}{2}\left\{{{\Xi}}_{\alpha (h,n,\ell )}^{(-)}{\left.\frac{\mathrm{d}}{\mathrm{d}\omega }\enspace \mathrm{ln}\enspace \varepsilon (\omega )\right\vert }_{{\omega }_{\alpha }}+{{\Xi}}_{\alpha (h,n,\ell )}^{(+)}{\left.\frac{\mathrm{d}}{\mathrm{d}\omega }\enspace \mathrm{ln}\enspace \mu (\omega )\right\vert }_{{\omega }_{\alpha }}\right\}\end{align} \tag{ 37a }$$

$$\begin{align}{\mathcal{N}}_{\alpha (e,n,\ell )}^{2}& =({\mu }_{\alpha }-1){\xi }_{n}^{2}({z}_{\alpha })+({\varepsilon }_{\alpha }-1)\left\{{[{\xi }_{n}^{\prime }({z}_{\alpha })]}^{2}+\frac{n(n+1)}{{\varepsilon }_{\alpha }}{h}_{n}^{2}({z}_{\alpha })\right\}\\ & \quad +\frac{{\omega }_{\alpha }}{2}\left\{{{\Xi}}_{\alpha (e,n,\ell )}^{(+)}{\left.\frac{\mathrm{d}}{\mathrm{d}\omega }\enspace \mathrm{ln}\enspace {\varepsilon }_{\alpha }(\omega )\right\vert }_{{\omega }_{\alpha }}+{{\Xi}}_{\alpha (e,n,\ell )}^{(-)}{\left.\frac{\mathrm{d}}{\mathrm{d}\omega }\enspace \mathrm{ln}\enspace {\mu }_{\alpha }(\omega )\right\vert }_{{\omega }_{\alpha }}\right\}\end{align} \tag{ 37b }$$

with the definitions:

$$\begin{equation}{{\Xi}}_{\alpha (h,n,\ell )}^{({\pm})}\equiv {\mu }_{\alpha }{\left[{\xi }_{n}^{\prime }({z}_{\alpha })\right]}^{2}+{\varepsilon }_{\alpha }{\xi }_{n}^{2}({z}_{\alpha })-\frac{n(n+1){h}_{n}^{2}({z}_{\alpha })}{{\mu }_{\alpha }}{\pm}{h}_{n}({z}_{\alpha }){\xi }_{n}^{\prime }({z}_{\alpha })\end{equation} \tag{ 37c }$$

$$\begin{equation}{{\Xi}}_{\alpha (e,n,\ell )}^{({\pm})}\equiv {\varepsilon }_{\alpha }{\left[{\xi }_{n}^{\prime }({z}_{\alpha })\right]}^{2}+{\mu }_{\alpha }{\xi }_{n}^{2}({z}_{\alpha })-\frac{n(n+1){h}_{n}^{2}({z}_{\alpha })}{{\varepsilon }_{\alpha }}{\pm}{h}_{n}({z}_{\alpha }){\xi }_{n}^{\prime }({z}_{\alpha }).\end{equation} \tag{ 37d }$$

For dispersionless media, light scattering becomes scale invariant and equation (37) simplifies to expressions which only depend on the size parameter of the mode, z_α , and the relative constitutive parameters, ɛ and μ (now real-valued) as was previously derived in reference [17]:

$$\begin{align}{\mathcal{N}}_{\alpha (e,n,\ell )}^{2}& \;\;\text{no}\;\text{disp.}{=}\;\;\enspace {\xi }_{n}^{2}({z}_{\alpha })(\mu -1)+\left\{{[{\xi }_{n}^{\prime }({z}_{\alpha })]}^{2}+\frac{n(n+1)}{\varepsilon }{h}_{n}^{2}({z}_{\alpha })\right\}(\varepsilon -1)\end{align} \tag{ 38a }$$

$$\begin{align}{\mathcal{N}}_{\alpha (h,n,\ell )}^{2}& \;\;\text{no}\;\text{disp.}{=}\;\;\enspace {\xi }_{n}^{2}({z}_{\alpha })(\varepsilon -1)+\left\{{\left[{\xi }_{n}^{\prime }({z}_{\alpha })\right]}^{2}+\frac{n(n+1)}{\mu }{h}_{n}^{2}({z}_{\alpha })\right\}(\mu -1).\end{align} \tag{ 38b }$$

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 equation (33) (after angular integration) for a lossless dielectric of ɛ = 16 as shown in figure 5. The convergence factor was simply set to η = 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₂, z₃, and z₄ RSs. For such low loss RSs, a numerical 'verification' of equations (33) or (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 equation (33) as one can guess by visual inspection of figures 5(b)–(d).

Zoom In Zoom Out Reset image size

The situation is quite different however for poorly confined RSs with large imaginary parts, like z₁ ≃ 1.04 − i0.5. As one can see in figure 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 equation (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 equations (32), (34) and (36) to z_α , 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 ${~{{\Psi}}}_{\alpha }(\boldsymbol{r})\equiv {{\Psi}}_{\alpha }(\boldsymbol{r})/{z}_{\alpha }^{1/2}$ for which the RS normalization reads:

$$\begin{equation}\left\langle {~{{\Psi}}}_{\alpha }\left\vert {\left[\frac{\mathrm{d}}{\mathrm{d}\omega }(\omega {\Gamma})\right]}_{{\omega }_{\alpha }}\right\vert {~{{\Psi}}}_{\alpha }\right\rangle =1,\end{equation} \tag{ 39 }$$

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_α [23, 24, 26, 39, 53]. 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 [54–56]. A typical definition for V_α derived along the LDOS Green function approach is:

$$\begin{equation}{V}_{\alpha }\mathrm{R}\mathrm{e}\mathrm{f}.[\mathrm{24}]{=}\frac{\langle \langle {\tilde {\mathbf{f}}}_{\alpha }\vert {\tilde {\mathbf{f}}}_{\alpha }\rangle \rangle }{\varepsilon ({\boldsymbol{r}}_{\text{c}}){\tilde {\mathbf{f}}}_{\alpha }^{2}({\boldsymbol{r}}_{\text{c}})},\end{equation} \tag{ 40 }$$

where r _c is a 'chosen' position (often the field maxima) while ${\tilde {\mathbf{f}}}_{\alpha }(\boldsymbol{r})$ is one of the commonly employed photonic notations [24] for the RS electric field, and the double '⟨⟨...⟩⟩' is an inner product notation reference [24]. Following the usual LDOS route to mode volume in the framework of this work leads to,

$$\begin{equation}{V}_{\alpha }\propto \frac{\left\langle {{\Psi}}_{\alpha }\right\vert {\left[\frac{\mathrm{d}}{\mathrm{d}\omega }(\omega {\Gamma})\right]}_{{\omega }_{\alpha }}\left\vert {{\Psi}}_{\alpha }\right\rangle }{\varepsilon ({\boldsymbol{r}}_{\text{c}}){\mathbf{E}}_{\alpha }^{2}({\boldsymbol{r}}_{\text{c}})},\end{equation} \tag{ 41 }$$

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

As discussed in the previous paragraph and the references therein, given that V_α is generally complex valued, there appears to be little utility in trying to interpret it as a true mode 'volume', and it is generally considered preferable to interpret it as the inverse normalized field squared. Mode volume assignments along the lines of equations (40) and (41) tell us that |V_α | is proportional ${\left\vert {\mathcal{N}}_{\alpha }\right\vert }^{2}$, so we generally expect small values of $\left\vert {\mathcal{N}}_{\alpha }\right\vert$ to correspond to 'small' mode volumes, or in other words, stronger interactions with an incident field. One indeed sees from table 2 that poorly confined, 'structural' modes like ${z}_{\ell =1}^{(e)}$ (cf figure 5(a)) have smaller $\left\vert {\mathcal{N}}_{\alpha }\right\vert$, than the more 'resonant' modes illustrated in the other graphics of figure 5.

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

R (nm)	ℓ	Disp.	λα (nm)	ɛα	${\mathcal{R}}_{\alpha }\equiv i/{\mathcal{N}}_{\alpha }^{2}$
100	1	(52b)	606.976 + i239.112	−13.7606 − i12.5419	−0.217 942 + i0.034 889
100	1	(53)	592.219 + i210.126	−5.904 24 − i12.9379	−0.309 192 − i0.005 709 58
80	1	(52b)	505.163 + i174.433	−9.480 69 − i7.588 28	−0.2111 − i0.003 857 08
80	1	(53)	547.731 + i79.8811	−3.642 00 − i3.522 01	0.017 306 − i0.100 09
100	1	(52a)	600.211 + i231.333	−11.1356 − i14.0631	−0.236 629 + i0.026 6621
⋅	2	⋅	290.678 + i33.6325	0.704 812 − i0.966 352	0.178 962 + i0.124 692
⋅	3	⋅	217.845 + i20.8224	2.580 05 − i0.449 824	0.013 9492 − i0.235 957
80	1	⋅	500.306 + i156.443	−6.780 66 − i7.895 19	−0.235 268 − i0.050 9908
⋅	2	⋅	281.965 + i50.9047	1.030 36 − i1.441 87	0.127 46 + i0.248 96
⋅	3	⋅	192.72 + i20.4675	3.110 19 − i0.394 101	0.103 72 − i0.210 644

Inspection of equation (38) appears to lead to the conclusion that weak constitutive contrasts lead to small values of $\left\vert {\mathcal{N}}_{\alpha }\right\vert$). However, the ${\mathcal{N}}_{\alpha }$ factors in equation (38) do not at all vanish in the limit of zero contrast since the vanishing (ɛ − 1) and (μ − 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.

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 [46] for the diagonal 'T-matrix' elements; T_h,n ≡ −b_n and T_e,n ≡ −a_n , and Ω_h,n = c_n and Ω_e,n = d_n for the internal field coefficients.

The meromorphic expansion of the T_n functions in terms of frequency, ω, have been determined to be [17, 18, 22, 37]:

$$\begin{equation}{T}_{q,n}(\omega ,R)=\frac{{A}_{q,n}{\text{e}}^{-2\text{i}z}-1}{2}+{\text{e}}^{-2\text{i}z}{A}_{q,n}\sum\limits _{\ell =-\infty }^{\infty }\frac{{r}_{\alpha }}{z-{z}_{\alpha }}\end{equation} \tag{ 43a }$$

$$\begin{equation}?{\simeq }\underset{\mathrm{n}\mathrm{o}\mathrm{n}-\mathrm{r}\mathrm{e}\mathrm{s}\mathrm{o}\mathrm{n}\mathrm{a}\mathrm{n}\mathrm{t}\enspace \mathrm{c}\mathrm{o}\mathrm{n}\mathrm{t}\mathrm{r}\mathrm{i}\mathrm{b}\mathrm{u}\mathrm{t}\mathrm{i}\mathrm{o}\mathrm{n}}{\underbrace{\frac{{A}_{q,n}{\text{e}}^{-2\text{i}z}-1}{2}}}+\enspace \underset{\text{Truncated}\;\text{RS}\;\text{expansion}}{\underbrace{{\text{e}}^{-2\text{i}z}{A}_{q,n}\sum\limits _{\ell =-{L}_{\mathrm{max}}}^{{L}_{\mathrm{max}}}\frac{{r}_{\alpha }}{z-{z}_{\alpha }}}},\end{equation} \tag{ 43b }$$

which takes the form of a non-resonant contribution to which one adds the spectral RS contributions. (z_α = k_α R with α(q, n, ℓ) being determined by multipole indices, q, n, and the discrete RS index, ℓ). The question mark on the truncated RS approximation in equation (43b) highlights the fact that the obligatory truncation of an infinite RS sum should be addressed with care (see the end of this section and section 4.4 for further discussion).

Recent works have proposed to derive the non-resonant contributions by summing static-mode contributions [42, 43], but we have found it more expedient to determine the non-resonant contributions analytically by invoking mathematical constraints derived in the context of the S/T-matrix response functions [17, 18]. An advantage of our formulation is that it holds true even for Drude-type constitutive parameters where static mode formulations do not appear to be applicable. The relationship between these approaches appears to merit further investigation [43].

The first, non-resonant, term on the right-hand side of equation (43) is a function of only the particle size parameter, z = kR = ωR/c. The A_q,n in equation (43), are multipole sign factors, which for q = 0 (magnetic) or q = 1 (electric), modes are:

$$\begin{equation}{A}_{q,n}\equiv {(-1)}^{n-q}.\end{equation} \tag{ 44 }$$

The residues, r_α , of the RSs in the RS expansion of the T-matrix in equation (43) are,

$$\begin{equation}{r}_{\alpha }=\frac{\text{i}{\text{e}}^{2\text{i}{z}_{\alpha }}}{{\mathcal{N}}_{\alpha }^{2}}\equiv {\text{e}}^{2\text{i}{z}_{\alpha }}{\mathcal{R}}_{\alpha },\end{equation} \tag{ 45 }$$

where the ${\text{e}}^{2\text{i}{z}_{\alpha }}$ phase factor does not affect the T-matrix residue equation (43) since it is canceled by the e^−2iz multiplicative factor so the RS residues are indeed ${\mathcal{R}}_{\alpha }=\mathit{\text{i}}/{\mathcal{N}}_{\alpha }^{2}$ as expected.

The internal field coefficients of Mie theory, Ω_q,n , are quite similar to the 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):

$$\begin{equation}{{\Omega}}_{q,n}(\omega ,R)={\text{e}}^{-2\text{i}z}{A}_{q,n}\sum\limits _{\ell =-\infty }^{\infty }\frac{{r}_{\alpha }^{({\Omega})}}{z-{z}_{\alpha }},\end{equation} \tag{ 46 }$$

As one could deduce from equation (7) and the RS expression in equation (16), the residues of internal RS coefficients, are related to the T-matrix residues in equation (46), as ${r}_{\alpha }^{({\Omega})}={\gamma }_{\alpha }{r}_{\alpha }$.

Since the non-resonant contributions were determined analytically here, we could expect to be able to retrieve the Mie coefficients from the RS spectral expansion in equation (43), but upon closer inspection we realize that there is no reason to presume that the neglected, high frequency, RSs will make negligible contributions to the Mie coefficients (which up to a sign difference are matrix elements of the T-matrix which in turn is a reformulation of the Green's function—as seen in equation (42)). In order to demonstrate this lack of 'convergence', we developed codes that determined the RSs in a given frequency window, and one of our initial 'solutions' to the convergence problem, for dielectrics, was to simply to include hundreds of RSs when reconstructing Mie coefficients using equation (43b).

We should keep in mind however that the initial interest of RSs was to characterize a physical system with a modest number of relevant RSs, and the necessity of including hundreds of RSs in spectral expansions would considerably diminish the utility of this formulation, particularly in the case of non-spherical objects where the asymptotics of high frequency RSs will be less straightforward. This convergence difficulty turned out to be even more problematic when treating dispersive scatterers as we will see in section 4.4.

Given the expressions for Ω_n , and T_n in section 3.1 and equation (46) in terms of the N_q,n and D_q,n functions defined in equations (86) and (87) of appendix D.2, the residues in equation (46) and section 3.1 can alternatively be calculated from Mie theory, as:

$$\begin{equation}{r}_{\alpha }^{({\Omega})}(R)=\frac{{\text{e}}^{2\text{i}{z}_{\alpha }}}{{\left.\frac{c}{R}\frac{\mathrm{d}}{\mathrm{d}\omega }{D}_{q,n}(\omega ,R)\right\vert }_{\omega ={\omega }_{\alpha }}},\qquad {r}_{\alpha }(R)=-\frac{{N}_{q,n}\left({\omega }_{\alpha },R\right){\text{e}}^{2\text{i}{z}_{\alpha }}}{{\left.\frac{c}{R}\frac{\mathrm{d}}{\mathrm{d}\omega }{D}_{q,n}(\omega ,R)\right\vert }_{\omega ={\omega }_{\alpha }}},\end{equation} \tag{ 47a }$$

which can be rewritten as,

$$\begin{equation}{r}_{\alpha }^{({\Omega})}(R)={\text{e}}^{2\text{i}{z}_{\alpha }}\underset{\omega \to {\omega }_{\alpha }}{\mathrm{lim}}\enspace \frac{z-{z}_{\alpha }}{{D}_{q,n}(\omega ,R)},\qquad {r}_{\alpha }(R)=-{N}_{q,n}({\omega }_{\alpha },R){r}_{\alpha }^{({\Omega})}.\end{equation} \tag{ 47b }$$

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

We remark from the second equality in equation (47b) that the relation between the T-matrix residues, r_n,α , 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 equation (14) such that the expressions of N_e,n and N_h,n can be rewritten:

$$\begin{equation}\begin{aligned}-{N}_{e,n}({\omega }_{\alpha },R)={z}_{\alpha }\frac{{\psi }_{n}^{\prime }({\rho }_{\alpha }{z}_{\alpha }){j}_{n}({z}_{\alpha })-{\varepsilon }_{\alpha }{j}_{n}({\rho }_{\alpha }{z}_{\alpha }){\psi }_{n}^{\prime }({z}_{\alpha })}{\text{i}{\rho }_{\alpha }}\to \frac{{\psi }_{n}^{\prime }({\rho }_{\alpha }{z}_{\alpha })}{{\rho }_{\alpha }{\xi }_{n}^{\prime }({z}_{\alpha })}=\frac{1}{{\gamma }_{e,n,\alpha }},\end{aligned}\end{equation} \tag{ 48a }$$

for the electric modes and,

$$\begin{equation}\begin{aligned}-{N}_{h,n}({\omega }_{\alpha },R)={z}_{\alpha }\frac{{\psi }_{n}^{\prime }({\rho }_{\alpha }{z}_{\alpha }){j}_{n}({z}_{\alpha })-{\mu }_{\alpha }{j}_{n}({\rho }_{\alpha }{z}_{\alpha }){\psi }_{n}^{\prime }({z}_{\alpha })}{\text{i}{\mu }_{\alpha }}\to \frac{{\psi }_{n}^{\prime }({\rho }_{\alpha }{z}_{\alpha })}{{\mu }_{\alpha }{\xi }_{n}^{\prime }({z}_{\alpha })}=\frac{1}{{\gamma }_{h,n,\alpha }},\end{aligned}\end{equation} \tag{ 48b }$$

for the magnetic modes. In both equations (48a) and (48b) we invoked the Wronskian relation,

$$\begin{equation}{\xi }_{n}^{\prime }(z){\psi }_{n}(z)-{\xi }_{n}(z){\psi }_{n}^{\prime }(z)=i,\end{equation} \tag{ 49 }$$

and recalled the definitions of the γ_α functions defined in equation (14).

Mie theory therefore analytically validates the relation,

$$\begin{equation}{r}_{\alpha }=\frac{{r}_{\alpha }^{({\Omega})}}{{\gamma }_{\alpha }},\end{equation} \tag{ 50 }$$

that we deduced after equation (46) 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}}_{\alpha }^{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.

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

Numerical results for a few low-order modes of a dielectric with no magnetic permeability contrast, μ = 1 and high, ɛ = 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_α 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 ℓ = 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_α and normalizations, ${\mathcal{R}}_{\alpha }=i/{\mathcal{N}}_{\alpha }^{2}$ 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.

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

• (a)  
  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 RS size parameters, z_α , by numerically solving the transcendental equation (14) (or equivalently equation (51)). The sufficient number of RSs depends on the particle radius R (or frequency of interest) and the dispersion relations of the constitutive parameters, but might 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 figures 6(b) and (e)).
• (b)  
  Calculate the square of the complex-valued RS norms, ${\mathcal{N}}_{\alpha (q,n,\ell )}^{2}$, using equation (37).
• (c)  
  Insert the values of ${\mathcal{N}}_{\alpha (q,n,\ell )}^{2}$ found in step (b) into equation (45) to determine the residues, r_α(q,n,ℓ), for the meromorphic expression of section 3.1 for the scattered field Mie coefficients, a_n = −T_e,n and b_n = −T_h,n . Cross section contributions, σ_n , can then be determined by the standard formulas recalled in equation (88) of appendix D.2.
• (d)  
  If one needs to determine the fields inside the sphere, use the values of r_α found in step (c) and the γ_α(q,n,ℓ) coefficients of equation (14) to calculate the residues, ${r}_{\alpha (q,n,\ell )}^{\left({\Omega}\right)}$, of equation (50) that appear in the meromorphic expression of equation (46), for the internal field Mie coefficients c_n = Ω_h,n and d_n = Ω_e,n (defined in equation (85b)).

The results displayed in table 1 of section 4.1 were obtained by applying steps (a) and (b) above. An illustration of carrying out step (c) above, with the same material parameters (ɛ = 16, μ = 1) used in figure 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 Im(z) on the vertical axis, while the associated contributions to the cross section efficiencies, Q ≡ σ/(πR²), are plotted on the positive vertical scale. In figures 6(c) and (f), we verified that this procedure reproduces Mie theory up to arbitrary computational accuracy for all multipole orders when large numbers of resonant states are included.

Zoom In Zoom Out Reset image size

The curves in figure 5 were generated both from RS calculations and Mie theory, but are identical since they agreed here up to six figure accuracy, but at the cost however of including hundreds of RSs in equation (43). 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.

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 the comprehensive collection due to Palik [64]. 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, ω, for gold [23] and silver [65] as:

$$\begin{align}\varepsilon (\omega )& ={\varepsilon }_{\infty }-\frac{{\omega }_{p}^{2}}{{\omega }^{2}+\text{i}\omega {{\Gamma}}_{D}}.\\ {\varepsilon }_{\infty }^{\text{Ag}}& =5,\qquad {\omega }_{p}^{\text{Ag}}=1.35{\times}1{0}^{16}\enspace ({\mathrm{s}}^{-1})\enspace (8.89\enspace \mathrm{e}\mathrm{V}),\qquad {{\Gamma}}_{D}^{\text{Ag}}=5.88{\times}1{0}^{13}\enspace ({\mathrm{s}}^{-1})\enspace (0.0387\enspace \mathrm{e}\mathrm{V})\enspace \end{align} \tag{ 52a }$$

$$\begin{equation}{\varepsilon }_{\infty }^{\text{Au}}=1,\qquad {\omega }_{p}^{\text{Au}}=1.26{\times}1{0}^{16}\enspace ({\mathrm{s}}^{-1})\enspace (8.29\enspace \mathrm{e}\mathrm{V}),\qquad {{\Gamma}}_{D}^{\text{Au}}=1.41{\times}1{0}^{14}\enspace ({\mathrm{s}}^{-1})\enspace (0.0928\enspace \mathrm{e}\mathrm{V}).\end{equation} \tag{ 52b }$$

The comparison with experimental data given in figure 13 of appendix G shows that the Drude model works relatively well for silver at wavelengths longer than around 0.30 μm, while the comparison of the Drude model for gold in figure 14 shows that the Drude model for gold only works well for wavelengths longer than 0.60 μ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 equation (52b). For example, Sikdar and Kornyshev [66] give the parameters for a model for gold taking into account inter-band transitions via two additional resonances:

$$\begin{equation}\varepsilon (\omega )={\varepsilon }_{\infty }-\frac{{\omega }_{pD}^{2}}{{\omega }^{2}+\text{i}\omega {{\Gamma}}_{D}}-{s}_{1}\frac{{\omega }_{p1,L}^{2}}{{\omega }^{2}-{\omega }_{p1,L}^{2}+\text{i}{{\Gamma}}_{1,L}\omega }-{s}_{2}\frac{{\omega }_{p2,L}^{2}}{{\omega }^{2}-{\omega }_{p2,L}^{2}+\text{i}{{\Gamma}}_{2,L}\omega }.\end{equation} \tag{ 53 }$$

Here the parameters are given in eV: _∞ = 5.9752, ℏω_pD = 8.8667 eV, ℏΓ_D = 0.037 99 eV, s₁ = 1.76, ℏω_p1,L = 3.6 eV, ℏΓ_1,L = 1.3 eV, s₂ = 0.952, ℏω_p2,L = 2.8 eV and ℏΓ_2,L = 0.737 eV. The conversion of ω to electron volts is achieved by dividing the frequency results by e/ℏ = 1.519 27 × 10¹⁵. As can be seen from figure 14, this model works well down to around 0.40 μm. An interesting feature of the dispersion model of equation (53) is that it supports bulk plasmons [67]. These occur for complex wavelengths or frequencies at which ɛ(ω) = 0. For the dispersion model of equation (53) they are at: λ_Lg = 0.257 778 + i0.020 6709 μm, λ_Lg = 0.395 618 + i0.053 6046 μm, λ_Lg = 0.502 518 + i0.048 786 μ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 any longitudinal 'RSs'.

Figure 7 shows the modulus of the dispersion equation of equation (51a) for n = 1 and for a gold sphere of radius 100 nm, both for the Drude model of equation (52b) and for the more elaborate Drude-Lorentz model of equation (53). The former has its minimum at λ_α = 0.606 976 + i0.239 112 μm, in good agreement with the value reported in Sauvan et al [23]. The latter has a slightly different value: λ_α = 0.592 227 + i0.210 097 μm. From table 2 the values of ɛ_α at the resonance and consequently the normalization factor disagree more significantly for the more elaborate model.

Zoom In Zoom Out Reset image size

Figure 8 shows the inverse modulus of the n = 1 dispersion relation of equation (51) for a gold sphere of radius 80 nm. The plot is more complicated for the more elaborate dispersion model of equation (53) 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 λ_α = 0.505 163 + i0.174 433 μm. With the model of equation (53), this becomes λ_α = 0.547 815 + i0.080 062 μm. Note the significant decrease in the imaginary part around this wavelength, as visible in figure 14, since the Drude model is only a crude approximation to the experimental data.

Zoom In Zoom Out Reset image size

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 80 nm. In order to calculate the mode normalization factors, the numerical differentiation, described in equation (47), is a valid alternative to equation (37) for finding the necessary dispersion derivatives. The data given in table 2 for the resonances of silver spheres of radius 100 and 80 nm shows complex wavelengths and normalization factors which are broadly similar to the corresponding values for gold.

We argued at the end of section 3.1 that one cannot necessarily expect the Mie coefficients of equation (43a), with their infinite sets of RS poles, to be accurately approximated by the directly truncated RS sums in equation (43b). In order to better understand this issue, let us first see what happens if we adopt the direct truncation of equation (43b) for a particularly simple physical example. Specifically, let us adopt the gold Drude model of equation (52b) and an R = 80 nm sphere, for which a well isolated 'principal' RS is illustrated graphically in figure 8 (and given numerically in table 2-line 3). This isolated, long wavelength RS has the promising property of lying near the broad, red-shifted plasmon 'resonance' peak in the electric dipole extinction cross section, ${\sigma }_{\mathrm{e}\mathrm{x}\mathrm{t},1}^{(\mathit{\text{e}})}$, calculated from Mie theory and plotted as the solid blue lines in figures 9(a) and (b).

Zoom In Zoom Out Reset image size

Despite the apparent simplicity of the physical situation, truncating the RS expansion in equation (43b) to a single independent RS state is far from being satisfactory (dotted green curve in figure 9(a): L_max = 1). If we try to 'improve' the approximation by adding additional RSs, as shown by a L_max = 4 calculation in figure 9(b) (dotted green curve) the agreement only worsens and quickly loses any semblance of physical relevance with predictions of large negative cross sections and unitarity violations. A solution to this problem will be described in the following paragraphs, and the results of this new, 'optimized', truncation technique are plotted as dashed red curves in figures 9(a) and (b) for L_max = 1 and L_max = 4 respectively. The word 'optimized' is used to indicate that this improved methodology tends towards the exact results with increasing values of L_max as desired.

The cause for the 'failure' of the direct spectral expansion was traced to the fact that the Mie coefficients resulting from the 'infinite' RS series, has an analytic structure that a direct truncation of the RS series in equation (43) (using analytic normalizations/residues) does not preserve. One can solve this problem by constraining the truncated RS spectral expansions to satisfy the same unitarity requirements as the exact T_q,n function deduced from Mie theory. This was done by appealing to the S-matrix, related to the T-matrix by $\mathbb{S}\equiv \mathbb{I}+2\mathbb{T}$, ($\mathbb{I}$ being the identity operator). The key observation is that the zeros of the S-matrix, denoted ${\omega }_{\alpha }^{(-)}$, satisfy the same equation as the RSs, i.e. equation (4), but with incoming boundary conditions replacing the outgoing RS boundary conditions (in equations where the RS eigenvalues also appear, the RS eigenvalues will be denoted ${\omega }_{\alpha }^{(+)}$). For lossless dielectrics, the ${\omega }_{\alpha }^{(-)}$ are simply the complex conjugates of the RS eigenstates, ${\omega }_{\alpha }^{(+)}$, but for dispersive materials the set of ${\omega }_{\alpha }^{(-)}$ values must be determined anew by solving the self-consistent eigenvalues with incoming boundary conditions (n.b. unlike the RSs, ${\omega }_{\alpha }^{(+)}$, which have only non-positive imaginary parts, $\mathrm{I}\mathrm{m}\left\{{\omega }_{\alpha }^{(-)}\right\}$ may be either positive or negative for dispersive scatterers).

The S-matrix can then be conveniently expressed in terms of eigenvectors of ${\omega }_{\alpha }^{(-)}$, and ${\omega }_{\alpha }^{(+)}$. For spherically symmetric particles, this leads to a particularly convenient, alternate expression for the response function residue, r_α , formulated entirely in terms of ${\omega }_{\alpha }^{(-)}$, and ${\omega }_{\alpha }^{(+)}$:

$$\begin{equation}\begin{array}{cc}{r}_{\alpha (q,n,\ell )}& =\frac{{r}_{\alpha (q,n,\ell )}^{(S)}}{2}=\frac{R}{2{c}_{b}}\frac{{\prod }_{{\ell }^{{{'}}}=-\infty }^{\infty }\left({\omega }_{\alpha (q,n,\ell )}^{(+)}-{\omega }_{\alpha (q,n,{\ell }^{{{'}}})}^{(-)}\right)}{{\prod }_{{\ell }^{{{'}}}=-\infty ;{\ell }^{{{'}}}\ne \ell }^{\infty }\left({\omega }_{\alpha (q,n,\ell )}^{(+)}-{\omega }_{\alpha (q,n,{\ell }^{{{'}}})}^{(+)}\right)}\\ & \simeq \frac{R}{2{c}_{\text{b}}}\frac{{\prod }_{{\ell }^{{{'}}}=-{L}_{\mathrm{max}}}^{{L}_{\mathrm{max}}}\left({\omega }_{\alpha (q,n,\ell )}^{(+)}-{\omega }_{\alpha (q,n,{\ell }^{{{'}}})}^{(-)}\right)}{{\prod }_{{\ell }^{{{'}}}=-{L}_{\mathrm{max}}\enspace ;{\ell }^{{{'}}}\ne \ell }^{{L}_{\mathrm{max}}}\left({\omega }_{\alpha (q,n,\ell )}^{(+)}-{\omega }_{\alpha (q,n,{\ell }^{{{'}}})}^{(+)}\right)},\end{array}\end{equation} \tag{ 54 }$$

where the poles (i.e. RS normalization as seen in equation (45)) can now be approximately determined as a function of a finite number of poles and zeros bounded by L_max (which we take to be the same as the truncation adopted in equation (43b)).

We numerically verified that for |ℓ| ≪ |L_max|, the approximate equality in equation (54) tends toward the analytic results of equation (45) for increasingly large values of L_max, but that for RSs near the cutoff (ℓ ∼ ±|L_max|), r_α differ significantly from the analytically calculated values. The advantage of this cutoff 'corrected' construction of the residues is that it enforces 'optimized' analytic behavior and static limits on the truncated Mie coefficients of equation (43b) at any truncation level, L_max, which thereby enables the RS expansion to converge to the Mie result. We plot the complex wavelengths of the first 25, positive wavelength, RSs of an R = 80 nm sphere in figure 10(a) for the simple Drude model of gold, equation (52b), and the more complex Drude–Lorentz model of equation (53), in figure 10(c). The results of the convergent RS expansion method with L_max = 25 are then plotted as red dashed curves in figures 10(b) and (d) respectively. This truncation level suffices to produce approximate curves that are almost indistinguishable from the exact (solid blue) curves. Similar agreements can be found for any multipolar order and thus allow Mie theory to be completely reconstructed within the RS framework for both dispersive and non-dispersive scatterers.

Zoom In Zoom Out Reset image size

Finally, it is crucial to remark that the convergence towards the Mie theory results was only possible because we had added the analytically determined non-resonant contribution to the RS spectral expansion in equation (43). Ignoring the non-resonant terms would again lead to catastrophic predictions (notably negative cross sections), as illustrated by the dotted green curves in figures 9(b) and (d) where we used the same convergent RS expansion as the red dashed curves, but where the non-resonant contributions were arbitrarily set to zero.

It is widely recognized that the plasmonic resonances can be viewed as a collective quasi-static excitation, so it therefore appears physically relevant that the non-resonant term, determined from quasi-static considerations, plays such a major role in metallic scattering. Such 'plasmonic' resonances are thereby rather different from their dielectric counterparts which tend to have their principal cross section peaks dominated by RS contributions (despite the continued importance of non-resonant contributions [17, 18, 22]).

Mie theory relies on developing fields in different regions in terms of VPWs, which are homogeneous media electromagnetic wave solutions in an angular momentum eigenstate basis. The excitation field, E _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),

$$\begin{equation}{\boldsymbol{E}}_{\text{e}}(k\boldsymbol{r})=E\sum\limits _{n=1}^{\infty }\sum\limits _{m=0}^{n}\left[{e}_{h,n,m}{\boldsymbol{M}}_{h,n,m}^{(1)}(k\boldsymbol{r})+{e}_{e,n,m}{\boldsymbol{N}}_{e,n,m}^{(1)}(k\boldsymbol{r})\right],\end{equation} \tag{ 84a }$$

where n is the electromagnetic angular momentum quantum number, m an absolute value of the angular momentum projection. The real parameter, E, determines the strength of the incident field. The real-valued multipole wave functions, N _q,n,m (k r ) and M _q,n,m (k r ), expressed in equation (67) are either of the magnetic, h(q = 0), type or electric, e(q = 1), type (denoted 'even' and 'odd' multipolar types in Bohren and Huffman [46]).

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ρ_ω using the notation of equation (11),

$$\begin{equation}{\boldsymbol{E}}_{\text{int}}(k\boldsymbol{r})=E\sum\limits _{n=1}^{\infty }\sum\limits _{m=0}^{n}\left[{s}_{h,n,m}{\boldsymbol{M}}_{h,n,m}^{(1)}(k{\rho }_{\omega }\boldsymbol{r})+{s}_{e,n,m}{\boldsymbol{N}}_{e,n,m}^{(1)}(k{\rho }_{\omega }\boldsymbol{r})\right].\end{equation} \tag{ 84b }$$

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,

$$\begin{equation}{\boldsymbol{E}}_{\text{s}}(k\boldsymbol{r})=E\sum\limits _{n=1}^{\infty }\sum\limits _{m=0}^{n}\left[{f}_{h,n,m}{\boldsymbol{M}}_{h,n,m}^{(+)}(k\boldsymbol{r})+{f}_{e,n,m}{\boldsymbol{N}}_{e,n,m}^{(+)}(k\boldsymbol{r})\right],\end{equation} \tag{ 84c }$$

which can be viewed as replacing the spherical Bessel functions in the definitions of the regular VPWs by outgoing Hankel functions as shown in equation (68).

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 Ω, such that, f = T.e and s = Ω.e. The fact that scattering by spheres does not change the angular momentum, both T and Ω become diagonal matrices for spherically symmetric scatterers and traditional 'Mie' theory provides algebraic expressions for these diagonal elements, henceforth denoted T_q,n and Ω_q,n . For the different multipole orders, one has,

$$\begin{equation}\begin{aligned}{f}_{q,n,m}& ={T}_{q,n}{e}_{q,n,m}=-\frac{{N}_{q,n}(\omega ,R)}{{D}_{q,n}(\omega ,R)}{e}_{q,n,m}\\ & {b}_{n}=-{T}_{h,n},\quad {a}_{n}=-{T}_{e,n},\end{aligned}\end{equation} \tag{ 85a }$$

where we specified the relation between T_q,n and the time honored Mie coefficients, traditionally denoted a_n and b_n . The 'numerator' and 'denominator' functions in equation (85a) are functions of the angular frequency, ω and particle radius, R are written out explicitly in equations (86) and (87).

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,

$$\begin{equation}\begin{aligned}{s}_{q,n,m}& \equiv {{\Omega}}_{q,n}{e}_{q,n,m}=\frac{1}{{D}_{q,n}(\omega ,R)}{e}_{q,n,m}\\ & \quad {c}_{n}={{\Omega}}_{h,n},\quad {d}_{n}={{\Omega}}_{e,n}.\end{aligned}\end{equation} \tag{ 85b }$$

The 'numerator' and 'denominator' functions in equation (85) are functions of the angular frequency, ω and particle radius, R:

$$\begin{equation}{N}_{h,n}(\omega ,R)\equiv z\frac{{\mu }_{\omega }{j}_{n}({\rho }_{\omega }z){\psi }_{n}^{\prime }(z)-{\psi }_{n}^{\prime }({\rho }_{\omega }z){j}_{n}(z)}{\text{i}{\mu }_{\omega }}\end{equation} \tag{ 86a }$$

$$\begin{equation}{N}_{e,n}(\omega ,R)\equiv z\frac{{\varepsilon }_{\omega }{j}_{n}({\rho }_{\omega }z){\psi }_{n}^{\prime }(z)-{\psi }_{n}^{\prime }({\rho }_{\omega }z){j}_{n}(z)}{\text{i}{\rho }_{\omega }},\end{equation} \tag{ 86b }$$

and,

$$\begin{equation}{D}_{h,n}(\omega ,R)\equiv z\frac{{\mu }_{\omega }{j}_{n}({\rho }_{\omega }z){\xi }_{n}^{\prime }(z)-{\psi }_{n}^{\prime }({\rho }_{\omega }z){h}_{n}(z)}{\text{i}{\mu }_{\omega }}\end{equation} \tag{ 87a }$$

$$\begin{equation}{D}_{e,n}(\omega ,R)\equiv z\frac{{\varepsilon }_{\omega }{j}_{n}({\rho }_{\omega }z){\xi }_{n}^{\prime }(z)-{\psi }_{n}^{\prime }({\rho }_{\omega }z){h}_{n}(z)}{\text{i}{\rho }_{\omega }},\end{equation} \tag{ 87b }$$

where we have again used z = kR = ωR/c and used the condensed notation of equation (11) for ɛ_ω , μ_ω , and ρ_ω , for the (possibly frequency dependent) constitutive parameters: an inspection of equations (86) and (87) 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 figure 6 are obtained from the Mie coefficients of equation (85a) via well-known relations:

$$\begin{equation}{Q}_{\text{ext}}\equiv \frac{{\sigma }_{\text{ext}}}{\pi {R}^{2}}=\frac{1}{\pi {R}^{2}}\sum\limits _{n=1}^{\infty }\left({\sigma }_{\text{ext},n}^{(h)}+{\sigma }_{\text{ext},n}^{(e)}\right)=-\frac{2}{{(kR)}^{2}}\sum\limits _{n=1}^{\infty }(2n+1)\mathrm{R}\mathrm{e}\left[{T}_{h,n}+{T}_{e,n}\right]\end{equation} \tag{ 88a }$$

$$\begin{equation}{Q}_{\text{scat}}=\frac{2}{{(kR)}^{2}}\sum\limits _{n=1}^{\infty }(2n+1)\left[{\left\vert {T}_{h,n}\right\vert }^{2}+{\left\vert {T}_{e,n}\right\vert }^{2}\right]\end{equation} \tag{ 88b }$$

$$\begin{equation}{Q}_{\text{abs}}={Q}_{\text{ext}}-{Q}_{\text{scat}}.\end{equation} \tag{ 88c }$$

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

$$\begin{equation}{z}_{\alpha (h,n,m,\ell )}={\int }_{{\mathcal{V}}_{\infty }}\mathrm{d}\mathbf{r}\left\{{\left.\frac{\mathrm{d}\left[\omega \varepsilon (\boldsymbol{r},\omega )\right]}{\mathrm{d}\omega }\right\vert }_{{\omega }_{\alpha }}{\boldsymbol{E}}_{\alpha }^{2}(\boldsymbol{r})-{\left.\frac{\mathrm{d}\left[\omega \mu (\boldsymbol{r},\omega )\right]}{\mathrm{d}\omega }\right\vert }_{{\omega }_{\alpha }}{\boldsymbol{H}}_{\alpha }^{2}(\boldsymbol{r})\right\}.\end{equation} \tag{ 91 }$$

For clarity, the overall magnetic multipole mode sign factor, (−1)ⁿ , resulting from the iⁿ phase factor in multipole fields (cf equations (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, ɛ(ω), 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,

$$\begin{align}& {\left.\frac{\mathrm{d}\left[\omega \varepsilon (\omega )\right]}{\mathrm{d}\omega }\right\vert }_{{\omega }_{\alpha }}{\int }_{0}^{R}\mathrm{d}\mathbf{r}\enspace {\boldsymbol{E}}_{\alpha (h,n,m,\ell )}^{2}(\boldsymbol{r})\\ & =\frac{{z}_{\alpha }^{3}{\gamma }_{\alpha (h,n)}^{2}}{{\mathcal{N}}_{\alpha }^{2}}{\left.\frac{\mathrm{d}\left[\omega \varepsilon (\omega )\right]}{\mathrm{d}\omega }\right\vert }_{{\omega }_{\alpha }}{\int }_{0}^{1}{~{r}}^{2}{j}_{n}^{2}({\rho }_{\alpha }{z}_{\alpha }~{r})\mathrm{d}~{r}=\frac{{z}_{\alpha }{\gamma }_{\alpha (h,n)}^{2}}{2{\mathcal{N}}_{\alpha }^{2}}{\left.\frac{\mathrm{d}\left[\omega \varepsilon (\omega )\right]}{\mathrm{d}\omega }\right\vert }_{{\omega }_{\alpha }}\\ & \quad {\times}\frac{{\left[{\psi }_{n}^{\prime }({\rho }_{\alpha }{z}_{\alpha })\right]}^{2}+{\psi }_{n}^{2}({\rho }_{\alpha }{z}_{\alpha })-n(n+1){j}_{n}^{2}({\rho }_{\alpha }{z}_{\alpha })-{j}_{n}({\rho }_{\alpha }{z}_{\alpha }){\psi }_{n}^{\prime }({\rho }_{\alpha }{z}_{\alpha })}{{\varepsilon }_{\alpha }{\mu }_{\alpha }}\\ & =\frac{{z}_{\alpha }}{2{\mathcal{N}}_{\alpha }^{2}}{\left.\frac{\mathrm{d}\left[\omega \varepsilon (\omega )\right]}{\mathrm{d}\omega }\right\vert }_{{\omega }_{\alpha }}\frac{{{\Xi}}_{h,n}^{(-)}}{{\varepsilon }_{\alpha }}=\frac{{z}_{\alpha }}{2}\frac{{{\Xi}}_{h,n}^{(-)}+{{\Xi}}_{h,n}^{(-)}{\omega }_{\alpha }\frac{\mathrm{d}}{\mathrm{d}\omega }{\left.\enspace \mathrm{ln}\enspace \varepsilon (\omega )\right\vert }_{{\omega }_{\alpha }}}{{\mathcal{N}}_{\alpha }^{2}},\end{align} \tag{ 92a }$$

where we used equation (14b) for γ_α(n,h) and the expression of equation (37c) for ${{\Xi}}_{n}^{(h,{\pm})}$.

Denoting by μ(ω), the magnetic permeability of the sphere divided by permeability of the external medium, the magnetic field integral of equation (91) is,

$$\begin{align}& \quad -{\left.\frac{\mathrm{d}\left[\omega \mu (\omega )\right]}{\mathrm{d}\omega }\right\vert }_{{\omega }_{\alpha }}{\int }_{0}^{R}\mathrm{d}\mathbf{r}\enspace {\boldsymbol{H}}_{\alpha (h,n,m,\ell )}^{2}(\boldsymbol{r})\\ & ={\left.\frac{\mathrm{d}\left[\omega \mu (\omega )\right]}{\mathrm{d}\omega }\right\vert }_{{\omega }_{\alpha }}\frac{{z}_{\alpha }{\gamma }_{\alpha (h,n)}^{2}}{{\mathcal{N}}_{\alpha }^{2}}\frac{{\varepsilon }_{\alpha }}{{\mu }_{\alpha }}{\int }_{0}^{1}\frac{n(n+1){j}_{n}^{2}({\rho }_{\alpha }{z}_{\alpha }~{r})+{\left[{\psi }_{n}^{\prime }({\rho }_{\alpha }{z}_{\alpha }~{r})\right]}^{2}}{{\rho }_{\alpha }^{2}}\mathrm{d}~{r}\\ & ={\left.\frac{\mathrm{d}\left[\omega \mu (\omega )\right]}{\mathrm{d}\omega }\right\vert }_{{\omega }_{\alpha }}\frac{{z}_{\alpha }{\gamma }_{\alpha (h,n)}^{2}}{{\mathcal{N}}_{\alpha }^{2}}\\ & \quad {\times}\frac{{\varepsilon }_{\alpha }}{{\mu }_{\alpha }}\frac{{\left[{\psi }_{n}^{\prime }({\rho }_{\alpha }{z}_{\alpha })\right]}^{2}+{\psi }_{n}^{2}({\rho }_{\alpha }{z}_{\alpha })-n(n+1){j}_{n}^{2}({\rho }_{\alpha }{z}_{\alpha })+{j}_{n}({\rho }_{\alpha }{z}_{\alpha }){\psi }_{n}^{\prime }({\rho }_{\alpha }{z}_{\alpha })}{2{\varepsilon }_{\alpha }{\mu }_{\alpha }}\\ & =\frac{{z}_{\alpha }}{2{\mathcal{N}}_{\alpha }^{2}}{\left.\frac{\mathrm{d}\left[\omega \mu (\omega )\right]}{\mathrm{d}\omega }\right\vert }_{{\omega }_{\alpha }}\frac{{{\Xi}}_{h,n}^{(+)}}{{\mu }_{\alpha }}=\frac{{z}_{\alpha }}{2}\frac{{{\Xi}}_{h,n}^{(+)}+{{\Xi}}_{h,n}^{(+)}{\omega }_{\alpha }{\left.\frac{\mathrm{d}}{\mathrm{d}\omega }\enspace \mathrm{ln}\enspace \mu (\omega )\right\vert }_{{\omega }_{\alpha }}}{{\mathcal{N}}_{\alpha }^{2}},\end{align} \tag{ 92b }$$

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

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

$$\begin{align}& {\int }_{R}^{\infty }\mathrm{d}\mathbf{r}\enspace {\boldsymbol{E}}_{\alpha (h,n,m,\ell )}^{2}(\boldsymbol{r})\rightarrow {z}_{\alpha }^{3}\enspace \underset{\eta \to 0}{\mathrm{lim}}{\int }_{1}^{\infty }{~{r}}^{2}{h}_{n}^{2}({z}_{\alpha }~{r}){\text{e}}^{-\eta {~{r}}^{2}}\enspace \mathrm{d}~{r}\\ & =\frac{{z}_{\alpha }}{2}\frac{n(n+1){h}_{n}^{2}({z}_{\alpha })-{\left[{\xi }_{n}^{\prime }({z}_{\alpha })\right]}^{2}-{\xi }_{n}^{2}({z}_{\alpha })+{h}_{n}({z}_{\alpha }){\xi }_{n}^{\prime }({z}_{\alpha })}{{\mathcal{N}}_{\alpha }^{2}},\end{align} \tag{ 93a }$$

and,

$$\begin{align}& \quad -{\int }_{R}^{\infty }\mathrm{d}\mathbf{r}\enspace {\boldsymbol{H}}_{\alpha (h,n,m,\ell )}^{2}(\boldsymbol{r})\rightarrow {z}_{\alpha }\enspace \underset{\eta \to 0}{\mathrm{lim}}{\int }_{1}^{\infty }\\ & \quad {\times}\left\{n(n+1){h}_{n}^{2}({z}_{\alpha }~{r})+{\left[{\xi }_{n}^{\prime }({z}_{\alpha }~{r})\right]}^{2}\right\}{\text{e}}^{-\eta {~{r}}^{2}}\enspace \mathrm{d}~{r}\\ & =\frac{{z}_{\alpha }}{2}\frac{n(n+1){h}_{n}^{2}({z}_{\alpha })-{\left[{\xi }_{n}^{\prime }({z}_{\alpha })\right]}^{2}-{\xi }_{n}^{2}({z}_{\alpha })-{h}_{n}({z}_{\alpha }){\xi }_{n}^{\prime }({z}_{\alpha })}{{\mathcal{N}}_{\alpha }^{2}},\end{align} \tag{ 93b }$$

and the sum of the above two integrals yields:

$$\begin{equation}\underset{\eta \to 0}{\mathrm{lim}}{\int }_{R}^{\infty }{\text{e}}^{-\eta {~{r}}^{2}}\enspace \mathrm{d}\mathbf{r}\left\{{\boldsymbol{E}}_{\alpha }^{2}(\boldsymbol{r})-{\boldsymbol{H}}_{\alpha }^{2}(\boldsymbol{r})\right\}={z}_{\alpha }\frac{n(n+1){h}_{n}^{2}({z}_{\alpha })-{\left[{\xi }_{n}^{\prime }({z}_{\alpha })\right]}^{2}-{\xi }_{n}^{2}({z}_{\alpha })}{{\mathcal{N}}_{\alpha }^{2}}.\end{equation} \tag{ 94 }$$

Finally, putting together the results of equations (92a), (92b) and (94) into equation (91), one finds the normalization factor of equation (37a) for magnetic type RSs for spherical scatterers with full temporal dispersion.

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

$$\begin{equation}\begin{aligned}{z}_{\alpha (e,n,m,\ell )}={\int }_{{\mathcal{V}}_{\infty }}\mathrm{d}\boldsymbol{r}\left\{{\left.\frac{\mathrm{d}\left[\omega \varepsilon (\boldsymbol{r},\omega )\right]}{\mathrm{d}\omega }\right\vert }_{{\omega }_{\alpha }}{\boldsymbol{E}}_{\alpha }^{2}(\boldsymbol{r})-{\left.\frac{\mathrm{d}\left[\omega \mu \left(\boldsymbol{r},\omega \right)\right]}{\mathrm{d}\omega }\right\vert }_{{\omega }_{\alpha }}{\boldsymbol{H}}_{\alpha }^{2}(\boldsymbol{r})\right\},\end{aligned}\end{equation} \tag{ 95 }$$

where the overall electric multipole mode sign factor, (−1)ⁿ⁻¹, will be henceforth suppressed for simplicity like we did above for the magnetic modes above.

The E-field integration inside the sphere is,

$$\begin{align}& {\left.\frac{\mathrm{d}\left[\omega \varepsilon (\omega )\right]}{\mathrm{d}\omega }\right\vert }_{{\omega }_{\alpha }}{\int }_{0}^{R}\mathrm{d}\boldsymbol{r}\enspace {\boldsymbol{E}}_{\alpha (e,n,m,\ell )}^{2}(\boldsymbol{r})\\ & =\frac{{z}_{\alpha }{\gamma }_{\alpha (e,n)}^{2}}{{\mathcal{N}}_{\alpha }^{2}}{\left.\frac{\mathrm{d}\left[\omega \varepsilon (\omega )\right]}{\mathrm{d}\omega }\right\vert }_{{\omega }_{\alpha }}{\int }_{0}^{1}\frac{n(n+1){j}_{n}^{2}({\rho }_{\alpha }{z}_{\alpha }r)+{\left[{\psi }_{n}^{\prime }({\rho }_{\alpha }{z}_{\alpha }r)\right]}^{2}}{{\rho }_{\alpha }^{2}}\mathrm{d}r\\ & =\frac{{z}_{\alpha }{\gamma }_{\alpha (e,n)}^{2}}{{\mathcal{N}}_{\alpha }^{2}}{\left.\frac{\mathrm{d}\left[\omega \varepsilon (\omega )\right]}{\mathrm{d}\omega }\right\vert }_{{\omega }_{\alpha }}\\ & \quad {\times}\frac{{\left[{\psi }_{n}^{\prime }({\rho }_{\alpha }{z}_{\alpha })\right]}^{2}+{\psi }_{n}^{2}({\rho }_{\alpha }{z}_{\alpha })-n(n+1){j}_{n}^{2}({\rho }_{\alpha }{z}_{\alpha })+{j}_{n}({\rho }_{\alpha }{z}_{\alpha }){\psi }_{n}^{\prime }({\rho }_{\alpha }{z}_{\alpha })}{2{\varepsilon }_{\alpha }{\mu }_{\alpha }}\\ & =\frac{{z}_{\alpha }}{2{\mathcal{N}}_{\alpha }^{2}}\frac{{{\Xi}}_{n}^{(e,+)}}{{\varepsilon }_{\alpha }}{\left.\frac{\mathrm{d}\left[\omega \varepsilon (\omega )\right]}{\mathrm{d}\omega }\right\vert }_{{\omega }_{\alpha }}=\frac{{z}_{\alpha }}{2}\frac{{{\Xi}}_{n}^{(e,+)}+{{\Xi}}_{n}^{(e,+)}{\omega }_{\alpha }{\left.\frac{\mathrm{d}}{\mathrm{d}\omega }\enspace \mathrm{ln}\enspace {\varepsilon }_{\alpha }(\omega )\right\vert }_{{\omega }_{\alpha }}}{{\mathcal{N}}_{\alpha }^{2}},\end{align} \tag{ 96a }$$

where we used equation (14a) for γ_e,n,α and used the definition of equation (37d) for ${{\Xi}}_{n}^{(e,{\pm})}$.

The H-field integral inside the sphere is,

$$\begin{align}& \quad -{\left.\frac{\mathrm{d}\left[\omega \mu (\omega )\right]}{\mathrm{d}\omega }\right\vert }_{{\omega }_{\alpha }}{\int }_{0}^{R}\mathrm{d}\boldsymbol{r}\enspace {\boldsymbol{H}}_{\alpha (e,n,m,\ell )}^{2}(\boldsymbol{r})\\ & ={\left.\frac{\mathrm{d}\left[\omega \mu (\omega )\right]}{\mathrm{d}\omega }\right\vert }_{{\omega }_{\alpha }}\frac{{z}_{\alpha }^{3}{\gamma }_{\alpha (e,n)}^{2}}{{\mathcal{N}}_{\alpha }^{2}}\frac{{\varepsilon }_{\alpha }}{{\mu }_{\alpha }}{\int }_{0}^{1}{~{r}}^{2}{j}_{n}^{2}({\rho }_{\alpha }{z}_{\alpha }~{r})\mathrm{d}~{r}\\ & =\frac{{z}_{\alpha }{\gamma }_{\alpha (e,n)}^{2}}{{\mathcal{N}}_{\alpha }^{2}}\frac{1}{{\mu }_{\alpha }}{\left.\frac{\mathrm{d}\left[\omega \mu (\omega )\right]}{\mathrm{d}\omega }\right\vert }_{{\omega }_{\alpha }}{\varepsilon }_{\alpha }\\ & \quad {\times}\frac{{\left[{\psi }_{n}^{\prime }({\rho }_{\alpha }{z}_{\alpha })\right]}^{2}+{\psi }_{n}^{2}({\rho }_{\alpha }{z}_{\alpha })-n(n+1){j}_{n}^{2}({\rho }_{\alpha }{z}_{\alpha })-{j}_{n}({\rho }_{\alpha }{z}_{\alpha }){\psi }_{n}^{\prime }({\rho }_{\alpha }{z}_{\alpha })}{2{\varepsilon }_{\alpha }{\mu }_{\alpha }}\\ & =\frac{{z}_{\alpha }}{2{\mathcal{N}}_{\alpha }^{2}}\frac{{{\Xi}}_{n}^{(e,-)}}{{\mu }_{\alpha }}{\left.\frac{\mathrm{d}\left[\omega \mu (\omega )\right]}{\mathrm{d}\omega }\right\vert }_{{\omega }_{\alpha }}=\frac{{z}_{\alpha }}{2}\frac{{{\Xi}}_{n}^{(e,-)}+{\omega }_{\alpha }{{\Xi}}_{n}^{(e,-)}{\left.\frac{\mathrm{d}}{\mathrm{d}\omega }\enspace \mathrm{ln}\enspace {\mu }_{\text{s}}(\omega )\right\vert }_{{\omega }_{\alpha }}}{{\mathcal{N}}_{\alpha }^{2}},\end{align} \tag{ 96b }$$

where ${{\Xi}}_{e,n}^{(-)}$ is defined in equation (37d).

The integrals outside the sphere are identical to those of equation (93) with the roles of electric and magnetic fields reversed, so we obtain the same result as equation (94) for the field integrals in a region exterior to a sphere of radius R. Inserting the results of equations (94), (96a) and (96b) into equation (95), one finds the normalization factor of equation (37b) for electric type RSs for spherical scatterers with full temporal dispersion.