A continuous family of Exact Dispersive Quasi-Normal Modal (DQNM) Expansions for dispersive photonic structures

In photonics, Dispersive Quasi-Normal Modes (DQNMs) refer to optical resonant modes, solutions of spectral problems associated with Maxwell's equations for open photonic structures involving dispersive media. Since these DQNMs are the constituents determining optical responses, studying DQNM expansion formalisms is the key to model the physical properties of a considered system. In this paper, we emphasize the non-uniqueness of the expansions related to the over-completeness of the set of modes and discuss a family of DQNM expansions depending on continuous parameters that can be freely chosen. These expansions can be applied to dispersive, anisotropic, and even non-reciprocal materials. As an example, we particularly demonstrate the modal analysis on a 2-D scattering model where the permittivity of a silicon object is drawn directly from actual measurement data.


Introduction
In photonics, the interaction between the electromagnetic field (light) and matter heavily relies on the concept of resonant modes of the structure, the privileged vibrational states of the optical system [1,2]. From a mathematical point of view, such resonant modes correspond to solutions of source-free Maxwell's equations, called Dispersive Quasi-Normal Modes (DQNMs). Under external excitation, these modes are initially loaded, then release their energy which, in turn, reflects the optical properties of the system. Therefore, it is believed that we can take advantage of these resonances (DQNMs) as building blocks in order to construct the physical characteristics of the given system. Since establishing the spectral representation of the diffracted field on a set of resonance-state basis can lead to a transparent interpretation of the numerical results, modal expansion formalisms have recently received a lot of attention [3][4][5][6][7].
The name Dispersive Quasi-Normal Modes (DQNMs) is derived from the fact that in practice, our media are highly dispersive and optical structures are located in unbounded media. As a result, computing these 'resonant' DQNMs implies solving complicated non-Hermitian (i.e. with complex eigenfrequencies) non-linear eigenvalue problems [8]. However, numerical computations of these non-linear eigenvalue problems can still be done thanks to the current development of a powerful software library in SLEPc [9]. Unfortunately, the final step to rigorously formalize modal expansions using these non-linear eigenmodes has still been a difficult mathematical task debated in recent literature [1,7].
In an earlier paper [10], we have introduced a simple approach where the exact DQNM expansion is based on the Keldysh theorem [11,12] and the concept of eigentriplets: The non-Hermitian property of dissipative systems requires the introduction of a pair of dual eigenvectors, i.e. the 'left' and 'right' eigenvectors associated to the same eigenvalue. In particular, we have shown the accuracy of our modal expansion on a wide frequency range in a very strongly dispersive case where the system is closed and the permittivity is given by an artificial Lorentz model. In the present work, we report a new insight on spectral theory when applied to Maxwell's equations: There exists not just one, but a continuous family of exact DQNM expansions (in particular formulas Eq. (9) and Eq. (14)). In order to shed light on this rather strange fact, we give an example where several expansions in the family are tested for both electric and magnetic fields. Not limited to a closed system as the previous paper, we also extend our computation to open structures. Moreover, the permittivity in our calculations is no longer limited to simple models (for example the Drude-Lorentz model) but is drawn directly from actual measurement data [13], which can open up many more possibilities in practical applications. The outline of the paper goes here: Section 2 introduces the general concepts of modal expansion and rational operators. It also reveals a very important spectral property of rational operators: There exists a continuous family of modal expansion formulas. Then, in section 3, we will derive the exact DQNM expansions for both electric and magnetic fields and look more closely at the non-uniqueness of these expansions. Finally, we illustrate some numerical examples in section 4 showing the effectiveness of our expansion formalism in applications with realistic materials.
It is worth noticing that the permittivity and permeability in the previous equations are expressed by bold notations ε ε ε(λ) and µ µ µ(λ) representing 'second-order tensors'. This implies the fact that there is no special restriction on the materials: they can be dispersive, anisotropic, and non-reciprocal.
The resonant modes, i.e. DQNMs, of our electromagnetic system are indeed eigensolutions of the following spectral problem: where v n are the 'right' eigenvectors associated to eigenvalues λ n . The main goal is to formalize a modal expansion equation which allows us to express the electromagnetic field u in terms of linear combination of v n : u = n a n (S)v n . This requires us to investigate the spectral properties of the operator M ξ ξ ξ,χ χ χ (λ). In particular, we will latter point out that since the time-dispersion of ξ ξ ξ(λ) and χ χ χ(λ), which refer to the permeability µ µ µ(λ) and permittivity ε ε ε(λ) in the case of electric fields and vice versa, can be represented by rational functions, M ξ ξ ξ,χ χ χ (λ) can be understood as a rational operator.

