Resonant-state expansion of dispersive open optical systems

# Resonant-state expansion of dispersive open optical systems

## Abstract

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 are generally complex, which physically reflects the fact that the energy leaks out of the system. The real part gives the position of the resonance, while the imaginary part gives its half width at half maximum, also determining the quality factor of the resonance as .

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 procedure 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); (6); (7); (8); (2). 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 ( is the position vector) and frequencies being the eigen-solutions of Maxwell’s wave equation

 Extra open brace or missing close brace (1)

and the electric fields satisfying the outgoing wave boundary conditions (7). The dispersive RP tensor of an unperturbed open optical system is taken in the form of a function in the complex frequency plane expressed as

 ^ε(r,ω)=^ε∞(r)+∑ji^σj(r)ω−Ωj, (2)

where is the high-frequency value of the RP and are the resonance frequencies (poles) of the RP determining the dispersion, with the weight tensors 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  (15). Therefore, for a physically relevant dispersion, each pole of the RP with a positive real part of has a partner at with , while poles with zero real part of have real . The Ohm’s law dispersion of the RP corresponds to the sum in Eq. (2) replaced by a single term with and being the dc conductivity tensor. The Drude model of metals consists of two poles with , , and . The Drude-Lorentz model introduces additional poles at with and complex conductivities . Fig. 1 provides an example of the Drude and Drude-Lorentz models approximating the measured complex refractive index 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  meV,  eV, and  (3). For the Drude-Lorentz model we used  meV,  eV,  eV,  eV,  eV,  eV, and  (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 representation

 ^Gω(r,r′)=c2∑nEn(r)⊗En(r′)2ωn(ω−ωn), (3)

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 to  (5)

 1+δ0,ωn = ∫VEn(r)⋅∂(ω2^ε(r,ω))∂(ω2)∣∣ ∣∣ωnEn(r)dr (4) +c22ω2n∮SV(En⋅∂Fn∂s−Fn⋅∂En∂s)dS,

where , is an arbitrary simply connected volume with a boundary surface enclosing the inhomogeneity of the system, and the derivative is taken along the outer surface normal. Note that for static modes (), 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 . Substituting Eq. (3) into Maxwell’s wave equation for the GF and using Eq. (1) we obtain

 ∑nω2^ε(r;ω)−ω2n^ε(r;ωn)2ωn(ω−ωn)En(r)⊗En(r′)=^1δ(r−r′) (5)

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

 1ωn(ω−ωn)(ω2ω−Ωj−ω2nωn−Ωj) =1ωn−Ωj+ωω−Ωj(1ωn−1ωn−Ωj) (6)

that Eq. (5) splits into the following closure relation

 12∑n^ε(r,ωn)En(r)⊗En(r′)=^1δ(r−r′) (7)

and sum rules

 ∑nEn(r)⊗En(r′)ωn−Ωj=0. (8)

Now, again using the algebraic identity

 1ωn(ω−ωn)−1ωωn+Ωjω2(ωn−Ωj)=Wjn(ω)ω(ω−ωn), (9)

where

 Wjn(ω)=ωnωω−Ωjωn−Ωj, (10)

and combining Eq. (3) with sum rules Eq. (8) for and according to the terms in Eq. (9), we find an additional spectral representation of the GF for each pole in the RP:

 ^Gjω(r,r′)=c2∑nWjn(ω)En(r)⊗En(r′)2ω(ω−ωn). (11)

The Ohm’s law dispersion introduces a pole in the RP which leads to the sum rule Eq. (8) corresponding to . This sum rule results in the representation of the GF given by Eq. (11) with . Note however, that the pole is implicitly present already in the non-dispersive system owing to the longitudinal modes (2). As a result, the sum rule Eq. (8) with holds also without dispersion (7); (8), due to the constant term in the RP, such as 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 replaced by a perturbed one, , with the perturbation of the form of Eq. (2) described by the perturbations of and of inside the unperturbed system. We find the electric field 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

 E(r) = −ω2c2∫^Gω(r,r′)Δ^ε(r′,ω)E(r′)dr′ (12) = −ω2c2∫[^G0ω(r,r′)Δ^ε∞(r′,ω) +∑j^Gjω(r,r′)iΔ^σj(r′,ω)ω−Ωj]E(r′)dr′.

This integral equation is then converted to a matrix equation by expanding the perturbed RS into the basis of unperturbed ones,

 E(r)=∑ncnEn(r), (13)

by using expansions Eq. (11) of the GF, and by equating the coefficients at different basis functions . The result is the eigenvalue equation

 (ωn−ω)∑m[2δnm+Vnm(∞)]cm=ωn∑mVnm(ωn)cm (14)

which is linear in , with the perturbation matrix

 Vnm(ω)=∫En(r)⋅Δ^ε(r,ω)Em(r)dr. (15)

This is the linear dispersive RSE. The perturbation matrix represents the change of the RP for any physical dispersion described by Eq. (2). In the absence of dispersion, , and Eq. (14) simplifies to which is the eigenvalue equation of the non-dispersive RSE (1); (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 non-zero 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 ( is the orbital number) of spheres made of a dispersive material (gold) and a non-dispersive material (sand, ) 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  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 of RSs for the RSE basis, including only RSs satisfying the condition . This excludes RSs having a wavevector in the medium above , which is the case of large or large with 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 , demonstrating a high accuracy given the strong perturbation. For the present geometry, is approximately proportional to , with for eV. The observed 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 has no dispersion, the only valid sum rule in Eq. (8) is the one with which provides only one alternative GF representation with in Eq. (11). Replacing by in Eq. (12) results in a nonlinear dispersive RSE

 2(ωn−ω)cn=ω∑mVnm(ω)cm, (16)

which is a direct generalization of the original, non-dispersive RSE (1); (2). For a finite number of poles in the RP, Eq. (16) can be written as a polynomial matrix equation. The order of the polynomial is given by the number of non-zero poles in Eq. (2), so for example, for the Ohm’s law model (linear matrix problem), for the Drude model (quadratic problem), and for the Drude-Lorentz model with 2 pairs of CPs. For any finite , such a polynomial eigenvalue problem can be solved by linearization (17), extending the basis of unperturbed RSs by a factor of .

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

###### Acknowledgements.
This work was supported by the Cardiff University EPSRC 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.

### References

1. E. A. Muljarov, W. Langbein, and R. Zimmermann, Europhys. Lett. 92, 50010 (2010).
2. M. B. Doost, W. Langbein, and E. A. Muljarov, Phys. Rev. A 90, 013834 (2014).
3. C. Sauvan, J. P. Hugonin, I. S. Maksymov, and P. Lalanne, Phys. Rev. Lett. 110, 237401 (2013).
4. Q. Bai et al., Opt. Express 21, 27371 (2013).
5. E. Muljarov, M. Doost, and W. Langbein, arXiv 1409.6877 (2014).
6. M. B. Doost, W. Langbein, and E. A. Muljarov, Phys. Rev. A 85, 023835 (2012).
7. M. B. Doost, W. Langbein, and E. A. Muljarov, Phys. Rev. A 87, 043827 (2013).
8. L. J. Armitage, M. B. Doost, W. Langbein, and E. A. Muljarov, Phys. Rev. A 89, 053832 (2014).
9. M. Doost, W. Langbein, and E. Muljarov, arXiv 1508.03851v1 (2015).
10. P. B. Johnson and R. W. Christy, Phys. Rev. B 6, 4370 (1972).
11. A. D. Rakić, A. B. Djuriić, J. M. Elazar, and M. L. Majewski, Appl. Opt. 37, 5271 (1998).
12. P. G. Etchegoin, E. C. Le Ru, and M. Meyer, J. Chem. Phys. 125, 164705 (2006).
13. P. G. Etchegoin, E. C. Le Ru, and M. Meyer, J. Chem. Phys. 127, 189901 (2007).
14. A. Vial and T. Laroche, J. Phys. D: Appl. Phys. 40, 7152 (2007).
15. L. D. Landau, L. P. Pitaevskii, and E. Lifshitz, Electrodynamics of Continuous Media, Vol. 8 of Course of Theoretical Physics, 2nd ed. (Elsevier Butterworth Heinemann, Oxford, 1984), iSBN-10: 0750626348.
16. W. Broer and B. J. Hoenders, J. Phys. A: Math. Theor. 42, 245207 (2009).
17. J. Higham, D. S. Mackey, N. Mackey, and F. Tisseur, SIAM J. Matrix Anal. Appl. 29, 143 (2006).
100960