On the pole expansion of electromagnetic fields

In several publications, it has been shown how to calculate the nearor far-field properties for a given source or incident field using the resonant states, also known as quasi-normal modes. As previously noted, this pole expansion is not unique, and there exist many equivalent formulations with dispersive expansion coefficients. Here, we approach the pole expansion of the electromagnetic fields using the Mittag-Leffler theorem and obtain another set of formulations with constant weight factors for each pole. We compare the performance and applicability of these formulations using analytical and numerical examples. It turns out that the accuracy of all approaches is rather comparable with a slightly better global convergence of the approach based on a formulation with dispersive expansion coefficients. However, other expansions can be superior locally and are typically faster. Our work will help with selecting appropriate formulations for an efficient description of the electromagnetic response in terms of the resonant states. © 2020 Optical Society of America under the terms of the OSA Open Access Publishing Agreement

Instead of using full-wave simulations, it has become quite popular to develop semi-analytical descriptions based on the resonant states, also known as quasi-normal modes [17,18]. The advantage of this approach is that it requires calculating only a finite and usually small number of resonant states of a system in the spectral region of interest, and to use this physically meaningful set of modes in order to describe both local excitation [19][20][21][22] and light scattering [23][24][25][26][27]. Several non-intuitive predictions of that theory have been validated recently [28,29]. If the number of resonant states is small, it is possible to derive semi-analytical expressions, e.g., for the performance of refractive index sensors [30][31][32][33][34]. Of particular interest is the interaction with quantum emitters and the calculation of related quantities such as mode volume, quality factor, and Purcell enhancement [19][20][21]35], which has been recently extended to a quantum-optical description via a quantization of the resonant states [36]. When using a large number of resonant states, an eigenvalue equation can be set up in order to derive the resonant states of a perturbed system via the so-called resonant-state expansion [30,37], where the perturbation can be extremely large. For instance, the resonant states of a metallic sphere can be calculated from those of a dielectric sphere [38]. Completeness is of course an important issue [39], and it might require including so-called static modes [40]. The resonant-state expansion has been also extended to guiding geometries [41][42][43][44], and so-called permittivity eigenmodes [45].
Several numerical approaches exist to calculate the resonant states, while analytical solutions are only known for few geometries such as planar slabs, see Fig. 1. An interested reader might consult the benchmark of different numerical methods in [46]. The principle concept is to search for solutions of Maxwell's equations with outgoing boundary conditions in the absence of sources.
However, the resulting field distributions are only determined up to a scalar complex number, so that the correct normalization for the resonant states has to be found by other means. Only recently, several approaches have been developed, including approximate far-field schemes [47], numerical schemes [32,48], integration over complex coordinates or perfectly-matched layers [20], and fully analytical approaches [18,37]. Fig. 1. Schematic of a homogeneous slab. The top panel displays transmittance spectra for a GaAs slab of thickness d = 1 µm and a refractive index of 3.5 in air at normal incidence from the top. Below the spectra, we depict the complex energy plane, with each resonant state corresponding to a pole marked by a red cross. The resonant field distributions inside the slab are indicated by red lines at the front of the slab.
Using an appropriate normalization scheme, it is then possible to expand the Green's dyadic inside a resonator geometry in terms of the resonant states [18], with the complex resonance frequencies as poles of the Green's dyadic. Thus, the local electromagnetic fields can be determined for a given source in terms of the resonant states. As noted in [17], this resonant expansion is not unique, and several equivalent approaches exist to calculate the pole contributions. Recently, the performance of various approaches has been compared numerically for Lorentzdispersive materials [49], with only minor differences between them. Furthermore, it has been shown that there exists in principle a continuous family of possible expansions [50]. However, all these formulations contain frequency-dependent weight functions for each pole. These weight functions consist of overlap integrals of resonant field distributions with the source terms.
Here, we approach the question of suitable field expansions via the Mittag-Leffler theorem. Thus, we arrive at formulations equivalent to those in [49] for arbitrary materials. As an extension to previous work, we obtain formulations that contain weight functions evaluated solely at the complex resonance frequencies. Hence, we do not have to evaluate overlap integrals repeatedly for each frequency of interest. The drawback of these formulations is a more complex background contribution. The performance of the different formulations is compared for the analytic example of a planar homogeneous and isotropic slab in air as well as for the numerical examples of a one-dimensional dielectric and a two-dimensional metallic grating.