Time-dispersive materials and rational operators
Studying the dissipation of the materials requires introducing the time-dispersion on the permittivity. This can be done by imposing a hypothetical mathematical model to the permittivity, for example, the Lorentz model [14]. However, it is more practical to assume that rational functions are a general representation of the frequency-dependence of the permittivities of media. In particular, given the polarization vector P, the causality principle can be imposed via the constitutive relation between P and electric field E, represented by the equation: It implies that the electric susceptibility χ χ χ e (iω), such that P(ω) = χ χ χ e (iω)E(ω), must be a rational function with respect to iω where all the coefficients are real (The a i and b i may be extended to tensors but with real coefficients); and so must the permittivity: ε ε ε(iω) = ε 0 (I + χ χ χ e (iω)) where I stands for the identity tensor.
The dispersion of permittivity can be obtained through an interpolation method [13] that is very accurate on a large range of frequencies and thrifty with the number of poles. The obtained rational functions are naturally causal (following Kramers-Kronig relations) and provide a natural analytic continuation of permittivities in the complex plane (including negative permittivity regions giving rise to surface plasmons).
Since the permittivity can be given a rational function, it is practically helpful to introduce the concept of rational operators: A rational operator R L (λ) is defined as follows: where L i are constant-coefficient partial differential linear operators acting on the electromagnetic field and R i (λ) stand for rational functions with respect to λ: The numerator and denominator of rational function R i (λ) are polynomials of degree deg(n i ) and deg(d i ) respectively.

Spectral property of rational operators
In this subsection, we expose the spectral property of rational operators, which leads us to the formalism of modal expansion. In this paper, the following convention is used: w, v = ∫ Ω w(r) · v(r)dΩ. The 'bra-ket' notation is also introduced: The 'ket' vector |v denotes a vector such that we can simply write v = |v when there is no ambiguity. On the other hand, the 'bra' w| represents the vector in dual space. The conjugate transpose, i.e. Hermitian conjugate, of a bra is the corresponding ket and vice versa: v| * = |v . In the finite-dimensional vector space, we have that v * = (v) where v denotes the transpose and v stands for the complex conjugate of v. Given a rational operator R L (λ), the eigentriplets (λ n , w n |, |v n ) of R L (λ) are defined as follows: where w n | and |v n are the 'left' and 'right' eigenvectors corresponding to the same eigenvalue λ n . We assume that the operator R L (λ) is diagonalizable and all the eigenvalues λ n are simple.
is the complex derivative of R L (λ) with respect to the complex variable λ. In practice, the complex derivative of operators is obtained by taking the complex derivative of coefficients (functions of λ ∈ C) in the operator, i.e. the rational functions R i in our case.
The rational operator can be multiplied by a N D -degree polynomial D(λ) that is divisible by all the polynomials in the denominators to get the N N -degree polynomial operator N L (λ) = R L (λ)D(λ). The degrees N N and N D are computed by the following equalities: Through the process of linearization, the polynomial operator N L (λ) can be expressed in terms of a system of linear problems. From there, we have proved that the solution u of the non-homogeneous problem R L (λ)u = S, which is also the solution of N L (λ)u = D(λ)S, can be expanded in the form of the following quasi-normal modal expansion formulas: (5)), we obtain:

The exact DQNM expansions for electromagnetic problem
Since M ξ ξ ξ,χ χ χ (λ) can be seen as a rational operator, it is possible to derive the DQNM expansion for an electromagnetic problem with dispersive media using the expansion formula Eq. (8). In order to do that, we firstly have to clarify the meaning of the 'left' eigenvalue problem of M ξ ξ ξ,χ χ χ (λ): w n |M ξ ξ ξ,χ χ χ (λ n ) = 0. For a deeper discussion of the 'left' eigenvalue problem and adjoint operator, we refer the reader to the appendix.

