Eigenstate normalization for open and dispersive systems

# Eigenstate normalization for open and dispersive systems

B. Stout, R. Colom, N. Bonod and R.C. McPhedran,
Aix-Marseille Univ, CNRS, Centrale Marseille, Institut Fresnel, 13397 Marseille, France.
Avignon Université, UMR 1114 EMMAH, Avignon Cedex 84018, France,
IPOS, School of Physics, University of Sydney, 2006, Australia
###### Abstract

Resonant States, also known as ‘Quasi-Normal’ Modes, are attractive basis sets for the spectral expansions of open radiative systems. Nevertheless, the exponential divergence of the resonant states in the far-field has long been an obstacle for defining viable scalar products. This theoretical work shows that the apparent ‘divergences’ of the resonant state scalar product integrals can be resolved analytically and that this leads to orthogonalization of the resonant states for dispersionless scatterers and finite normalizations for the resonant states of spherical scatterers with electric and/or magnetic temporal dispersion. The correctly normalized expansions are shown to decouple the power of the iconic ‘Mie’ theory by re-expressing the exact Mie coefficients as meromorphic(spectral) functions of frequency. The methods and formulas demonstrated here can be generalized to more complex geometries, and provide new physical insight into the nature of light-matter interactions.

## 1 Introduction

The eigenstates of open systems, known as Resonant States (RSs), can be viewed as a means of simultaneously describing both the spatial and temporal behavior of the oscillatory ‘modes’ of an open system. As such, RSs have been a highly fruitful, albeit problematic, concept in many branches of physics for over a hundred years, [1] where they have been called by a wide variety of names like: quasi-normal modes, transient modes, leaky waves, etc.. The RS concept was widely developed in quantum mechanics [2] after its pioneering success by Gamow for describing alpha decay [3], and later by Sieggert for characterizing nuclear reactions [4]. Throughout the 20 century, quantum RS descriptions continued to draw the attention of prominent physicists like Wigner [5], Peierls [6], and Zel’dovich [7, 8] who described the mathematical difficulties associated with RSs as follows:

An exponentially decaying state that describes, for example, the phenomenon of a decay, is characterized by a complex value of the energy, the imaginary part of the energy giving the decay probability. The wave function of this state increases exponentially in absolute value at large distances, and therefore the usual methods of normalization, of perturbation theory, and of expansion in terms of eigenfunctions do not apply to this state.

The utility of RSs is not limited to the quantum realm however, and they are equally applicable to a host of ‘classical field’ systems existing in 1, 2, or 3 spatial dimensions and which are governed by a wide variety of differential equations. It thus appears that RSs are a nearly universal tool for describing open systems in essentially arbitrary spatial dimensions, but we will henceforth restrict attention to 3D systems for the purpose of discussion.

Combining the assumption that energy is ‘radiated’ into a ‘background media’ in a form governed by a Helmholtz type equation with outgoing boundary conditions (to describe openness), leads to frequency domain RSs being characterized by complex eigenvalue wavenumbers, , associated with a far-field radial dependence proportional to , where is the phase-velocity of the radiation in the background media. The imaginary part of the RS wavenumber, , is negative due to causality considerations, which leads to a desirable exponential decay of an excitation in the time domain, but an inspection of the signs in the far-field behavior shows that temporal decay is necessarily accompanied by an exponential increase in the RS amplitude at large distances, . An example of the generality of the RS approach is found in electromagnetic scattering problems [9, 10, 11], where one encounters both vectorial and scalar RSs problems. Mathematical descriptions of vector RSs are more complicated than those for scalar fields, and this was a major motivation for this work (as well as ref. [12]: wherein the necessary mathematical groundwork is developed in detail but which is not currently intended for publication in journal form).

This spatial divergence of the RSs doesn’t invalidate their utility, since causality considerations regularize the exponential divergence when one transforms from frequency domain descriptions back into the time domain (cf. [1, 13, 14, 15, 16]). Nevertheless, such considerations fail to alleviate the mathematical and conceptual difficulties posed by the ‘exponential catastrophe’ [17] of spectral theories in the frequency (energy) domain. Although a variety of papers have studied and more or less ‘resolved’ this problem, either analytically or numerically, by use of the method of Perfectly Matched Layers (PML) to tame the ‘QNM’ scalar product integrals [18, 19, 20, 21], a general, conceptually simple, and mathematically rigorous approach still appeared to be lacking. This work helps to fill this gap with a general, mathematically rigorous, ‘benchmark’ approach to RS normalization/orthogonalization.