Theoretical approach
For the sake of brevity, we formulate Maxwell's equations in frequency domain [time dependence exp(−iωt), Gauss units, wavenumber k = ω/c] in terms of a matrix-operator formulation [18]: Here,M = kP −D, witĥ where ε, µ, ζ, and ξ are the permittivity, the permeability, and possible bi-anisotropic contributions, respectively. The six-dimensional supervector F comprises electric and magnetic fields: If we consider a region in space that is free of sources, i.e., the right-hand side of Eq. (1) is zero, the electromagnetic fields usually consist of incoming and outgoing parts. As described in [17,27], the incoming field F in may be included for a scattering geometry with background material distributionP BG and local material change ∆P =P −P BG as the incoming part of a background field F BG = F in + F out BG , which is a solution of the following Maxwell's equations: whereM BG = kP BG −D. Thus, the total field supervector can be expressed as a superposition of the background and a scattered field as F tot = F BG + F scat . Although we consider here only homogeneous and isotropic background material distributionsP BG , it is in general possible to introduce more complex background systems such as planar interfaces between two materials to account for a different substrate material. After some algebra, we obtain [17,27]: In the latter equation, the second term is equivalent to a source current J scat emitting the scattered field. Knowing the Green's dyadicĜ of Eq. (1) witĥ where I is a six-dimensional unit matrix, we can calculate the scattered field as Based on the Mittag-Leffler theorem [51], the Green's dyadic can be expanded inside a resonator in terms of the resonant states [18]: These resonant states constitute a discrete number of solutions of Maxwell's equations for J = 0 at complex wavenumbers k that satisfy purely outgoing boundary conditions: This equation defines the resonant field distributions F n only up to a complex scalar. In order to use them in Eq. (8), an appropriate normalization scheme has to be applied. As mentioned in the introduction, several approaches have been developed [18,20,32,37,47,48]. We are using here the analytic approach described in [18]. Note that we neglect in Eq. (8) possible cut contributions [52] and the superscript R refers to the reciprocal conjugate of the electric and magnetic fields of the resonant states [53,54]. Reciprocal conjugate resonant states are solutions at the same resonant energies but for reciprocal boundary conditions. For planar geometries, reciprocal boundary conditions are given by an inversion of the inplane wavevector components. Combining Eq. (8) with Eq. (7), we thus can construct the total field as where we have introduced the overlap integral I n , which is given in the absence of bi-anisotropic contributions (ζ = ξ = 0) by with ∆ε = ε scat − ε BG and ∆µ = µ scat − µ BG . The latter terms are in general k dependent 3 × 3 tensors describing the difference of permittivity and permeability between the background material distribution and the scatterer. The pole expansion of the total field given by Eq. (10) is broadly used as a semi-analytical method to expand the near fields for a given resonator system [17,27,48]. Henceforth, we refer to it as formulation with dispersive expansion coefficients. When approaching the resonant expansion of the Green's dyadic via the wave equation for the electric fields for µ = 1, a slightly different expansion of the field is obtained: This formulation is used, e.g., in [23,55]. The only difference is the factor k/k n in the pole contributions. Using the sum relation n E n ⊗ E R n /k n = 0 [21], the equivalence of Eq. (10) and (12) can be shown when considering the complete set of resonant states and a constant µ = 1. Both formulations are compared in [49] for a finite number of resonant states, with a negligible advantage of the expansion via Eq. (12).
It should be noted that the integrals I n (k) in Eqs. (10) and (12) have to be evaluated at every wavenumber k. This complicates the application of this type of expansion, since overlap integrals have to be calculated repeatedly. To circumvent this drawback, it is possible to consider the total field as an analytic function of k with a countable number of poles and to apply the Mittag-Leffler theorem to the total field in order to simplify the k dependence of Eq. (10). However, the application of the Mittag-Leffler expansion to the field is in general not as simple as that for the Green's dyadic. Depending on the asymptotic behavior of the field for k → ∞, a suitable order of the Mittag-Leffler expansion has to be chosen.
In general, the Mittag-Leffler theorem states that a complex function f (z) that is analytic except for a countable number of poles a n with residues b n and the asymptotic behavior lim can be expanded as [51] Here, f (m) (0) denotes the m th -order derivative of f at z = 0. For p = 0, we obtain which can be used for the Green's dyadic, while we obtain for the first-order (p = 1): In both cases, a rather simple pole contribution as a sum over b n /(z − a n ) arises, with a nontrivial background f (z) + n b n /a n for the first-order Mittag-Leffler expansion. In general, a p th -order Mittag-Leffler expansion can be reformulated as That means that for any order p, we obtain a pole contribution of the form n b n /(z − a n ), accompanied for p>0 by a polynomial of order p − 1.
If the p th -order Mittag-Leffler expansion is applicable, an expansion of order p with p >p is possible, too. Assuming that the Mittag-Leffler expansion is unique, we thus infer sum rules which must hold for any order m>p. In principle, the sum rules in Eq. (17) could be used in order to estimate the accuracy of a pole expansion when considering only a finite set of resonant states. However, this is beyond the scope of this work. Furthermore, it should be noted that a finite number of poles in the expansion results in an error that can be compensated locally by a polynomial correction. For instance, in Ref. [27], a polynomial fit has been carried out to approximate the optical scattering matrix based on a finite number of resonant states. From the above analysis, we deduce that it is irrelevant if the resulting polynomial correction originates from a nonresonant background or from neglecting spectrally-distant resonant states.
In the case that the background field is free of poles, the p th -order Mittag-Leffler expansion of the total field yields where F (m) tot denotes the m th derivative of F with respect to k. Evidently, the pole contribution in Eq. (18) is much simpler than that in Eq. (10), since it requires calculating the overlap integral in Eq. (11) only once at the complex wavenumbers k n of the poles. The drawback is that we need to account for a more complex background. Henceforth, we refer to the pole expansion based on Eq. (18) as the p th -order Mittag-Leffler expansion of the total field, with p = 0, 1 in our examples.
It should be noted that the same arguments hold for the scattering matrix, too. More specifically, in [26], a k-dependent numerator occurs in the pole expansion of the scattering matrix that is accompanied by a trivial background term. In [27], an alternative derivation leads to a k-independent numerator in the pole expansion at the cost that the background term of the scattering matrix becomes less trivial. In some cases, a combination of trivial background and pole contributions might be assumed [24], but this does not hold in general.
As mentioned above, the pole expansion of the Green's dyadic is only valid inside the resonator. Outside, the field divergence of the resonant states towards the far-field region prevents an accurate description of the electromagnetic response. However, resonant states can be regularized in the exterior [56,57]. Similarly, reciprocity principle can be used to derive the analytic continuation of the near to the far fields [27,58]. Here, we are using the latter approach to derive the pole expansion of the electromagnetic fields in the far-field region. More specifically, if {I N } and {O N } constitute sets of incoming and outgoing basis functions, respectively, that are solutions of Maxwell's equations in the exterior of the resonator with N as the quantum number to distinguish the different modes, we may expand any total field in terms of these basis functions as F = N α in N I N + α out N O N . Via reciprocity and orthogonality, we can then derive sets of reciprocal-conjugate probe functions {I R N } such that Details on the derivation of this relation can be found in [27], including appropriate basis functions for certain geometries. In summary, by knowing the fields at the outermost interfaces of a given resonator structure, e.g., based on the pole expansion of the fields via Eqs. (10) or (18), we can thus construct the scattered field in the exterior at arbitrary wavenumbers. Note that when using Eqs. (18) and (19), we fix the phase of the incoming and outgoing basis functions to be zero across the normal direction of ∂V, since Eq. (19) must result in the decomposition of a resonant state in terms of outgoing basis functions when approaching its complex wavenumber.