The exact DQNM expansions for electric fields
Equipped with the eigentriplets of the operator M ξ ξ ξ,χ χ χ (λ), our next step is to apply Eq. (8) in order to obtain the DQNM expansion for electric fields. Given the operator M µ µ µ,ε ε ε with the eigentriplets (λ n , E ln |, |E r n ), the modal expansion of the solution E of wave equation M µ µ µ,ε ε ε E = S appears to be: where the inner product in the denominator can be computed explicitly as follows: where the prime denotes the complex derivative of functions: The final step is to identify the value of degree ρ in Eq. (9) by finding out the value of N D and N N . Since the frequency-dependence of permittivity can be efficiently represented by a rational function [13], it is reasonable to consider that the orders of polynomials in the numerator and denominator of such rational functions are equal. Indeed, the permittivity is given by ε ε ε = ε 0 (I + χ χ χ e ) where ε 0 is the electric permittivity of free space. And at very high frequencies, the susceptibility χ χ χ e decreases asymptotically as χ χ χ e ∝ (1/ω 2 )I and ε ε ε → ε 0 I. Thus, we can assume that ε ε ε is given by rational functions whose numerator and denominator are polynomials of the same degree. Then, with the same argument for permeability, Eq. (6) implies that In order to understand the reason behind Eq. (10), let's assume our media are isotropic. Then, our Maxwell operator can be rewritten as follows: where I stands for identity operator. It is easy to recognize M µ µ µ,ε ε ε (λ) in the previous equation as a rational operator with deg(n 0 ) = deg(d 0 ) and deg(n 1 ) = deg(d 1 ) + 2. Using Eq. (6), it is easy to verify that N N − N D = 2. It is worth pointing out that Eq. (10) holds for all dispersive, anisotropic and even non-reciprocal materials.
As a result, f ρ (λ) in Eq. (9) can be any arbitrary polynomial up to degree 1. This means f ρ (λ) would take the form f ρ (λ) = α + λβ with ∀α, β ∈ C and |α| + | β| 0. As a consequence, there exists a continuous family of expansion formulas for the electric field E of the operator For example, if we set f ρ (λ) = 1, the modal expansion is similar to the formalization established using the Keldysh theorem in our earlier paper [10]. On the other hand, with the choice f ρ (λ) = λ, Eq. (9) becomes: which indeed recovers other expansion formulas previously proposed in the literature [1,3]. As an example, we recall equations (4) and (7) from reference [3] using our new notation: where the source S is a dipole located at the point r 0 and f n (ω) is a nonresonant background that is negligible, according to [3]. In the case where the media are reciprocal, i.e. the 'left' eigenvectors are the complex conjugate of their 'right' counterparts (as proved in the appendix), we can prove that Eq. (12) implies Eq. (13) considering λ = iω, S = λJ, and J = −iωpδ(r − r 0 ).

The exact DQNM expansions for magnetic fields
Similarly as in the case of electric fields, a DQNM expansion can be obtained for magnetic fields. Given the operator M ε ε ε,µ µ µ with the eigentriplets (λ n , H ln |, |H r n ), the modal expansion of the solution H of wave equation M ε ε ε,µ µ µ H = S appears to be: where the inner product in the denominator can be computed explicitly as follows: H ln , M ε ε ε,µ µ µ (λ n )H r n = ∫ Ω H ln · ∇ × ((ε ε ε −1 (λ n )) ∇ × H r n ) + H ln · (λ 2 n µ µ µ(λ n )) H r n dΩ.