It is known since at least the work of Zel’dovich [7], that scalar products of scalar RS theories should be defined without complex conjugating the adjoint, ‘bra’, fields (cf. section 2.2, and refs. [8, 17, 22, 23]). With this convention, even though the scalar product integrand amplitude still diverges, it also oscillates as , thus rendering it possible to regularize the scalar product integrals to obtain finite values. Zel’dovich [7] proposed to regularize the scalar products by first introducing a Gaussian ‘killing’ function, into the scalar product integral and then taking the limit after the integrations. Analytic expressions for such integrals were studied in a 1992 paper by one of the present authors and two colleagues [24].

Recent extensions of these techniques for 3D electromagnetic scattering applications were referred to somewhat picturesquely (or melodically) as “Killing Mie Softly” in an arXiv paper detailing the distribution theory mathematics by two of the present authors [12]. The present article demonstrates that the implementation of these analytical methods leads to RS scalar products that are unambiguously defined in terms of spatial integrals over all space. This allows spectral expansions of open systems to access to the entire panoply of methods that have encountered such sweeping success in the description of closed (Hermitian) systems. A hint of the potential power of this technique, is provided by retrieving and extending the results of the iconic ‘Mie’ theory for light scattering.

This article is organized as follows: Section 2 introduces the RS concept in the context of electromagnetic field equations. It proceeds to provide analytic expressions for the RSs of isotropic spherical scatterers. A new result of section 2.2 is to give an example of how the results of ref. [12] confirm RS orthogonality for dispersion-free spherical scatterers. Subsection 2.3 presents one of the principal derivations of this paper, which are analytic formulas for the RS normalizations for the case of spheres composed of media with full electric and magnetic dispersive constitutive parameters. This section also strives to give some intuitive understanding of the mathematical regularization techniques by means of graphic plots.

Section 3 translates the RS eigenvalues and normalizations into meromorphic expansions of the Mie theory coefficients. This section is important since it shows how the RS formulation can reproduce exact Mie theory results up to essentially arbitrary accuracy. Secondly, one sees in this section that RS expansions yield meromorphic expressions valid for all frequencies, whereas traditional Mie theory formulas must be reevaluated for each frequency under consideration. Section 3 ends with a brief ‘user’s’ guide on how to implement the results of this work in the case of spherical scatterers. This is intended to highlight the simplicity of the final RS physics as well as serve as a ‘guide’ for those chiefly interested in applications.

Finally, section 4 gives a number of numerical results in the form of two tables for the RS eigenvalues together with their normalization factors for different dispersion models from the optical literature (including Drude type dispersion models and its generalizations). These can serve as benchmark values for those using these methods or elaborating purely numerical models for nanophotonic and metamaterial problems. This section is intended to give a brief glimpse into the rich physical insight that RSs can bring to the analysis of complex physical processes. Section 5 constitutes a brief conclusion to this paper. Appendices B-E, elaborate on mathematical definitions and derivations presented in the main text.

## 2 Resonant states

When formulating Resonant States (RSs) in electromagnetism, it is convenient to express the electromagnetic field ‘ket’ state, and appropriately dimensioned source currents, , as multi-component fields (6 components),

 (1)

where and are respectively proportional to the electric and magnetic fields with identical units, here m, as this is a convenient choice for treatments of LDOS [25], energy density, Green’s functions, etc.

Material response is modeled as finite size ‘scatterers’ with sharp boundaries and composed of materials with spatially local but temporally dispersive constitutive parameters, , and/or , immersed in an isotropic non-dispersive background medium described by real-valued constitutive parameters or , with the background media wave velocity given by , where is the vacuum speed of light. Henceforth, all scatterer constitutive parameters will be normalized with respect to the background medium properties, , and .

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

 ωcbΓ(ω)|Ψ⟩=L|Ψ⟩+1i|J⟩, (2a) where the medium metric, Γ(ω), and the linear differential operator, L, are symmetric 6×6 matrix operators, Γ(ω)≡[ε(r,ω)00−μ(r,ω)]andL≡i[0∇×∇×0], (2b)

Resonant states are defined as source-free eigensolutions of Eq.(2) obeying outgoing boundary conditions,

 kαΓ(ωα)|Ψα⟩=L|Ψα⟩, (3a) where kα=ωα/cb is the complex eigenvalue wave-number in the background medium, kα≡k′α+ik′′α. (3b)

We point out that the medium metric, , is not entirely an arbitrary function of frequency since we require it to obey Kramers-Kronig causality relations, with the real-valued nature of the fields in the time domain imposing, . This has the consequence that even though negative frequency RSs must be included in our analysis, they are not independent of the positive frequency solutions. One should remark that Eq.(3a) is a self-consistent eigenvalue equation in terms of the complex wavenumber (frequency), which explains why finding the RS eigenvalues is a procedure that must almost always be carried out numerically.

It is also important to remark that RSs occur for both positive and negative real frequency parts. One readily verifies however that these RSs are not independent. Notably, taking the complex conjugate of Eq.(3) for a RS with positive real part () and remarking that , one obtains,

 (4)