Analytic example
Let us now consider the simple example of a planar isotropic and homogeneous dielectric slab with refractive index n 2 in a symmetric environment with refractive index n 1 . At normal incidence, the Fresnel equations for the amplitudes of reflection ρ and transmission τ at the top interface of the slab (from medium 1 to 2) are given by ρ = (n 1 − n 2 )/(n 1 + n 2 ) and τ = 2n 1 /(n 1 + n 2 ), respectively. In that case, the total electric field is aligned parallel to the slab. Without the loss of generality, we chose x-polarized fields with E tot = E tot x, where x is the unit vector along x direction. Thus, fields inside and outside the structure can be written as where E 0 is the amplitude of the incidence plane wave. Note that our z axis is oriented from top to bottom parallel to the wavevector of a normally incident field, as indicated in Fig. 1. The coefficients r and t are the reflection and transmission amplitudes of the dielectric slab, with r(k) = ρ e 2in 2 kd − 1 1 − ρ 2 e 2in 2 kd , t(k) = τ 2 n 2 n 1 e ikn 2 d 1 − ρ 2 e 2in 2 kd .
The complex wavenumbers and internal field distributions of the resonant states in the slab are given at normal incidence by [17,59] where n is an integer number. For E BG (z; k) = E 0 exp[ikn 1 (z + d/2)], the overlap integral in Eq. (11) yields: Here, sinc x = sin x/x. By employing Eq. (24) in Eq. (10), the total field can thus be expressed by the formulation with dispersive expansion coefficients.
For the Mittag-Leffler expansion of the total field, we need to investigate the asymptotic behavior of Eq. (20). We can distinguish two cases: The total field has a non-vanishing upper bound at the position of the top interface, so that the global behavior of the fields can be captured only in a first-order Mittag-Leffler expansion. This agrees with [39], where the pole expansion of the fields in a slab is derived by a different approach. In order to obtain the final Mittag-Leffler expansion of the total field, we have to evaluate now Eq. (24) at k n and insert it together with Eq. (23) and E tot (z; k = 0) = E 0 into Eq. (18). In the columns, the pole expansion of the field has been calculated using the formulation with dispersive expansion coefficients (left) as well as the first-(middle) and zeroth-order (right) Mittag-Leffler expansion of the total field.
As test system, we consider a 1 µm thick slab made of GaAs (n 2 = 3.5) and surrounded by air. Figure 2 Expectedly, the agreement between analytic results and the different expansions becomes better when increasing the number of poles. Almost no deviations can be seen inside the slab for 33 poles. Furthermore, we observe that the zeroth-order expansion of the total field does not agree at the top half space. This is not surprising, because the zeroth-order Mittag-Leffler expansion is not applicable at the top interface due to the asymptotic behavior at this interface, see Eq. (25a).
Let us examine the differences between the three approaches in more detail. For that purpose, we calculated the deviation between the pole expansions and the analytic results inside the slab at 650 meV and display the results in Fig. 3(a) in a logarithmic scale. As expected, the largest deviations are observed at the top interface when using the inappropriate zeroth-order expansion of the total field (blue solid line). The fields based on the dispersive expansion coefficients (green solid line) converge best at the top interface, but exhibit larger deviations at the bottom interface, where the first-order results (red dotted line) agree best.  Figure 3(b) contains the Fourier transform of the deviation between exact results and their pole approximation, which is normalized by the Fourier transform of the analytic field. We observe a similar behavior for zeroth-and first-order expansion of the fields. The agreement to exact results is best for small spatial frequencies. For large spatial frequencies, there is a saturation, with the deviation being larger for the zeroth-order expansion. The reason for the saturation is that the finite number of resonant states in the pole expansion captures only a finite number of low spatial harmonics. For the results obtained by the formulation with dispersive expansion coefficients, the agreement also first decreases up to the spatial frequencies captured by the corresponding number of poles, but it then decreases significantly to very low values. A possible reason is that the k-dependent background field allows for a better description of the background contributions to the fields.
Finally, we plot in Fig. 3(c) the convergence of the different approaches in dependence of the number of poles for the resonant expansion. More specifically, we calculate the L 2 error as with the integration being carried out over the entire slab region. It can be seen that the results obtained by the formulation with dispersive expansion coefficients converge best, followed by the first-order Mittag-Leffler expansion of the total field. The worst convergence is obtained when using the zeroth-order Mittag-Leffler expansion of the total field.