A continuous family of exact DQNM expansions for electromagnetic wave
Thus far, we have already mentioned the existence of a continuous family of modal expansions for electromagnetic fields, for example, Eq. (9) for electric fields. This property results from the over-completeness of the set of eigensolutions of non-linear in frequency operators in general and Maxwell operators in particular. Still, it is easy to emphasize the counter-intuitiveness of the non-uniqueness of DQNM expansion. For example, it is possible to discard the contribution of any resonant mode from the DQNM expansion of the Maxwell operator.
In order to fully understand the previous property, consider the DQNM expansion for electric fields Eq. (9) with f ρ (λ) = α + λβ, with ∀α, β ∈ C and |α| + | β| 0. Then, among n eigenvalues λ n , let's choose specifically an eigenvalueλ with a corresponding eivenvectorẼ. The issue arises if we choose the function f ρ (λ) such that α = −λβ. Then the factor f ρ (λ n ) will become 0 when λ n =λ. This implies that at the eigenvalueλ, the resonant mode, i.e. the eigenvectorẼ will be excluded in the expansion of the electric field E. This is an extremely unexpected spectral property which has never been discovered by any modal expansion formalism in the literature! Moreover, the non-uniqueness of DQNM expansions also raises a question of which formula of f ρ (λ) we should choose to perform the expansion. We notice that the expansion for electric fields Eq. (9) will explode at the rate 1 λ−λ n when λ moves close to λ n . In practice, this singularity never occurs because λ n are complex resonant frequencies while λ only accepts imaginary values. Unfortunately, the problem will raise if we choose, for example, f ρ (λ) = λ −λ 1 whereλ 1 is in the vicinity of one of resonant modes λ n . Then the expansion Eq. (9) will explode at the rate 1 (λ−λ 1 ) 2 when λ moves close toλ 1 . In fact, we will numerically demonstrate in the next section that the expansion Eq. (9) is no longer accurate in the vicinity of one of the resonant modes λ n if f ρ (λ) = λ −λ 1 .

Numerical examples
We will illustrate numerical results in the geometry of an object shaped like an ellipse Ω 1 inside a perfectly conducting vacuum square Ω 0 (see Fig. 1). The parameters of the structure are chosen in such a way that the material and geometric resonances highly interact with each other. In particular, we want to exhibit the problem where the geometric resonances are in the vicinity of the electric poles of the permittivity.
The elliptic scatterer is made of silicon whose dispersive permittivity is given by the multi-pole 'Lorentz' model: where N p is the total number of poles in our model and ω ε i stands for the complex values of electric poles. For instance, the data of silicon can be extracted from the following table (The frequency unit is ×10 14 (rad.s −1 )): In fact, by adding more poles to Eq. (15) up to N p = 4, we can draw an exact formula of permittivity for any realistic materials from measurement data [13,15]. For example, we plot in Fig. 2 the real and imaginary part of the permittivity of silicon calculated based on Eq. (15) using different numbers of poles. According to Fig. 2, the more number of poles, the better resembling realistic data our model of permittivity. At the same time, the relative permittivity of vacuum is set to be constant ε vac = 1. Then, the whole structure will be illuminated by the Dirac delta source S = δ(r S ) whose coordinates is given by r S = (−2.4, 0.8)(×10 −1 µm), see Fig. 1. The maximum element size is set to be 0.03 µm, in comparison to the smallest wavelength in vacuum 0.2355 µm (equivalent to the highest frequency of the spectrum).