which is of the same form as Eq.(3), so for any RS eigenvalue, , there is an associated RS with eigenvalue, . The RSs with negative real parts should not be confused with incoming field solutions as they are still associated with outgoing fields. The presence of these RSs (mirror symmetric with respect to the imaginary axis) insures that fields are real valued in the time-domain, but they can also be viewed as a consequence of time reversal symmetry. In view of the above, it is convenient to designate RSs with positive (negative) real parts with positive (negative) integer indices, , which are related via, . In the cases considered here, there will generally be at most one RS eigenvalue on the negative imaginary axis for a given multipole order, so we can denote such states when they occur (even though multiple RSs on the imaginary axis are not excluded in general).

For scatterers of arbitrary forms, there is little choice but to solve the eigenvalue equation of Eq.(3) numerically. However, given that the equations of motion in Eq.(2) are separable in spherical coordinates, one can obtain analytical expressions for the RSs for an isotropic sphere in spherical coordinates. These analytical expressions for the RS of spherical scatterers are established in the next section.

### 2.1 Resonant states for spherically symmetric scatterers

In the region outside an isotropic spherical scatterer, the resonant states of a given multipolar order are required to be proportional to the standard vector spherical partial waves (VPWs) described in classic texts, traditionally denoted and , (and written out explicitly in section B.3). For the sake of clarity, the radial distance will henceforth be normalized with respect to the sphere radius (i.e. sphere surface at ). Using these conventions and notations, the RS kets respectively of the electric-type, and magnetic-type, , can be expressed in the region outside the sphere, , as:

 (5a) where the kα is the complex RS wavenumber of Eq.(3), and zα≡kαR, the corresponding RS size parameter, with M(+)(e,o),n,m and N(+)(e,o),n,m, the dimensionless ‘even’ and ‘odd’ vector partial waves (VPWs) that satisfy outgoing boundary conditions (cf. section B.3 and ref. [26]). The k3/2α factor in the RS wave-functions are a direct consequence of choosing field dimensions of m−3/2, and is typical of field theoretic approaches.

For the region inside the sphere, , one can analogously deduce that the electric-type and magnetic-type RSs inside the sphere are proportional to ket expressions in terms of the regular VPWs, and with arguments scaled in terms of the wavenumber of the internal medium:

 (5b)

where and are the normalized constitutive parameters evaluated at the resonant state frequency, , and corresponds to a normalized refractive index, also evaluated at the RS frequency, i.e.,

 εα≡εs(ωα)εb,μα≡μs(ωα)μb,ρα=√εs(ωα)μs(ωα)εbμb, (5c)

where we are considering the most general case of dispersive media. For non-dispersive media, all three of these constitutive parameters reduce to frequency independent real-valued numbers. Concerning notation, we remark that the spherical symmetry partially diagonalizes the RS basis in terms of their angular momentum and electric/magnetic natures. So in principal each discrete index is determined by the discrete indices, , and , as well as an additional discrete index, ranging from to , say for an electric type RS, but in the interest of simplicity, we here and henceforth simply write .

The complex multiplicative factors and in Eq.(5b) must be chosen in order to satisfy the two boundary conditions associated with the continuity of the components of the electric and magnetic fields that are parallel to the sphere surface, given the choice of Eq.(5a) for the external outgoing fields. These boundary conditions impose that each factor must simultaneously satisfy two distinct conditions:

 γe,n,α =ραεαhn(zα)jn(ραzα)=ραξ′n(zα)ψ′n(ραzα) (6a) γh,n,α =hn(zα)jn(ραzα)=μαξ′n(zα)ψ′n(ραzα), (6b)

where are spherical Bessel functions, while are the outgoing spherical Hankel functions with and their respective Ricatti-Bessel function counterparts (cf. section B.3). For arbitrary size parameters and frequencies, the respective expressions for each in Eq.(6) are not equal, and their equality defines the occurrence of a resonant state size parameter, . The equalities in Eq.(6) also correspond to poles of the electric and magnetic Mie coefficients recalled in Eq.(33) of Appendix B.

It is important to remark that the RSs of Eq.(2.1) are only determined up to the (complex) numerical proportionality factor, , whose value needs to be determined. Analytic expressions for the are presented in Eq.(14) below and derived using the analytical techniques developed in ref. [12] and illustrated in appendix A. This non-trivial normalization is one of the principal results of this work.

### 2.2 Scalar products and orthogonality of resonant states

As already mentioned in the introduction, we make the usual treatment of the RS adjoint space by defining a ‘bra’ state, , without complex conjugation of the fields. With this convention, the scalar product is simply,

 ⟨Ψα|Ψβ⟩≡∫∞dr{eα(r)⋅eβ(r)+hα(r)⋅hβ(r)}. (7)