Numerical examples of periodic structures
As first numerical example, we consider a one-dimensional dielectric grating with a periodicity of 0.8 µm that consists of rectangular ZnSe rods in air with a refractive index of 2.6, a width of 0.53 µm, and a thickness of 0.1 µm. A schematic is shown in Fig. 4. Henceforth, we investigate s-polarized incident plane waves with in-plane wave vector components k x = 0.5 µm −1 and k y = 0.7 µm −1 . Figure 5(a) displays the transmitted (blue) and reflected (red) intensity as a function of photon energy. The numerical results have been obtained using our in-house implementation of the Fourier modal method [60,61]. In the given range, the system exhibits three resonant states at complex eigenenergies 1892 − 15i meV, 2027 − 150i meV, and 2370 − 29i meV. Their resonant field distributions are displayed in Fig. 4(b). Based on these resonant states, we then construct the pole expansion of the fields via the formulation with dispersive coefficients and the Mittag-Leffler expansion of the total field. For the latter, we consider only the first-order Mittag-Leffler expansion, since we have shown already that a zeroth-order expansion is in general inappropriate. Furthermore, we do not calculate the background field in any expansion exactly, but we fit a polynomial of quadratic order to the deviation between the pole contributions of the three resonant states and numerically exact results at energies of 1700 meV, 2050 meV, and 2400 meV. Thus, we capture both the impact of the background field as well as spectrally distant poles. The latter is particularly necessary even for the formulation with dispersive expansion coefficients, because there are several resonances at lower energies as well as the lowest-order Rayleigh anomaly that we do not account for in our example. The relative difference of both Mittag-Leffler expansions to the numerically exact results is displayed in Fig. 5(b). More specifically, we calculate the L 2 error according to Eq. (26) by integrating over one unit cell of the structure along x and over the thickness of the grating plus 50 nm of superstrate and substrate. Then, we normalize this absolute error by the integral over |E ref | 2 of the numerically exact results. It can be seen that we achieve an overall good agreement below 1 %. As expected, the deviation even drops to zero at the three energies of the background fit. The difference between the two pole expansions is negligible. Panels (c-h) in Fig. 5 contain field magnitudes at energies of 1870 meV (left column) and 2225 meV (right column), which are close to the two resonant states at 1892 − 15i meV and 2370 − 29i meV. The rows display (from top to bottom) results based on full-wave simulations, the first-order Mittag-Leffler expansion of the total field, and the formulation with dispersive expansion coefficients. The small differences are hardly visible, while the calculation time for obtaining 150 field distributions in the spectral range of Fig. 5(a) is twice as long with the full numerical simulations as with the pole expansion approaches in our implementation, with the expansion based on constant coefficients being slightly faster.
As second numerical example, we consider a two-dimensional metallic grating consisting of periodically-arranged rectangular holes filled with air inside a slab made of gold, see Fig. 6(a). The periodicity of the structure is 1.2 µm in x and y direction, and the slab is 50 nm thick. The holes are 300 nm by 960 nm wide. The refractive index of the surrounding is assumed to be 1, the permittivity of gold is calculated using the critical point model [62]. For more complex materials, a pole expansion can be fitted to tabulated data [63]. In the following, we account for one resonant state at the complex eigenenergy 543 − 135i meV, see panel (c) in Fig. 6. The structure is excited from the top by a plane wave at normal incidence with the electric field polarized along the x direction. Figure 6(b) displays the transmittance (blue), reflectance (red), and absorbance (black) as a function of photon energy. As in the previous example, the numerical results have been obtained using the Fourier modal method, including the field distribution in Fig. 6(d), which is derived at 537 meV inside the structure 12 nm below the top interface.
Using the aforementioned resonant state, we then calculated the electric field by the formulation with dispersive expansion coefficients (e) and the first-order Mittag-Leffler expansion of the total field (f). The background field and the spectrally distant poles have been fitted via polynomials of first order accounting for the deviation between the pole contributions of the resonant states and numerically exact results at energies of 300 meV and 800 meV. The fastest approach is based on the Mittag-Leffler expansion of the full field, but the fields based on the formulation with dispersive expansion coefficients are more accurate, as can be seen by the small deviations between panels (e) and (f) of Fig. 6.
Finally, it should be noted that we do not include any results of Eq. (12) here, because those are comparable to that of Eq. (10) for the numerical examples, while Eq. (12) does not work properly for the dielectric slab. This is probably due to a slow convergence of the sum rule mentioned after Eq. (12).

Conclusion
We have compared different formulations of the pole expansion of electromagnetic near fields for the analytical example of a homogeneous and isotropic slab as well as for the numerical examples of a dielectric and a metallic grating. The most common formulation based on the formulation with dispersive expansion coefficients yields globally the best agreement. Locally, the first-order Mittag-Leffler expansion of the total field, which we developed here, can be superior. Moreover, this approach includes constant weight factors for the pole contribution, so that it is typically faster. This means that the formulation based on the first-order Mittag-Leffler expansion is preferable whenever it is necessary to calculate the spectral response at many different wavenumbers. However, a zeroth-order expansion of the total field is inappropriate for the pole expansion of the fields.