A family of DQNM expansion
The aim of this subsection is to exemplify the non-uniqueness of our DQNM expansion for the geometry in Fig. 1. In particular, we want to compare the results obtained by solving the direct electric scattering problem M µ µ µ,ε ε ε E = S in TE polarization and the field E constructed using DQNM expansion Eq. (9) with different functions f ρ (λ). For the sake of clarification, let's begin with a 1-pole model of permittivity, i.e. N p = 1. The 3-D electrodynamic eigenvalue computations require genuine edge elements to avoid spurious modes but, in our 2-D case, we use longitudinal fields E z e z or H z e z and the associated edge elements reduce to the Lagrange basis elements, here second-order, for the (scalar) component E z or H z .
The eigenvalue problem is numerically solved thanks to recent versions of SLEPc library [9,16] available in the Finite Element Method (FEM) open source package ONELAB/Gmsh/GetDP [17,18] that we are using to implement our models.
The complex resonances are shown in Fig. 3. We emphasize the existence of an accumulation point in the vicinity of the electric pole ω ε 1 (red cross) where ε Si → ∞. The modes around the pole concentrate inside the scatterer and have spatial frequency tending to infinity (for instance mode 2 in Fig. 3), which distinguishes them from conventional modes with eigenfield located in the free-space background (see mode 1 and mode 3 in Fig. 3). It is worth reminding that in the case of closed structures, the spectrum of eigenfrequencies is indeed symmetric through the imaginary axis. Since the contribution of the modes on the left half of the complex plane is numerically insignificant (the factor 1 λ−λ n of these modes is relatively small compared to their counterparts on the right half of the complex plane), in this paper, we do not include them.  Fig. 4); the real and imaginary part of the electric field E 1 calculated at the detector point in Fig. 1 (represented by the middle-left and bottom-left subfigures of Fig. 4). The DQNM expansion in the cases where f ρ (λ) = 1 (blue lines) and f ρ (λ) = λ (red lines) demonstrates excellent agreement with green dots obtained by directly solving the scattering problem except in the vicinity of (ω ε 1 ). In addition, the norm of the field ∫ Ω 1 |E|dΩ reconstructed using the DQNM expansion formula with f ρ (λ) = λ − λ 0 and f ρ (λ) = λ 2 are shown in the top-right subfigures of Fig. 4. The value λ 0 = iω 0 is selected such that ω 0 = 31.628 − 1.478i (×10 14 rad.s −1 ), which is almost coincide with an eigenfrequency. With f ρ (λ) = λ − λ 0 (purple lines), we notice discrepancies around the frequency (ω 0 ), which confirms our prediction in subsection 3.3 about the singularity at the roots of f ρ (λ). As a result, it is recommended to choose f ρ (λ) = α + λβ such that the value −α/β is far away from our domain of interest. Unsurprisingly, when the degree of the polynomial f ρ (λ) is higher than 1, i.e. f ρ (λ) = λ 2 (orange line), we see the less accurate in the numerical results of the DQNM expansion since it is not an appropriate formula.
Finally. it is worth noting that if we choose to derive our DQNM expansion formulas from Eq. (7), the degree of the polynomial g σ (λ) can be set to be larger than 1 (as demonstrated by the bottom-right subfigure of Fig. 4). In particular, the more poles the permittivity model has, the higher the degree of the polynomial g σ (λ) can be.

Multi-pole model of permittivity
This subsection is intended to illustrate the difficulty when increasing the number of poles in the permittivity model Eq. (15). The spectrum of complex resonances with the 4-pole permittivity model N p = 4 is shown in Fig. 5. We especially emphasize the 4 poles (red crosses) which divide the complex plane in several sub-regions. In those sub-regions, we can distinguish two types of eigenfrequencies: Some are distributed separately (for example mode 1, 2, and 4 in Fig. 5) which are responsible for the physical properties of the system; some gather into separate clusters (for example mode 3 in Fig. 5) whose eigenfields oscillate with very high spatial frequencies. These 'cluster' modes, which result from the accumulation points around the poles of the permittivity, affects the performance of the DQNM expansion around the corresponding frequencies (Fig. 6). Indeed, the expansion shows good agreement with respect to the direct data (green dots) except around the frequencies of the 'cluster' modes and the permittivity poles. It is also worth noting that when f ρ (λ) = λ − λ 0 , the expansion Eq. (9) endures discrepancies around (ω 0 ) where ω 0 = 29.633 − 0.295i (×10 14 rad.s −1 ).  To sum up, by raising the number of poles, the permittivity model fits better with realistic materials. On the other hand, a large number of poles will fragment the complex plane of eigenfrequencies and limit the effective frequency range of our DQNM expansion.

Magnetic fields and Plasmonic resonances
The 'cluster' modes around the poles of the permittivity model are not the only problem regarding the numerical implementation of DQNM expansion. In this subsection, we will discover the existence of the plasmonic resonances and their effects on DQNM expansion in the case of magnetic fields in TM polarization.
In particular, for the spectrum of eigenfrequencies of magnetic fields H, there exists the second kind of accumulation points (see sidebar A, B, and C in Fig. 7 for more details), which locate around the plasmon branch where ε(ω 2 ) = −1 (green crosses in Fig. 7). These modes are indeed plasmonic resonances which distribute on the interface between Ω 0 and Ω 1 with the spatial frequencies tending to infinity (for example mode 1, 3, and 4 in Fig. 7). They must be distinguished from the 'cluster' modes caused by permittivity poles (mode 2 in Fig. 7). It is worth pointing out that the locations of these plasmonic resonances do not conspicuously converge (see inset C in Fig. 7) to the analytical position ω 2 where ε(ω 2 ) = −1. In the case of polygonal sign-changing interfaces, the use of structured symmetric mesh stabilizes the numerical discretization of the plasmonic accumulation point [16,19]. However, the symmetry requirements with respect to a curved boundary remain to be elucidated. Unsurprisingly, the plasmonic accumulation points exacerbate the error of the DQNM expansion Eq. (14) in the frequency range nearby them (see Fig. 8).