where the integral is carried out over all space. Given the structure of the field equations, in Eq.(2), the general orthogonality conditions for the RSs require using the medium metric, , as an ‘adjoint’ operator, like that commonly used in relativistic field theory so that orthogonality for a dispersionless medium takes the form of a weighted scalar product, , which written out explicitly yields,

 (8)

and the first line serves as a reminder that the matrix can be here associated with either the ‘bra’ or ‘ket’ state. The electric and magnetic field components are not independent, and orthogonality is then applicable to either electric or magnetic fields for dispersion free media and .

The example of the RS orthogonality for the RSs of dispersionless spherical scatterers is found by integrating the magnetic part of the RS scalar product for two electric type RSs in Eq.(2.1) (type with ). The volume scalar product integral inside the sphere is,

 ∫r<1drμh(e)n,m,α(r)⋅h(e)n,m,β(r)=−z3/2αz3/2βγe,n,αγe,n,βN(e)n,αN(e)n,βε{∫10r2jn(ρzαr)jn(ρzβr)dr} (9a) =−z3/2αz3/2βγe,n,αγe,n,βN(e)n,αN(e)n,βε{jn(ρzα)ψ′n(ρzβ)−jn(ρzβ)ψ′n(ρzα)ρ2z2α−ρ2z2β} (9b) =−z3/2αz3/2βN(e)n,αN(e)n,βzβhn(zα)h′n(zβ)−zαhn(zβ)h′n(zα)z2α−z2β, (9c) where Eq.(9a) is obtained after analytical integration of the angular variables, while term in brackets of (9b) is a direct application of the bracketed integral in Eq.(9a) as given by Eq.(85a) of [12] (cf. also Watson ref. [27], and Eqs.(28) and (31b)). The final result of Eq.(9c) exploits the expressions of Eq.(6a) for γe,n,α, and is only true provided that both zα and zβ are RS (size parameter) eigenvalues.

One completes the integration of the magnetic field scalar product by carrying out the volume integral in the region outside the scatterer (i.e. for ),

 ∫r>1drh(e)n,m,α(r)⋅h(e)n,m,β(r)=−z3/2αz3/2βN(e)n,αN(e)n,βlimη→0∫∞1r2hn(zαr)hn(zβr)e−ηr2dr (9d) =z3/2αz3/2βN(e)n,αN(e)n,βzβhn(zα)h′n(zβ)−zαh′n(zα)hn(zβ)z2α−z2β≡Θ(e)h,n(zα,zβ), (9e)

where we used Eq.(103a) in ref. [12] to evaluate the integral in the first line of this equation (reproduced in Eq.(31h)). This integral would be ill-defined had we not invoked the ‘Killing’ regularization, described in ref. [12] and appendix A, which tames the integrand’s exponentially divergent amplitude with a factor in order to yield the finite analytic expression of Eq.(9e) in the limit. For later convenience, we defined the right hand side of Eq.(9e) as a magnetic scalar product ‘surface term’, , evaluated at the sphere radius.

Simply adding Eqs.(9c) and (9e) immediately leads to the orthogonality of the electric type RSs when integrating their magnetic fields over all space,

 ∫∞drμ(r)hα(r)⋅hβ(r)=0, (10a) where we dropped the (e) superscript in Eq.(10a) since an analogous orthogonality holds for (h) type RSs (orthogonality holds for the RSs of any dispersion-free media).

Numerical plots of the integrands in Eq.(9a) and (9d) can help us appreciate the complexity of what is accomplished in Eq.(10a). The real and imaginary parts of the magnetic field integrands, including the multiplicative constants, of Eq.(9a) () and Eq.(9d) () are plotted in fig.1(a) for up to for the first two electric dipole RSs at and of a high index dielectric sphere (cf. Table 1). The integrands are continuous at the sphere surface as one would expect for magnetic field energy. The graphs are then extended in Fig. 1(b) to with a Gaussian killing factor of . A high accuracy numerical integral over the entire axis with a sufficiently small value of does indeed tend towards zero in both the real and imaginary parts as is plausibly indicated by visual inspection of Fig. 1(b) and 1(d) .

An analogous result holds for orthogonality of the scalar product in terms of the electric field even though the situation is a bit more complicated due to the discontinuous contributions of the electric field normal to the scatterer surface. Nevertheless, using the analytical results of section C, one again finds,

 ∫∞drε(r)eα(r)⋅eβ(r)=0, (10b)

a result which is illustrated graphically by plotting the electric fields integrands in Eq.(45a) and (45d) for up to and fig.1(d) for up to . The mathematics for the electric field integration was somewhat more difficult for the electric case since the required integrals could not be found in standard references like that of Watson [27], but they are provided in reference [12].

The challenge of a purely numerical evaluation of field orthogonality is now clear: one can look to enhance numerical accuracy by choosing a small factor, but this will simultaneously increase the amplitudes outside the particle, thus rendering numerical accuracy all the more difficult. Given modern computation abilities, this is of course possible, at least in simple cases illustrated here, but clearly numerically intensive as discussed in appendix A. Note however, that treatments of ‘orthogonality’ in dispersive media and resonant state normalization will generally require the use of the full electromagnetic field vector, , as discussed in section 2.3.

