Resonant-state expansion of dispersive open optical systems

A resonant-state expansion (RSE) for open optical systems with a general frequency dispersion of the relative permittivity, described by a finite number of simple poles, is presented. As in the non-dispersive case, the RSE of dispersive systems converts Maxwell's wave equation into a linear matrix eigenvalue problem in the basis of unperturbed resonant states, in this way numerically exactly determining all relevant eigenmodes of the optical system. This dispersive RSE is verified by application to the analytically solvable system of a sphere in vacuum, with a dispersion of the dielectric constant described by the Drude and Drude-Lorentz models. We calculate the change of the optical modes when converting the sphere material from gold to non-dispersive silica and back to gold, and evaluate the accuracy using the exact solutions.

Any optical system is characterized by its resonances which are a cornerstone of physics. The concept of resonant states (RSs) is a mathematically rigorous way of treating the resonances. Formally, RSs are the optical eigenmodes of the system, i.e. the eigen-solutions of Maxwell's wave equation, which satisfy the outgoing wave boundary conditions. In open optical systems the RS eigenfrequencies ω n are generally complex, which physically reflects the fact that the energy leaks out of the system. The real part Re(ω n ) gives the position of the resonance, while the imaginary part Im(ω n ) gives its half width at half maximum, also determining the quality factor of the resonance as Q n = |Re(ω n )/[2 Im(ω n )]|.
We have recently developed the resonant-state expansion (RSE), a rigorous perturbative method for calculation of RSs, which is treating perturbations of open optical systems of arbitrary strength and shape [1]. We have shown its advantages over established computational methods in electrodynamics, such as finite difference in time domain (FDTD) and finite element method (FEM), in terms of accuracy and efficiency [2]. Specifically we note that the RSE (i) uses the natural discretization in the frequency domain provided by RSs, (ii) reduces the solution of Maxwell's wave equation to a linear matrix eigenvalue problem, and (iii) produces all RSs originating from the basis states in a single calculation, avoiding spurious solutions. This enables the RSE to determine numerically exactly all the RSs in a frequency range of interest, with the accuracy limited by the basis truncation only.
Other methods use artificial discretization in space and/or in time/frequency domain and the approximation imposed by perfectly matched layers (PMLs) at the system boundaries, both giving rise to issues. FEM, for example, determines RSs one by one, iteratively solving a nonlinear equation with unknown analytics -it is therefore impractical if not impossible to verify that all RSs within a complex frequency area have been found. In FDTD, RSs can be found by fitting the calculated time evolution by a sum of RSs. Only RSs which have been excited in the simulation are visible, and the fitting pro-cedure does not uniquely determine the number of RSs. Additionally, the spatial discretization and PMLs give rise to spurious solutions.
The reason why the RSE was not available until recently is in the fact that RSs with complex eigenfrequencies have wave functions which are exponentially growing in space away from system, and the proper general normalization of such RSs was not known. The issues with the normalization has been discussed recently, e.g. in Refs. [3] and [4] where different numerical procedures were suggested. At the same time, the correct normalization corresponding to the spectral representation of the Green's function (GF) of Maxwell's wave equation in terms of RSs is at the heart of the RSE. The correct analytic normalization is contained in our first work on the RSE [1]. This normalization was recently generalized to an arbitrary surface of integration and to optical systems with dispersion which allowed for an exact theory of the Purcell effect [5], almost 70 years after its discovery.
So far the RSE has been applied to non-dispersive systems of different dimensionality and geometry [1,2,[6][7][8]. However, almost all realistic systems, even dielectrics such as glass, have a frequency dispersion of the refractive index. We have recently found [9] that the direct substitution of an Ohm's law dispersion into the non-dispersive RSE maintains its linearity. The Ohm's law dispersion can be a reasonable approximation for some materials whose relative permittivity (RP) is mainly determined by their dc conductivity or when the dispersion can be approximated by a term linear in the light wavelength over the frequency region of interest. However, metals are better described by the Drude model [10], and a significant improvement is achieved by adding Lorentzian terms [11], which is further refined by using complex weights (residues) of the frequency poles called critical points (CPs) of the RP [12][13][14].
In this Letter we generalize the RSE to the case of a frequency dispersion of the RP with a countable number of poles, suited to describe the RP of any physical material. We verify it on the exactly solvable system of spherical metal and dielectric nano-particles. We start with a dispersive basis of RSs with the wave functions E n (r) (r is the position vector) and frequencies ω n being the eigen-solutions of Maxwell's wave equation and the electric fields E n (r) satisfying the outgoing wave boundary conditions [7]. The dispersive RP tensor ε(r, ω) of an unperturbed open optical system is taken in the form of a function in the complex frequency plane expressed asε whereε ∞ (r) is the high-frequency value of the RP and Ω j are the resonance frequencies (poles) of the RP determining the dispersion, with the weight tensorsσ j (r) corresponding to generalized conductivities of the medium at these resonances. The Lorentz reciprocity theorem requires that all tensors in Eq. (2) are symmetric, and the causality principle requires thatε * (r, ω) =ε(r, −ω * ) [15]. Therefore, for a physically relevant dispersion, each pole of the RP with a positive real part of Ω j has a partner at Ω −j = −Ω * j withσ * −j =σ j , while poles with zero real part of Ω j have realσ j . The Ohm's law dispersion of the RP corresponds to the sum in Eq. (2) replaced by a single term with Ω 0 = 0 andσ 0 (r) being the dc conductivity tensor. The Drude model of metals consists of two poles with Ω 0 = 0, Ω 1 = −iγ, and σ 1 (r) = −σ 0 (r). The Drude-Lorentz model introduces additional poles at ω = Ω j with j = ±2, ±3, . . . and complex conductivitiesσ j . Fig. 1 provides an example of the Drude and Drude-Lorentz models approximating the measured complex refractive index n r (ω) = ǫ(ω) of gold [10], with the parameters taken from Refs. [3] and [13], respectively, and used in the following for illustration of the dispersive RSE. In particular, we used for the Drude model Ω 1 = −92.8 i meV, σ 1 = −744 eV, and ε ∞ = 1 [3]. For the Drude-Lorentz model we used .20 e iπ/4 eV, and ε ∞ = 1.54 [13].
The GF of Maxwell's wave equation has an infinite countable number of simple poles in the complex frequency plane, and therefore has the spectral represen-tationĜ where the sum is taken over all RSs, and ⊗ denotes the dyadic product of vectors. The spectral representation Eq. (3) requires that the RSs are normalized according (green error bars) and approximated by the Drude model [3] (thick blue lines), and by the Drude-Lorentz model with two pairs of CPs [13] (thin red lines).
to [5] 1 where F n = (r · ∇)E n , V is an arbitrary simply connected volume with a boundary surface S V enclosing the inhomogeneity of the system, and the derivative ∂/∂s is taken along the outer surface normal. Note that for static modes (ω n = 0), the volume of integration can be extended to the entire space and the surface term vanishes, since the electric field in such modes decay or just vanishes outside the system. All other modes instead have an electric field exponentially growing outside the system, so that the normalization has to be evaluated for a finite V . Substituting Eq. (3) into Maxwell's wave equation for the GF and using Eq. (1) we obtain which has to be satisfied for any ω. For the dispersion of the RP in the form of Eq. (2), we find with the help of the algebraic identity that Eq. (5) splits into the following closure relation and sum rules Now, again using the algebraic identity where and combining Eq. (3) with sum rules Eq. (8) for Ω 0 = 0 and Ω j = 0 according to the terms in Eq. (9), we find an additional spectral representationĜ j ω of the GF for each pole in the RP: The Ohm's law dispersion introduces a ω = 0 pole in the RP which leads to the sum rule Eq. (8) corresponding to Ω 0 = 0. This sum rule results in the representation G 0 ω (r, r ′ ) of the GF given by Eq. (11) with W 0 n (ω) = 1. Note however, that the ω = 0 pole is implicitly present already in the non-dispersive system owing to the longitudinal ω n = 0 modes [2]. As a result, the sum rule Eq. (8) with Ω 0 = 0 holds also without dispersion [7,8], due to the constant term in the RP, such asε ∞ (r) in Eq. (2). This explains why Ohm's law does not need any significant reformulation of the RSE compared to the non-dispersive case and does not require an extension of the basis of RSs [9].
We now solve a perturbed Maxwell's wave equation equivalent to Eq. (1), with the unperturbed RPε(r, ω) replaced by a perturbed one,ε(r, ω) + ∆ε(r, ω), with the perturbation ∆ε(r, ω) of the form of Eq. (2) described by the perturbations ∆ε ∞ (r) ofε ∞ (r) and ∆σ j (r) ofσ j (r) inside the unperturbed system. We find the electric field E(r) and the eigenfrequency ω of a perturbed RS by using the unperturbed GF in different representations Eq. (11) for the corresponding terms of the RP, yielding This integral equation is then converted to a matrix equation by expanding the perturbed RS into the basis of unperturbed ones, by using expansions Eq. (11) of the GF, and by equating the coefficients at different basis functions E n (r). The result is the eigenvalue equation This is the linear dispersive RSE. The perturbation matrix V nm (ω) represents the change ∆ε(r, ω) of the RP for any physical dispersion described by Eq. (2).
The linear dispersive RSE is suited for both dispersive and non-dispersive unperturbed systems with perturbations which do or do not increase the number of poles in the RP. When increasing the number of non-zero poles of the RP, the number of poles of the GF is increased too [16], and the size of the RSE basis is extended by an additional countable infinite number of RSs for each nonzero pole of the RP, with the RS eigenfrequencies asymptotically approaching this pole. Poles of the RP with finite weight in the perturbed system but zero weight in the unperturbed system are included in the basis by taking the limit of the pole weight tending to zero. In this limit, the pole-related RSs have frequencies converging to the pole but refractive indices taking separate discrete values, as detailed below.
To illustrate the method and evaluate its convergence, we show in Figs. 2 and 3 the transverse magnetic (TM) eigenmodes with l = 1 (l is the orbital number) of spheres made of a dispersive material (gold) and a non-dispersive material (sand, n r = 1.5) in vacuum, and perturbations which transform gold to sand in Fig. 2 and sand to gold in Fig. 3. The eigenmodes of the sand and gold spheres in vacuum were taken in the analytic form [2] and normalized according to Eq. (4), see Ref. [5] for explicit analytic expressions. The radius of the sphere R = 200 nm is chosen such that both Drude and Drude-Lorentz approximations of the gold dispersion shown in Fig. 1 are valid for the frequency of the fundamental surface plasmon (SP) mode shown in Figs. 2 and 3 by arrows.
We select a finite number N of RSs for the RSE basis, including only RSs satisfying the condition |n r (ω n )ω n | < ω max . This excludes RSs having a wavevector in the medium above ω max /c, which is the case of large ω n or large |n r (ω n )| with ω n close to the poles of the dispersion. This basis selection can be optimized in the future. The RSE results for the perturbed eigenmodes are compared with the analytic solutions, and the relative errors are shown in Figs. 2(b) and 3(b) for different ω max , demonstrating a high accuracy given the strong perturbation. For the present geometry, N is approximately proportional to ω max , with N = 456 for ω max = 200 eV. The observed 1/N 3 convergence to the exact solution is comparable to that of the non-dispersive RSE [1,2].
Going from gold to sand [ Fig. 2(a)] the RSE reproduces the RSs of the non-dispersive sand sphere, and additionally produces a number of quasi-degenerate RSs at the Drude and Lorentz poles. These RSs are present in the system since in the linear RSE the same poles of the dispersion are present before and after perturbation. Poles which have zero weight in a system still lead to an infinite series of RSs, with frequencies at the pole position, but corresponding to different refractive indices, as exemplified in the inset of Fig. 3(a). For the sphere geometry, they can be calculated analytically by taking the limit of the pole weight to zero in the secular equation. A perturbation which creates a finite weight of the pole lifts the degeneracy of these RSs as exemplified in Fig. 3. We now compare this result with an alternative dispersive RSE approach which uses a non-dispersive system as basis and creates the additional RSs due to the poles of the dispersion via the nonlinearity of the resulting generalized eigenvalue problem. Assuming that the unperturbedε(r) has no dispersion, the only valid sum rule in Eq. (8) is the one with Ω j = 0 which provides only one alternative GF representationĜ 0 ω (r, r ′ ) with W 0 n (ω) = 1 in Eq. (11). ReplacingĜ j ω byĜ 0 ω in Eq. (12) results in a nonlinear dispersive RSE which is a direct generalization of the original, nondispersive RSE [1,2].  [17], extending the basis of unperturbed RSs by a factor of M . We illustrate this alternative method for the Drude dispersion of the perturbed system, for which Eq. (16) is a quadratic matrix problem. For the same basis cut-off ω max as used for the linear dispersive RSE, the energies of the Fabry-Pérot RSs are reproduced with a similar accuracy, see Fig. 3(b). However, the SP mode has about 2 orders of magnitude larger error and modes around the Drude pole are also having orders of magnitude larger error as shown in the inset of Fig. 3(b). This can be understood considering that in the nonlinear RSE the basis does not contain the pole RSs, and is therefore less suited to describe the RSs close to the RP poles.
In conclusion, the presented generalization of the RSE to materials with an arbitrary frequency dispersion of the relative permittivity, described by a finite number of simple poles, is extending the applicability of the RSE to all relevant open optical systems, paving the way to its widespread use in electromagnetic simulation. This work was supported by the Cardiff University EP-SRC Impact Acceleration Account EP/K503988/1, and the Sêr Cymru National Research Network in Advanced Engineering and Materials. The authors acknowledge discussions with M. D. Doost and H. Sehmi.