Open structure and Perfectly Matched Layer (PML)
In the previous computations, we have demonstrated the efficiency of the expansion formalisms in bounded structures. In practice, the electromagnetic system is often open, which makes the leaky resonant modes grow exponentially in space at infinity [20]. A solution is to use the Perfectly Matched Layers (PML) to truncate and damp the electromagnetic fields in free space [21]. Thus, in the last numerical example of this paper, we will mention the open structure as well as the effect of PML modes on the DQNM expansion. The PML layer Ω PML is introduced as seen in Fig. 9. We follow [22] in replacing the initial material properties ε ε ε and µ µ µ in the PML domain Ω PML (vacuum in this case) by equivalent material ε ε ε s and µ µ µ s given by the following rule: where J s is the stretched Jacobian matrix such that: where s x = s y = σ exp(iφ) with σ = 1 and φ = π/10. The regions Ω x PML , Ω y PML , and Ω xy PML are depicted in Fig. 9.
We solve the eigenvalue problem of the electric field in TE polarization and the position of eigenfrequencies is shown in Fig. 10. It is easy to see that the theoretical continuous spectrum, which is supposed to be located on R + , is rotated of an angle θ = −arg(2.8 + 5.2 exp(iφ)) ≈ −0.20456 rad (For detailed instructions of calculation of the angle θ, we refer the reader to [23]). This results in a large number of so-called discretized Bérenger 'PML' modes [21,24], whose eigenfield is concentrated in the PML region Ω PML and far away from the scatterers (as can be seen from the field map of mode 1 in Fig. 10). It should be noticed that the neatly aligned points corresponding to these modes are clearly becoming numerically unstable further on the curve as explained by the pseudo-spectrum theory of L. Trefethen [25]. We can see that the optical properties of the given open structure are fully captured by our DQNM expansion technique at low frequencies (see Fig. 11). When the frequency is larger, the discrepancies of the norm of electric fields inside the scatterer become noticeable since there is almost no field in the domain Ω 1 (Pay attention that the figure is drawn in logarithmic scale and the values are really low). Numerical experiences show that the instability of PML modes (at high frequencies) may add noise to the DQNM expansion. In addition, we notice that the expansion with f ρ = λ provides a better approximation of the field inside the scatterer comparing to the case f ρ = 1. Fig. 11. Electric field obtained by expansion for different functions f ρ = 1 (blue curves) and f ρ = λ(red curves) or by solving a direct problem classically (green dots) for the unbounded structure.

Conclusion
In this paper, the non-uniqueness of the DQNM expansion formalisms is tackled: We systematically discussed the existence of a family of expansion equations, which are determined by the frequency-dependent factor f ρ (λ n ) f ρ (λ) . By modifying this factor, we have been able to recover our previous result in [10] as well as other expansion formulas proposed in [3]. Three examples of expansion formulas of the family are verified by numerical computation. They suggest that the function f ρ (λ) should be chosen in such a way that its roots do not overlap with our domain of interest. Moreover, we also demonstrate the possibility to use the multi-pole model of permittivity to impose the dispersion of realistic materials in the calculation. Even with great potential in practice, we still have to be extremely careful when applying modal expansion with this multi-pole model, since some discrepancies may occur in some particular frequency ranges: • In the vicinity of 'cluster' modes, as the consequence of poles of the permittivity model.
• At high frequencies due to the numerical instability of PML modes as explained by the pseudo-spectrum theory.
Finding a way to overcome these mentioned difficulties will guide the next development of the research on DQNM expansion. Meanwhile, the work will be extended to a wider range of materials as well as 3-D structures in the future.

Funding
This research was supported by ANR RESONANCE project, grant ANR-16-CE24-0013 of the French Agence Nationale de la Recherche.