An important consequence of the above results is found by combining the electric and magnetic integrals of Eq.(C) and Eq.(2.2) to find the electromagnetic scalar product of Eq.(8), which now reads,

 ⟨Ψα|Γ∣∣Ψβ⟩=∫Vrdr[ε(r)e(e)α(r)⋅e(e)β(r)−μh(e)n,m,α(r)⋅h(e)n,m,β(r)]+⎧⎪ ⎪ ⎪⎨⎪ ⎪ ⎪⎩r2k3/2αk3/2βN(e)n,αN(e)n,βhn(rkα)ξ′n(rkβ)rkβ−ξ′n(rkα)rkαhn(rkβ)kα−kβ⎫⎪ ⎪ ⎪⎬⎪ ⎪ ⎪⎭∂Vr=0, (11a) where the volume Vr was derived as the sphere surface, but it can be directly generalized to be any volume r/R>1, containing the scatterer, with ∂Vr designating the surface of Vr. The term in curly brackets can be formulated directly in terms of the surface RS fields, {…}∂Vr=−i∣∣r2(h(e)n,m,β(rkβˆr)×e(e)n,m,α(rkαˆr)+e(e)n,m,β(rkβˆr)×h(e)n,m,α(rkαˆr))⋅ˆr∣∣r=Rkα−kβ, (11b)

which allows the volume/surface integrals to be generalized to non-spherical volumes provided they enclose the scattering system. Analogous results hold for the magnetic mode scalar products. The result of Eq.(11) is applicable to RSs for arbitrary scattering geometries.

### 2.3 Resonant state normalizations

When RSs are introduced in terms of the equations of motion as done in Eq.(3), there is an overall normalization degree of freedom as seen in Eq.(2.1). The physically correct value of normalization will be fixed by linear response functions like Green functions or the Mie theory coefficients of section 3, but we can foresee the natural normalization condition, which for dispersionless media reads,

 (12)

where the structure of the field equations of Eq.(2), requires that the medium metric, , be included in the definition of scalar product as already remarked in Eq.(8). A number of authors have proposed to normalize to unity in analogy with quantum mechanics, but the normalization proposed here adopts a field theory standpoint where states are normalized with respect to their energy. For example, the electromagnetic ‘Lagrangian’ energy of the RS is given by .

For temporally dispersive media, the RS normalization condition must be generalized to read,

 ⟨Ψα|ddω(ωΓ)|Ψα⟩=∫∞dr{d[ωε(r,ω)]dω∣∣∣ωαe2α(r)−d[ωμ(r,ω)]dω∣∣∣ωαh2α(r)}=zα, (13)

where the factors and are required in order to account for the energy associated with material degrees of freedom [28] such that the Lagrangian energy, , of both electromagnetic and material degrees of freedom is (provided certain physical assumptions concerning the material response functions). One should remark that Eq.(13) is the only normalization condition required, since Eq.(12) is simply the reexpression of Eq.(13) under the assumption of frequency independent (i.e. lossless) constitutive parameters.

From the expressions for the RSs in Eqs.(2.1) it is relatively straightforward, but tedious, to obtain the finite values of the integrals in Eq.(13) thanks once again to the ‘Killing’ regularizations of [12]. Details of the calculation are relegated to appendix D where one finds:

 (14a) [N(h)n,α]2=(εα−1)ξ2n(zα)+(μα−1){[ξ′n(zα)]2−n(n+1)h2n(zα)μα}+ωα2{Ξ(h,−)nddωlnε(ω)∣∣∣ωα+Ξ(h,+)nddωlnμ(ω)∣∣∣ωα}, (14b) where we defined: Ξ(e,±)n ≡εα[ξ′n(zα)]2+μαξ2n(zα)−n(n+1)h2n(zα)εα±hn(zα)ξ′n(zα) (14c) Ξ(h,±)n ≡μα[ξ′n(zα)]2+εαξ2n(zα)−n(n+1)h2n(zα)μα±hn(zα)ξ′n(zα). (14d)

By inspection of Eqs.(14), we see that the normalization factors depend on both RS frequency and particle size. For dispersionless media, Eqs.(14) simplify to expressions which only depend on the size parameter, and the normalized, real-valued, constitutive parameters of the sphere, and  [29]:

 [N(e)n,α]2no disp.=ξ2n(zα)(μ−1)+{[ξ′n(zα)]2+n(n+1)εh2n(zα)}(ε−1) (15a) [N(h)n,α]2no disp.=ξ2n(zα)(ε−1)+{[ξ′n(zα)]2+n(n+1)μh2n(zα)}(μ−1). (15b)

Therefore, in the dispersionless case, the problem becomes scale invariant, and the RS expansion is a complete solution for all size scales in addition being a solution for all frequencies.

As in the case for RS orthogonalization integrals, we can gain physical insight into integrals of Eq.(12) by plotting the radial dependence of the normalized Lagrangian integrands (after the analytic angular integration). We simply set in these graphs because the exponential divergence only occurs for values far greater than the extent of the plots for the RSs , , and . For such low loss RSs, a numerical ‘verification’ of Eq.(12) or Eq.(13) is not difficult provided that one judiciously truncates integrations before the onset of exponential divergence. This is not at all the case however for the RS , which has a much larger imaginary part. As partially visible in the plot of Fig. 2(a), the exponential divergence sets in immediately outside the sphere so it is impossible to verify the normalization of Eq.(12) without regularizing the integral. Similar behavior is observed for the RS of lossy spheres studied in section 4.2.

## 3 Resonant state expansions of Mie theory

### 3.1 Spectral expansions

In the interest of completeness, essential elements of Mie and linear response theory are recalled in appendix B (cf. [16, 29, 30, 31, 32, 33] for more details and applications). Since Mie theory predates the matrix response formulations, it remains common in the literature to keep the traditional notations for the diagonal ‘-matrix’ elements as and , and for internal field coefficients, and .

The internal field coefficients, , can be expressed as meromorphic functions of the complex frequency, , with all the poles located in the lower half-plane of complex frequencies,[29, 33]

 Ω(e,h)n(ω,R)=e−2iz∞∑ℓ=−∞r(Ω)n,αz−zα, (16)

where the RS indices, are functions of the multipole indices and the discrete RS index (see the discussion following Eq.(2.1) and Fig. 3). For dimensional reasons, the frequency appears in the pole expansion on the right hand side of Eq.(16) via the size parameter . One should take care however that this development is uniquely a function of , only in the case of dispersionless media. In the presence of dispersion, the values, have a non trivial dependence on the sphere radius, , which generally can only be determined numerically (by solving Eq.(6) for example). The functions have a similar meromorphic pole expansions in terms of frequency , with the exception of an additional non-resonant term, [16, 29, 33]

 T(e,h)n(ω,R)=Ane−2iz−12+e−2iz∞∑ℓ=−∞rn,αz−zα, (17a) where An are sign factors which for electric and magnetic modes respectively are, A(e)n≡(−)n+1A(h)n≡(−)n, (17b)

where again is a function uniquely of for dispersionless media. We also remark that the Mie coefficients are not entirely determined by the resonant eigenvalues, , but that they also require the presence of a non-resonant term and multiplicative factors invoking an phase factor. [29, 33, 16]

### 3.2 Spectral residues

Given the expressions for , and in Eqs.(17) and (16) in terms of the and functions (defined in section 33 of the appendix), the residues in Eqs.(16) and (17) can be calculated from Mie theory either analytically or numerically as,

 r(Ω)n,α=e2izαcRddωDn(ω,R)∣∣ω=ωα,rn,α=−Nn(ωα,R)e2izαcRddωDn(ω,R)∣∣ω=ωα, (18a) which can be rewritten as, r(Ω)n,α=e2izαlimω→ωαz−zαDn(ω,R),rn,α=−Nn(ωα,R)r(Ω)n,α. (18b)

We remark from the second equality in Eq.(18b) that the relation between the -matrix residues, , and the internal field residues, , is purely analytic. This relation can be rewritten in a more transparent form by invoking the respective RS conditions of Eqs.(6) such that the expressions of and can be rewritten:

 −N(e)n(ωα,R)=zαψ′n(ραzα)jn(zα)−εαjn(ραzα)ψ′n(zα)iρα→ψ′n(ραzα)ραξ′n(zα)=1γe,n,α, (19a) for the electric modes and, −N(h)n(ωα,R)=zαψ′n(ραzα)jn(zα)−μαjn(ραzα)ψ′n(zα)iμα→ψ′n(ραzα)μαξ′n(zα)=1γh,n,α, (19b)

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

 ξ′n(z)ψn(z)−ξn(z)ψ′n(z)=i, (20)

and recalled the definitions of the functions defined in Eq.(6). Using the results of Eq.(19) in Eq.(18), we find that for both electric and magnetic modes,

 rn,α=r(Ω)n,αγn,α, (21a) which is a result that could have been anticipated in view of linear response theory and the role of the γn,α functions in expressing the RS wave functions in Eq.(2.1).

With the newly derived result of Eq.(21a) in hand, one can find the formula for by carrying out a numerical or analytical evaluation of Eq.(18b). This won’t be necessary however, because the connection between the analytic results of Mie theory are obtained by remarking that since all the resonant states were constructed to have a amplitude dependence in the far field, the -matrix residues should be equal to the inverse square of the complex normalization factor up to a phase factor which we find as:

 rn,α=e2izαiN2n,α≡e2izαRn,α, (21b)

where we defined the Mie coefficient residue, . In the case of Mie theory, the validity of Eq.(21b) can be verified either analytically or numerically, but one should remark that Eq.(21b) is a general relation for -matrix residues for objects of arbitrary form whose factors have been chosen such that the general normalization condition of Eq.(13) is satisfied. It is interesting to remark that the simplicity of the relation between the Mie residue, , and RS normalization, in Eq.(21b), is a consequence of the judicious RS normalization adopted in Eq.(13). Had we imposed a normalization other than that of Eq.(13), we would have to corrected for it at this stage of the calculations.

### 3.3 Implementation: practitioners guide

The steps to numerically implement the results of this paper in the case of spherical scatterers are as follows:

1. For each required multipole order , and mode type , the only non-trivial step is to find a sufficient number of the resonant state size parameters, , by numerically solving the transcendental Eq.(6), and store its solutions along with the parameters associated with these solutions. What a sufficient number of RSs is, depends on the particle radius , the dispersion relations of the constitutive parameters, and the degree of numerical accuracy desired, but can be as low as 1 for dispersive materials. When only a few RSs and multipole orders are required, graphical solutions of with numerical refinements may suffice (see section 4 for examples, but one should exercise caution since important RSs can occur far from the real axis).

2. Calculate the square of the complex-valued RS norms, , using Eq.(14).

3. Insert the values of found in step 2 into Eq.(21b) to determine the residues, , for the meromorphic expression of Eq.(17) for the scattered field Mie coefficients, and . Cross sections contributions, , can then be determined by the standard formulas recalled in Eq.(37) of the appendix B.2.

4. If one needs to determine the fields inside the sphere, use the values of found in step 3 and the coefficients found in step 1 in Eq.(21a) to calculate the residues, , of the meromorphic expression of Eq.(16), for the internal field Mie coefficients and (defined in Eq.(33b)).

An example of the results of this procedure is given in Fig. 3which graphically illustrates the respective electric and magnetic modes with multipoles, and . The first few RSs for each mode are given as blue dots in the lower half plane with on the vertical axis, while the associated contributions to the cross section efficiencies, , are plotted on the positive vertical scale for a spherical scatterer with relative permittivity of , with respect to the background medium. One remarks that the determination of the RS complex eigenvalues in step 1 is the only information that has to be obtained numerically. The normalized ‘strengths’ of the modes (and thereby the cross sections) are determined by the analytic formulas of step 2 and 3. Since the problem is scale invariant, the RSs and cross sections only depend on the size parameter . Similar calculations can be performed for dispersive media, but in that case, the problem is no longer scale invariant and the position of the RS will depend on the size of scatterer as discussed in section 4.

## 4 Mode Normalization: Numerical implementation

We now provide a few concrete numerical examples of RS’s values, and their normalization factors for different materials of interest in nanophotonics and plasmonics. One begins by carrying out step 1 of section 3.3 for determining the eigenvalues by respectively rewriting the magnetic and electric RSs eigenvalue conditions of Eqs.(6). One interesting way to express the mode conditions is in terms of reduced logarithmic derivatives of the Ricatti-Bessel functions:

 [xln′ψn(x)]x=ραzα−εα[xln′ξn(x)]x=zα =0 (22a) [xln′ψn(x)]x=ραzα−μα[xln′ξn(x)]x=zα =0, (22b)

where we recall from Eq.(5c) that , , and are in general frequency dependent constitutive parameters evaluated at the RS frequency such that the left hand sides of Eqs.(22) have a non-trivial dependence on . A variety of numerical techniques exist for solving Eqs.(22), but it can be instructive to start with graphical solutions as will be illustrated in Figs (4) and (5) below.

An advantage of comparing with Mie theory is that a formal or numerical evaluation of the Mie theory predictions using Eqs.(18) allows one to confirm the normalization values predicted by Eq.(3.2) and Eq.(14) from the evaluation of the RS normalization integrals. In all approaches, results agree up to reasonable computational accuracy.

### 4.1 Resonant states of high index dielectric spheres

In recent years, high index dielectrics are being considered in nano-optics as an alternative to plasmonics [34]. In contrast to the highly dispersive and lossy metals used in plasmonics, the high index materials may be treated as non-dispersive and lossless over significant wavelength ranges. Making the approximation of a non-magnetic material and real non-dispersive permittivity, we can use the simplified expressions of Eq.(15) for the RS normalization coefficients.

Numerical results are shown in Table 1 for a few low-order modes with the same parameters as studied in Fig. 3.More detailed analysis of such modes may be found in [33]. There are an infinite number of modes for each value of , and they follow an asymptotic pattern resembling that of zeros of Bessel functions. Beginning with electric modes, Table 1 gives the values of and the normalization factors for the first two dipole and quadrupole modes. An interesting feature of the RSs in this case is that for even (odd), there is one electric (magnetic) mode with purely imaginary , labeled with . Table 1 shows that for imaginary RS eigenvalues, the interaction ‘strengths’ (residues), , have purely imaginary values.

### 4.2 Resonant states of ‘Drude’ model conductors

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

Data on the complex refractive index of gold as a function of frequency are available in Volume 1 of he comprehensive collection due to Palik [35]. 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 as [18],

 ε(ω)=ε∞−ω2pω2+iωγD, (23)

where , , and . The comparison given in Fig. 9 shows that the Drude model works well only for wavelengths longer than .

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

 ε(ω)=ε∞−ω2pDω2+iωγD−s1ω2p1,Lω2−ω2p1,L+iγ1,Lω−s2ω2p2,Lω2−ω2p2,L+iγ2,Lω. (24)

Here the parameters are given in eV: , , , , , , , and . The conversion of to electron volts is achieved by dividing the frequency results by . As can be seen from Fig. 9, this model works well down to around . An interesting feature of the dispersion model of Eq.(24) is that it supports bulk plasmons. [37] These occur for complex wavelengths or frequencies at which . For the dispersion model of Eq.(24) they are at: , , . 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, , of zero for these longitudinal ‘resonant states’.

Silver adheres more closely to the Drude model than does gold in the visible and near infrared. Yang et al [38] suggest the following fit:

 εAg(ω)=ε∞−ω2pω2+iωγD (25)

where , , and . The comparison given in Fig. 8 shows that the Drude model works well for wavelengths longer than around .

Fig. 4 shows the modulus of the dispersion equation of Eq.(22a) for and for a gold sphere of radius , both for the Drude model of Eq.(23) and for the more elaborate model of Eq.(24). The former has its minimum at , in good agreement with the value reported in Sauvan et al [18]. The latter has a slightly different value: . From Table 2, the values of at the resonance and consequently the normalization factor disagree more significantly for the more elaborate model.

Fig. 5 shows the inverse modulus of the dispersion relation of Eq.(22) for a gold sphere of radius . The plot is more complicated for the more elaborate dispersion model of Eq.(24) 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 . With the model of equation (24), this becomes . Note the significant decrease in the imaginary part- as Fig. 9 shows, around this wavelength, the Drude model is only a crude approximation to the experimental data.

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

## 5 Conclusions

This work shows that the analytic regularization of the scalar products of resonant state wave functions provides a polyvalent tool for both theoretical and numerical analysis of open systems. Although the derivations of the expressions employed here required rather extensive investigations, employing the methods of generalized function theory, their implementation are readily accessible. Among our demonstrations, we have derived closed-form analytic expressions for the normalization of the light scattering resonant states of spherical scatterers with arbitrary temporal dispersion including magnetic degrees of freedom. Spectral expansions, rendered possible by this method, open up numerous perspectives for wide-band and time-domain investigations of open systems. Illustrative figures were presented to help the reader appreciate the interest of these formulas to replace/enhance purely numerical treatments. The tabular results we have presented should prove useful as a reference point for other researchers. These results also reveal insights into investigations of radiation-matter interactions and radiative decay in general.

## 6 Acknowledgements

Research conducted within the context of the International Associated Laboratory for Photonics between France and Australia (LIA ALPhA). This work has been carried out thanks to the support of the A*MIDEX project (no. ANR-11-IDEX-0001-02) funded by the Investissements d’Avenir French Government program, managed by the French National Research Agency (ANR).

## Appendix A Killing Mie softly

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

We now give an example of the results in the references [12, 24]. The integral chosen is that over the product of the spherical Bessel functions and , where and are (possibly complex) wavenumbers: see equation (LABEL:snr18).

 Ijy(n,K,k,η) = ∫∞0x2exp(−ηx2)jn(Kx)yn(kx)dx = π2{exp[−(K2+k2)/(4η)]2πη[−H(n+1/2,k,K,η)+(n+1/2)h−1,n+1/2(−Kk2η)]}.

The analytic result for this integral from [24] is given by the rightmost expression in (LABEL:snr18). Here denotes an associated Bessel function [39], and is the following finite-range integral:

 H(b,k,K,η)=∫K/k1u(b−1)exp[Kk4η(u+1/u)]du. (27)

Using expansions given by Luke [39] and evaluating the finite integral in (27) by direct numerical integration, the expression (LABEL:snr18) may be verified for arbitrary choices of the parameters , and .

The asymptotic treatment given in [24] which takes the limit as is quite lengthy and complicated. However, it yields a simple result, which is in keeping with the well-known result from Watson [27]:

 ∫zzCμ(kz)Dμ(lz)dz=z{kCμ+1(kz)Dμ(