Quasi-normal mode theory of the scattering matrix, enforcing fundamental constraints for truncated expansions

We develop a quasi-normal mode theory (QNMT) to calculate a system's scattering $S$ matrix, simultaneously satisfying both energy conservation and reciprocity even for a small truncated set of resonances. It is a practical reduced-order (few-parameter) model based on the resonant frequencies and constant mode-to-port coupling coefficients, easily computed from an eigensolver without the need for QNM normalization. Furthermore, we show how low-$Q$ modes can be separated into an effective slowly varying background response $C$, giving an additional approximate formula for $S$, which is useful to describe general Fano-resonant phenomena. We demonstrate our formulation for both normal and fixed-angle oblique plane-wave incidence on various electromagnetic metasurfaces.


I. INTRODUCTION
Scattering phenomena in all areas of wave physics are well described by the universal S-matrix operator. As the resonant (quasi-normal) modes (QNMs) of a system heavily determine its scattering response and coincide with the poles of S, numerous works [1][2][3][4][5] have focused on expressing S as an expansion over QNMs, calculated via eigensolvers. In this paper, we present a QNM theory (QNMT) for multiport lossless scatterers that simultaneously satisfies all fundamental physical constraints of reciprocity, energy conservation and time-domain realness even for the practical case of a small truncated QNM set (in contrast to previous formulations [3,4]) and without the need for the intricate normalization of the divergent QNMs [4,5]. Weak absorption or gain can then be easily incorporated as a perturbation. Furthermore, by explicitly separating a slowly varying effectivebackground response C, we provide a novel additional formula for S, approximate but very convenient to design Fano-scattering systems [6] such as even-order elliptic filters [7]. This C is calculated without resorting to any type of fitting [4,8,9] and without having to choose a specific background scattering medium [5]: We simply use a subset of low-Q modes of the entire system. We then build useful intuition for how various low-Q-mode configurations shape the background C. We demonstrate the accuracy of our QNMT for plane-wave incidence on several electromagnetic (microwave and photonic, 2-port, and 4-port) metasurfaces. In particular, we solve a nonlinear eigenproblem with a complex Bloch wave vector to calculate S with QNMT for a fixed angle of incidence (instead of a fixed transverse wave vector [10,11]).
The resonant modes of open physical systems are often called "quasi-normal modes" (QNM), as they are not square-integrable and exponentially diverge outside the resonator. The QNM eigenvalues correspond to the poles of physical quantities describing the system response. Based on a pole expansion of a desired such quantity, QNMT is usually concerned with identifying the pole residues/coefficients and any additional background terms. These expansions can allow a fast approximate solution for scattering and emission problems, while providing physical understanding and good spectral accuracy around sharp resonances, in contrast to direct numerical methods using frequency/time discretization [12][13][14][15][16]. A quantity often studied is the Green's function [17], which is an infinite-dimensional operator that can be used to construct any solution of the system. However, for problems with a finite number of scattering ports, the scattering S matrix is usually the desired system descriptor, whose calculation is the goal of this paper. It is a simpler finite-dimensional operator which can be computed without requiring the full Green's function. For lossless 1-port systems, the numerator of the scalar S is trivial, since its zeros simply coincide with the conjugates of the poles [18], while loss can be simply treated either by perturbation or by directly computing the zeros [1,19]. In the general multiport case, one recent approach used the full field equations to compute frequency-dependent S-expansion coefficients as volume integrals involving the QNMs and the excitation port fields to achieve good accuracy [5]. Other QNMT formulations have shown how to project the Green's function onto the scattering ports to obtain the S expansion, which is further given as a reducedorder model with frequency-independent residues, an advantage for simplicity and easier interpretation [2,4]. Most of these derivations based only on field-equations solutions require the QNMs' normalization, which can be accomplished only by intricate techniques with increased computational complexity due to the QNM farfield divergence [20]. To avoid normalization, a new phenomenological approach was proposed in Ref. 3, starting from the coupled mode theory (CMT) equations [8,21] and changing basis to the QNMs (a rather confusing approach, as the uncoupled orthonormalized modes required by CMT are ambiguous for arbitrary scatterers). However, for lossless reciprocal systems, these existing formulations do not guarantee energy conservation for a small truncated expansion, but presumably only in the infinite limit. While QNMTs with frequency-dependent expansion coefficients are expected to converge faster towards satisfying this important physical constraint, those with constant residues exhibit large errors for a practical small number of QNMs (Ref. 4 mentions that they violate energy conservation visibly even with 301 modes and we show examples where Ref. 3 violates it by 50% with few modes). In this paper, we consider an S-matrix expansion over the QNMs and directly derive conditions for it to satisfy the necessary physical constraints. We calculate the QNM-to-ports coupling matrix D from simple surface integrals of the fields without need for QNMnormalization, and then finetune D to prioritize and impose these conditions. In this way, rather than the expansion parameters being fixed from the field equations (e.g., if calculated via the Green's function), they are adjusted as more modes are included in order to enforce reciprocity and energy conservation for any finite sum. The final result is a simple equation for S [Eq. (7)] using only the eigenfrequencies and the finetuned D [Eq. (11)] (Sec. II). We confirm the improved accuracy of our QNMT using 2-port and 4-port electromagnetic metasurface examples, and with excitation of both normally and obliquely incident plane waves. For the latter, most previous approaches [10,11] imposed a fixed incidence transverse wave vector (k ⊥ = ωsinθ/c), so that the angle θ changed with frequency ω (given the constant wave speed c). Instead, QNMT can be used to compute S(ω) for fixed θ by evaluating the relevant QNMs involving eigenfrequencydependent complex Bloch wave vectors, formulated as a generalized linear eigenproblem in Ref. 22, while we directly solve the nonlinear eigenproblem here (Sec. III).
In addition to providing a fast computational tool, QNMT (like CMT) has the advantage of offering a simple analytical model for gaining physical insight into resonant systems and for designing practical resonant devices. One interesting example is the case of Fanoresonant shapes [6] emerging from the interplay between a high-Q resonance and a slowly varying background response, useful for sensors and filters. This background scattering is usually described by a separate matrix C, which previous works have almost always estimated only by fitting it a posteriori to the total S, either with a polynomial approximation [4] or an effective averaged structure [8,9]. Recently, an exact volume-integral formula was alternatively derived by factoring out a choice of physical background [5] (but may require further development to handle certain boundary conditions, such as perfect electric conductors in electromagnetism). While it is understood that this background is related to the low-Q modes of the actual structure, a detailed system-atic prescription to compute it directly from them and its relation to the final S are lacking. In Sec. IV, we extend our QNMT to non trivial direct-scattering pathways, by showing that a slowly varying C can be calculated with our general recipe using only the system low-Q modes and by placing the high-Q modes into a different matrix S, in order to then obtain a good approximation S =SC [Eq. (15)]. We then analyze simple low-Q pole configurations corresponding to different physical interpretations of C, such as a desired background transmission or group delay. Finally, we demonstrate this additional formulation for the electromagnetic-metasurface examples mentioned above.

A. Formulation
We consider a general scattering problem of an arbitrary linear time-independent scatterer, coupled to incoming (excited) and outgoing (scattered) radiation via several physical linear ports. At the frequency ω of excitation, the scatterer has a total of P "coupling port modes" (CPM) of radiation, which can be either single propagating modes of P different physical ports or several propagating modes of fewer ports (while all other port modes are either evanescent, or of incompatible symmetry, or their coupling is simply too small at ω). Let CPM p propagate with wave vector k p (ω) and field φ p (ω, r) = φ ⊥ p (ω, r ⊥ p )e ikpr p , separable in the propagation (r p ) and transverse (r ⊥ p ⊥ k p ) coordinates. In the most common case of reciprocal lossless physical ports, the CPMs at ω are orthogonal under the standard (conjugated) "power" inner product (a cross-sectional overlap surface integral) and can be normalized to carry unit power φ ⊥ p |φ ⊥ q = δ pq [21]. Then, for the P pairs of incident and scattered CPM waves, if the vectors s + and s − denote, respectively, their amplitudes at specific reference cross sections r p = z p of their physical ports, |s ±p | 2 equals the power carried by the ±p wave, s † ± s ± is the net incident/scattered power, and the system scattering matrix S at these reference cross sections is defined by Since the scattering system is open (coupled to radiation), its Hamiltonian H is non-Hermitian, so it supports a set of resonant modes [with resonant frequencies ω n and fields ψ n (r) namely H (iω n ) ψ n = iω n ψ n ]. Causality and stability [23] (or simply passivity [18]) imply that the system response is analytic in the upper half of the complex ω plane, namely ω n must lie in the lower half plane. ψ n are linearly independent but quasi-normal and nonorthogonal under the standard (conjugated) "energy" inner product (a volume integral). However, when the system is reciprocal, H is complex symmetric, so its QNMs are orthogonal under the non-conjugated inner product {ψ n |ψ l } = 0 for n = l [13,24]. We consider N such modes, whose normalization we leave unspecified, and denote by the diagonal matrix Ω their complex frequencies and by the vector a their amplitudes upon excitation. Moreover, it is usually assumed there are also pathways other than the resonant QNMs for direct scattering of input to output CPM waves, through the background medium without the scatterer, quantified by a separate scattering matrix C [8,21].
The part of the scattered field not due to direct pathways (s Q − = s − − Cs + ) can be written outside the scatterer as an expansion over the complete set of port modes (propagating CPMs and evanescent). QNMT makes the approximation that it can be written also within the volume V of the scatterer as a linear combination of the N QNMs (an assumption also used in CMT [8]): a n (ω)ψ n (r); r ∈ V.
(1) By inserting the second line into the exact equation for the field inside the scatterer and by mode matching the two lines on the cross-section z p where the p port meets the scatterer boundary, one respectively gets the final two QNMT equations, which, with exp (−iωt) notation, takes the form (for a rigorous derivation see, e.g.,, Ref. 5): where and z p − z p is the distance of the p-port reference crosssection from the boundary of the scatterer. This z p crosssection choice for the calculation of the D overlaps is further justified in Sec. IV C. The P × N matrices K and D quantify the couplings of the QNMs to the input and output CPM waves respectively, and they are generally frequency dependent. Equations. (2) with Ω diagonal (also known as state-space representation in diagonal canonical form in circuit theory [25]) constitute the basis of QNMT and their solution for the scattering matrix S is given by Although a general C (ω) was included in Eqs.
(2) and (4) to align with literature, we rely on the presumed completeness of the resonant QNMs to stipulate that each incident CPM wave is scattered to other CPMs only due to resonances. Therefore, for now, we take C as a diagonal phase matrix. Later, we will show that, indeed, FIG. 1. A 2-port scattering system with important QNMT parameters. At excitation frequency ω, the coupling port modes (CPMs) with transverse fields φ ⊥ p have input and output amplitudes respectively s±p, related by the S matrix through s− = Ss+. The open scattering system supports quasi-normal modes (QNMs) with complex frequencies ωn and fields ψn, which have amplitudes an upon excitation. The CPM-to-QNM coupling coefficients are Dpn, with ratios σn = D2n/D1n. Imposing realness, unitarity, and symmetry constraints on S allows us to compute it as a function of only ωn and σn [Eq. (7)]. Optionally, by separating low-Q modes ψ C n , we can also construct a slowly-varying background matrix C, which can give a physical intuition about the scattering response and help in specific scattering designs [7].
low-Q QNMs can be combined to write a general effective background C (ω); in particular, a fully-transmissive C comes from a zero-frequency mode with an infinite radiative rate (Sec. IV B).
Shifts of ports' reference cross sections-When S(ω) is a meromorphic function, it is useful to employ the Weierstrass factorization theorem and arguments from causality to write S = e iτ ω S e iτ ω , where τ is a constant diagonal matrix with real positive elements and S is a "proper" rational function (Appendix A). Eq. (4) then implies that C(ω) = e i2τ ω C , D(ω) = e iτ ω D and K(ω) = e iτ ω K , where C is a diagonal constant phase matrix that can be taken equal to −I (as we justify later and is often used in CMT [21]) without loss of generality, and D , K are now constant matrices. This is the case for CPMs with fields transverse to their direction of propagation (φ p ·k p = 0), such as plane waves or dual-conductor TEM microwave modes, which will be the focus of this paper. They have k p = ω/c p (where c p the wave velocity) and φ ⊥ p (r ⊥ p ) independent of ω, so Eq. (3) suggests that D pn (ω) = D pn e iω(zp−z p )/cp with D pn in fact constant. Thus, in all structures simulated in this paper, we remove this linear phase to compute S , referenced at the new port cross sections z p on the scatterer boundary. In practice, τ pp may be slightly larger from (z p − z p )/c p , adding a small constant group delay just to the phases of S, so it is of no concern for applications dependent only on their amplitudes, such as amplitude filters [7].
Therefore, hereafter we drop the and consider S to be such a proper rational function that can be decomposed as (5) For other types of CPMs, S(ω) may not be meromorphic, for example when higher-order CPMs have a cutoff frequency, which appears as a branch-point [and frequency dependent D (ω), K (ω)]. Similarly to previous QNMT formulations with constant coefficients, such systems are not investigated in this work.
Normalization independence-Recall that measurable physical quantities (such as the S matrix) do not depend on the choice of normalization for the QNM amplitudes a. For example, it is easily seen from the QNMT Eqs. (2) that, for the fixed normalization of s ± , two different sets (a n , K qn , D pn ) and a n , K qn , D pn scale as a n /a n = K qn /K qn = D pn /D pn , hence S pq in Eq. (5) is unchanged. Thus an overall scaling factor can be chosen arbitrarily for each QNM.
For some physical quantities analytically computed via the QNM fields (such as the Green's function), this normalization independence is typically ensured by dividing with the volume-integral norm {ψ n |ψ n } of these fields. However, since they are non-integrable, regularizing this norm is a procedure that adds complexity (e.g., choice of method and dependence on outgoing boundary condition) [20] to the QNMT formulations based only on calculations from the field equations [2,4,5]. In contrast, as highlighted in Ref. 3 and we show also later, this phenomenological (relying on physical constraints) QNMT formulation does not require such norm evaluation, so it is much simpler.

B. Physical constraints
We now proceed by imposing constraints on S, based on the physical properties of the system.
Realness-For many physical systems, real input fields lead to real output fields so that the system response S(t) must also be real. In frequency domain, this realness of S is stated as S * (iω) = S(−iω * ). For such systems, the same relation holds for their Hamiltonian H satisfying H (iω) ψ = iωψ, so every QNM solution (ω n , ψ n ) is paired with another QNM (ω n , ψ n ) = (−ω * n , ψ * n ). Similarly, the CPMs satisfy φ * Equation (3) thus shows that D pn = D * pn . Then, to satisfy S * (iω) = S(−iω * ), Eq. (5) requires also K pn = K * pn . Note, S realness implies that not only poles but also zeros of any S pq appear in pairs (ω o , −ω * o ), and that S pq is a rational function of iω with real coefficients. Note that our QNMT is still applicable to systems that do not satisfy realness (for example [26]), where simply the QNMs included in the S-expansion do not appear in pairs.
Energy conservation-In absence of absorption or gain, energy conservation implies that the S matrix is [21]. In Appendix B, we prove that a necessary and sufficient condition is given by thus S can be written as In Appendix C, we show that this Eq. (6) choice of K satisfies the realness requirement, namely, for D pn = D * pn , we get K pn = K * pn . In Appendix D, we rearrange Eq. (7) to show that S is fully and uniquely determined by the resonant frequencies ω n and the ratios σ r,pn = D pn /D rnn (for some chosen port r n for each mode n). These quantities can be readily calculated using any appropriate eigenmode solver, where D pn is determined by the surface integral in Eq. (3) and the ratios σ r,pn remove the ψ n -normalization-dependent scaling-factor. Therefore, as promised, S in Eq. (7) does not require computing the volume-integral norms of the QNMs. Note also that, as more modes are included, the residue coefficients automatically update themselves through M in order to satisfy energy conservation for the entire set, in contrast to other formulations based on the exact field equations, where these residues are cosntant [4,5]. Finally, Eq. (6) is different from the usual CMT expression of energy conservation D † D = 2 |Im {Ω}|, associated with modes orthonormal under the standard "energy" inner product, which does not hold for QNMs.
Reciprocity-Reciprocity implies that the S matrix is symmetric (S = S t ) [21]. From Eq. (5), we can see that this is equivalent to having and therefore for some arbitrary diagonal matrix Λ with entries λ n [with the only restriction that λ n = λ * n = 0, so that Eq. (8) is compatible with realness], where a specific choice of λ n fixes the ψ n -normalization-dependent scaling-factor.
For any such choice, as mentioned earlier, the numerator of pole n for S pq in Eq. (9) stays the same, however, its calculation requires evaluation/regularization of divergent volume integrals involving the QNM fields (including their norm {ψ n |ψ n }), which we try to avoid here. Note also that, since Λ can be arbitrary, the usually assumed condition of reciprocity K = D (corresponding to the specific normalization choice Λ = I) is not necessarily true.
D-optimization-In order to satisfy both energy conservation and reciprocity (unitarity and symmetry of S), the input coupling coefficients K must satisfy both Eqs. (6) and (8) simultaneously, namely the output coupling coefficients D must satisfy (10) Here again, for each resonance n, either D rnn (for one port r n ) or λ n can be chosen arbitrarily. In practice, this reciprocity condition Eq. (10) is a set of P N equations. Let D c be the coupling coefficients computed from the eigenmode solver. In most cases, as we see in numerical examples, it turns out that D c are very close to satisfying this required condition, but they do not satisfy it perfectly, since the finite set of chosen resonances is not truly complete. This is why in Ref. 3, Eq. (10) was not enforced exactly, rather the N coefficients of Λ were chosen as λ j = [X † X] jj /[X † D] jj with X = D * (M t ) −1 as one way to minimize its error (note [27]), and their final S matrix, which was formulated as in Eq. (9), was reciprocal but not necessarily unitary. Instead, here, we give priority to exactly satisfying these physical properties of the actual system, if we want our model to be a physically realistic and thus, as we will show, a more accurate description. Therefore, we finetune the P N coefficients D by using a constrained optimization procedure, with the goal of exactly satisfying Eq. (10) while staying as close as possible to the computed system performance: where f (D, D c ) is some penalty function to ensure that D stays close to D c (note [28]). An obvious represents the residues of the scattering matrix expansion given in Eq. (7) and {p, q} is a chosen subset of indices. When summing over p = q, the global optimum is reached for D satisfying R pq (D) = R qp (D) = [R pq (D c ) + R qp (D c )] /2, but such a set of D is not guaranteed to exist. When it does exist, directly solving this system of equations has given the best results in our experience with actual 2-port systems.

C. Properties of 2-port systems
Since many common applications of scattering theory involve lossless reciprocal 2-port systems, we discuss some properties of their S matrix.
(13) As a consequence, when all σ n are ±1 (as for a symmetric structure, whose modes must be even or odd [29]), S 11 = S 22 and S 21 = S 12 , so the 2-port system is immediately reciprocal. [Obviously, σ n = ±1 is not necessary for the reciprocity condition Eq. (10) to hold.] Note that energy conservation alone implies that |S 12 | = |S 21 |, so reciprocity in lossless 2-ports is mostly a statement on the transmission phase responses.
Finally, in Appendix E, it is shown for a real lossless reciprocal 2-port system that (i) the zeros of S 21 can only appear as complex quadruplets is a zero pair of S 22 . The restrictions (i) on the S 21 zeros imply that its numerator is a polynomial of ω 2 with real coefficients, optionally with multiplicative iω factors. We emphasize that, for any lossless system with more than one port, the zeros of the S coefficients are different from the "S-matrix zeros", where det(S) = 0 and which always coincide with the conjugates of the poles [19].

D. Absorption and gain
In the presence of small absorption loss, S is no longer unitary, but it can be calculated perturbatively when all relevant modes have high Q. In particular, the denominators of Eq. (7) must obviously use the polesΩ of the actual absorptive system, however, if the QNM loss rates are split into radiative (Γ r ) and non-radiative (Γ nr ) parts, the numerators of Eq. (7) scale as DM −1 D ∼ Γ r (1 + const · Γ nr ) [no scattering when Γ r → 0]. We then see that, for high-Q modes (where both Γ r and Γ nr are much smaller than the frequencies), these numerators are equal to those of the lossless case to first order in Γ r , Γ nr , namely absorption only affects the poles. (A similar argument is often implicitly used in CMT, where it is typically assumed that Γ = D † D + Γ nr with D not changing in the presence of Γ nr because Γ itself is small [3], or is explicitly used to argue that the coupling coefficients to different CMT channels can be determined independently [29].) Therefore, in practice, to compute S for absorptive or active scatterers, we first calculate the QNMs (ω n , D pn ) of the lossless (radiative only) structure and use them to evaluate the numerators of Eq. (7). Then, we turn on absorption or gain mechanisms (adiabatically if needed for QNM-tracking purposes) to get the exact denominator poles (ω n ). Since the lossless-case D was finetuned for reciprocity, the non-unitary S will still be symmetric. The perturbation argument assumes high-Q modes but does not restrict the relative strength between radiation and absorption/gain rates. Indeed, as we see in the next examples, our QNMT gives quite accurate predictions even in the presence of modes with Γ nr Γ r and, in fact, even when relatively low-Q modes are present.

III. EXAMPLES IN ELECTROMAGNETISM
Our QNMT for the S matrix is applicable to all kinds of wave physics, such as acoustics [30][31][32], electromagnetics [3][4][5]7], and quantum mechanics [2,33]. Therefore, in our derivation, we used general physics-agnostic notation to render our results usable for any wave-scattering problem. In this section, to examine the accuracy of our QNMT, we study multiple examples in electromagnetism.

A. Normal incidence on microwave metasurface
We now study scattering of an electromagnetic plane wave with frequency f = ω/2π normally incident on the metasurface depicted in the inset of Fig. 2(a). It consists of alternating dielectric (green) and metallic (grey) layers, where the latter have been etched out to form two-dimensional square periodic lattices (of period a) of thin metallic crosses, whose centers are the same for all patterned layers. A square air-hole has also been etched throughout the entire thickness d of the metasurface in the region between the crosses. The metal thickness is 18 µm (corresponding to 0.5 oz copper). We study for frequencies below the first diffraction cutoff (f cut = c/a at normal incidence), so only transmission and reflection need to be considered. Moreover, there is C 2v symme- try, so the response for normal incidence is independent of the polarizationê and only 2 ports are needed. Numerical computation of the "exact" frequency response (S-matrix) for plane-wave excitation as well as of the eigenmodes for use in our QNMT is carried out using COMSOL Multiphysics [34]. Complex eigenfrequencies ω n are immediately obtained from the eigensolver, while the coupling coefficients D 1n , D 2n are computed from Eq. (3) as "power" inner-product surface integrals at the two (left/right for p = 1, 2) external boundaries of the metasurface between the QNM field ψ n = (E n , H n ) and the coupling port where, as emphasized earlier, only their ratio σ n is needed. In Appendix F, we provide further details and guidelines for the numerical simulations (mainly how to compute very-low-Q modes), and within the Supplemental Material [35]we give tables with the calculated QNMs (ω n , σ n ) for every structure presented in the paper.
For the asymmetric structure of Fig. 2, the parameters were chosen arbitrarily to test a very general response, with σ n departing substantially from ±1. Even so, we see a very good match between the exact numerical computation (black lines) and the QNM expansion of Eq. (7) (red lines), for both the amplitude of S 21 = |S 21 |e iφ21 (a) and its time delay , in both cases of lossless (solid lines) and lossy (dashed lines) structures. In the lossless case, we emphasize again that, due to our symmetrization procedure of Eq. (11), the QNM expansion we obtained is both unitary and symmetric. This is why the zeros of S 21 = S 12 are either real or complex conjugate pairs, and the zeros of S 11 , S 22 complex conjugates of each other, as they should [ Fig. 2(c)]. The response with copper and dielectric losses mostly maintains the same overall features, and merely exhibits reduced transmission at high frequencies and "superluminal" (0 ≤ τ 21 < 1) or negative group delay (τ 21 < 0) around transmission zeros. The latter does not violate causality [36], instead it has been shown to necessarily occur at peaks of absorption [37,38], and thus is indeed typically associated with lossy bandstop ("notch") transmission responses (such as zeros) [39]. Our QNMT correctly predicts even these unusual phenomena.
To quantify the benefits of our QNMT, we calculate the errors associated with not exactly enforcing reciprocity or energy conservation, for the same QNMs of the lossless structure. If the D coefficients are not finetuned with Eq. (11), S from Eq. (7) is not exactly symmetric, so Fig. 2(d) shows the resulting error in |S 21 − S 12 | 2 (orange curve). It is relatively small (although increasing at higher frequencies), indicating that Eq. (7) is already a good approximation. In contrast, the QNMT of Ref. 3 [analogous to our Eq. (9)] is reciprocal but violates energy conservation by large amounts, leading to non physical "absorption/gain". As indeed shown in Fig. 2(d), for this lossless 2-port, the sum of transmission and reflection |S pp | 2 +|S 21 | 2 (p = 1, 2) deviates from 1 by almost ±0.4 at some frequencies (blue curves)! It turns out that these large errors exhibit themselves mostly in the reflection coefficients. Additionally, in both cases, the violation of a physical constraint leads to errors also in the groupdelay prediction (such as negative group delay, which is impossible for a lossless 2-port). However, these delay errors are of less importance, since they usually appear around transmission zeros and they are mitigated when a finite-bandwidth pulse is considered [40], as detailed in Appendix G.
One key advantage of the QNMT method is that it resolves spectra around very-high-Q modes with perfect detail, while a frequency simulation requires a very dense uniform frequency grid to resolve them, being ignorant of their location. This, in turn, leads to a stark benefit in speed for QNMT. For this example, on the same machine and finite-element mesh, the QNMT calculation took an average of ∼ 60 secs per mode (×10 modes in Fig. 2), while the frequency-domain calculation an average of ∼ 100 secs per point (×600 points in Fig. 2).

B. 4-port metasurface via coupled polarizations
We now consider the 4-port system described in Fig. 3 which consists of a microwave metasurface with three dielectric layers sandwiching two metallic sheets, with patterned arrays of rotated cross-like apertures. The ports correspond to the two polarizations on the left (1,2) and right (3,4) sides of the structure. This system does not have the required symmetry for the normally incident plane-wave polarization to be conserved, instead the two orthogonal polarizations on each side cross-couple in both reflection and transmission, so a 4-port system is needed.
In Fig. 3(a), we plot the cross-polarization transmission |S 41 | 2 and again find a good agreement with our QNMT. In the Supplemental Material [35], we show the other components of the S-matrix. In Fig. 3(b), we again show for comparison the errors of QNMT without symmetry or unitarity. For the non-unitary QNMT of Ref. 3, the quantity 1 − p |S p1 | 2 reaches values as low as −0.5 (blue curve), in stark contradiction with energy conservation. The asymmetric QNMT using Eq. (7) without the D correction of Eq. (11) has errors in |S pq −S qp | 2 as large as 0.3 (orange curves). Moreover, in this example, it also has nonzero |S pq | 2 − |S qp | 2 with an error up to ∼ 0.015. This happens, because, for P -port systems with P > 2, unitarity alone does not guarantee |S pq | = |S qp | anymore, so our method of fine-tuning D to also enforce symmetry [Eq. (11)] corrects errors not just in the scattering phase, but in the amplitudes too. C. Oblique incidence on 2d photonic metasurface When a plane wave is incident on a metasurface at an angle θ, its transverse wave vector component at frequency ω is k ⊥ = ω sinθ/c. Phase matching then imposes that this must also be the Bloch wave vector within the metasurface. QNMT modeling with such excitation may seem intractable at first glance, if one tries to obtain a full band diagram to apply QNMT at each fixed real k ⊥ . However, we calculate here the relevant QNMs from a nonlinear eigenproblem, where the phase matching condition is imposed at the complex eigenfrequency by analytic continuation (thus giving a complex Bloch wave vector k ⊥ n = ω n sinθ/c). To find such unusual resonances, we developed software for two-dimensional (2d) dielectric structures, whose geometry can be split into uniform layered sections: at any complex ω, these complex Bloch modes are calculated within each section with a T -matrix formulation, then matched at interfaces between sections, and finally radiation conditions are applied to find the resonances (similar to CAMFR [41] and other interface mode-matching analyses [42,43]).
We then study scattering of a plane wave incident at a 30 • angle on a 2d photonic grating with its E-field transverse to the plane [inset of Fig. 4(a)]. For frequencies below the first diffraction cutoff f cut = c/a(1 + sinθ), the system again has only 2 ports. In Fig. 4(a,b), we show transmission, calculated both exactly (black curves) and with our QNMT (red curves). The agreement is indeed very good all throughout the range. The QNMTs without unitarity or symmetry reach errors ∼ 0.5 and ∼ 0.25 respectively [ Fig. 4(c)]. We highlight that QNMTs are expected to improve as more (higher-frequency) modes are included, however, the non-unitary formulation [3] exhibits large errors even down to very low frequencies.
Transmission is also plotted in the case of strong dielectric losses and it highlights that our perturbative approach works very well, even though now (due to absorption) several modes have decay ratesΓ n more than an order of magnitude larger than their rates Γ n for the lossless structure (see QNMs within Supplemental Material [35]).

IV. BACKGROUND SCATTERING REPRESENTATION BY LOW-Q MODES
In some design situations [7], one needs an effective slowly-varying background response, which has to be designed collectively and in conjunction with the high-Q modes. This background response (responsible also for the concept of Fano resonances [6]) is usually modeled in standard CMT via the direct-coupling matrix C in Eq. (2). In most cases, researchers have approximated C by simulating an effective background structure derived by some type of topology averaging which removes the high-Q resonances, with parameters chosen a posteriori for a best fit [4,8,9]. Here, we show how C can be calculated using the actual structure under study by appropriately combining its low-Q modes, providing also intuition for this dependency and its physical interpretation. Given the background C, we also derive how the total S can be computed and we test it for the electromagnetic examples of Sec. III.
A. QNMT with background C matrix Consider a lossless system supporting some high-Q resonant modes (ω n = Ω n − iΓ n , D n ), while the rest (ω C n = Ω C n − iΓ C n , D C n ) have much smaller Q (Fig. 1). Starting from Eq. (5), we combine these low-Q-mode terms within the sum to define C (ω) ≡ S {ω C n ,D C n } (ω), so that S takes the form of Eq. (4), with Ω including only the high-Q modes. In the limit Γ C → ∞, C becomes a frequency-independent matrix, which is unitary, since S(ω → ∞) = C. Therefore, we can use Eqs. (7) to calculate C and Eq. (11) to guarantee its symmetry. In this process, M C nl ≈ P p=1 D C pl D C * qn / Γ C l + Γ C n , so the Γ C n → ∞ cancel out between M C and the C-denominator poles, leading to a constant C not necessarily equal to −I. In the next subsection, we calculate this limiting C for some simple but useful pole configurations. Now, for the total S-matrix of Eq. (4), following the same procedure of Appendix B but with −I replaced by a unitary symmetric constant C, one can easily see that S † S = I is now equivalent to K t = −M −1 D † C. Therefore, the background C can be factored out of Eq. (4) to write: Here,S has the form of a separate scattering matrix, which itself also satisfies realness, unitarity and the properties of Eqs. (12,13). The condition of Eq. (10), to additionally impose exact symmetry on S, is modified to through which a finetuned D can be evaluated by an optimization procedure, as before, to obtain the final S. As we discuss in Appendix F, such extremely-low-Q modes are very difficult to locate numerically. Thankfully, in practice, real systems usually have modes with reasonably low Q, which are easier to find. However, the problem then is that their associated C (ω) is slowly varying instead of constant and is not necessarily unitary (as it does not describe a physical system by itself). Nevertheless, we can still use Eq. (7) for those low-Q modes to obtain a unitary (thus approximate) C, symmetrize it with Eq. (11), and then use Eq. (15) with the high-Q modes insideS to get S simply as a best-effort approximation to the actual system response. ThisSC construction still guarantees that S will also satisfy realness [S * (ω) =S * (ω)C * (ω) =S(−ω * )C(−ω * ) = S(−ω * )] and unitarity (S † S = C †S †S C = I), but does not guarantee reciprocity (S t = C tSt =SC = S) even though C t = C [at least, unitarity implies |S 12 | = |S 21 | for a 2-port]. Attempts to symmetrize S (for example, [44]) are not expected to have much success: if it were possible to exactly satisfy also reciprocity for a varying C(ω), one would then be able to build a unitary and symmetric S by multiplying unitary and symmetric C matrices formed by individual modes, which obviously is not possible. Regardless, as we show later in examples, for many physical systems, S =SC is a good enough approximation, which we use in separate work [7] to design accurate metasurface standard (e.g., elliptic) filters with a non-trivial background. In such design situations where a specific background is desired, C =S −1 S can alternatively be used to estimate and then design C(ω c ) at the target frequencies ω c without having to calculate any low-Q modes, rather by using S(ω c ) from a direct simulation andS(ω c ) from QNMT using only the high-Q modes.

B. C matrix due to Γ → ∞ modes in 2-port systems
In this subsection, we study some very basic configurations of high-Γ modes (ω C n , σ C n ) in 2-port systems to build intuition on their influence in shaping the background scattering. (For notational simplicity, here we drop the superscript " C ".) Single zero-frequency mode (−iΓ o , σ o )-Consider a 2-port system supporting a zero-frequency mode ω o = −iΓ o . As explained above, when Γ o → ∞, this mode can be factored out of S into a unitary symmetric background response C. The symmetry of C dictates that the modal ports-coupling ratio σ o is real. Using Eq. (D3) for just this mode, we find M = (1 + σ 2 o )/2Γ o and then which can give any "reflection matrix" in the orthogonal group O(2). A given transmission t is achieved for In particular, we obtain r = 1 for σ o = 0, while r = −1 for σ o → ∞. On the other hand, σ o = ±1 gives a fullytransmissive background with t = ±1. A concrete example which confirms this last result comes from considering a uniform material slab with thickness d and refractive indexñ. Its "Fabry-Perot" modes are an equispaced spectrum given by [13] ω n = 2c nd Consider now appropriately largeñ, so that atanh (1/ñ) π/2 ⇔ Γ o Ω 1 . Then, at frequencies of interest 0 < ω Ω 1 , the modal contribution (∼ Γ o /Ω n ) is negligible for all n = 0. In the limit d → 0, indicating the absence of slab, the only relevant n = 0 system mode has Γ o → ∞. Therefore, perhaps counter-intuitively, full transmission can be seen as equivalent to such a mode at 0 − i∞ with σ 0 = 1.
This result of t = 1 for a zero-thickness slab is a consequence our initial phase choice C = −I in Eq. (5). If we had instead chosen C = +I, we would have obtained t = −1 for zero thickness, which would be a valid but awkward phase convention.
When the structure has localized resonant elements but does not exhibit very strong overall reflection, the background response is commonly assumed to arise from the averaged geometry. When this is a simple uniform material slab, it corresponds to a Fabry-Perot system with the infinite set of equispaced alternating-symmetry modes given in Eq. (19) and a "comb" transmission response (see e.g., Fig. 5 blue circles and dashed lines). In the limitñ → 1, the system approaches a slab of free space of thickness d and its array of modes (equi-spaced by πc/d) is shifted down towards −i∞. The corresponding limiting value of C is simply the scattering matrix for propagation through d, namely |C 21 | → 1 and dφ C 21 /dω → d/c.
One-sided equispaced modes [nω o − iΓ o , 0 or ∞]-When instead the structure has strongly reflective components separating the two ports, the background response will comprise low-Q modes with fields mostly localized on either of the two port sides, thus having |σ n | 1 or |σ n | 1 respectively, adding up to |C 11 |, |C 22 | ≈ 1. If the average geometry on one side is a uniform material slab of thickness d/2 on a perfect mirror, its modes will again be Eq. (19) but for only odd or only even n. Asñ → 1 and Γ o → ∞, its round-trip reflection phase approaches dφ C 11 /dω → d/c.
From the last two examples, note that, although a conjugate-mode pair gives C ≈ −I in the large-Γ limit, when considering many such modes, their small pole contributions (deviations from −I), ∝ i2ω/Γ o from Eq. (20), add up to give a non-trivial phase term (in transmission or reflection).    which first touches the scatterer boundary [Eq. (3)]. A reasonable question is whether such a choice is always appropriate, especially in unusual geometries where a thin "needle" sticks out from the scatterer or when very-lowindex materials surround the scatterer. In such scenarios, the QNMs are expected to be localized close to the center of the scatterer and likely are already exponentially increasing inside the suggested outermost boundary. Here, we study a simple such photonic structure to explain the physics that come into play to render our boundary choice indeed suitable and we show how the S =SC formulation can give an interesting physical interpretation.
A plane wave at frequency f = ω/2π is normally incident on a uniform material slab of refractive indexñ s = 3 and thickness d that is attached to another slab of the same thickness but withñ = 1.05 (Fig. 5 inset). As a reference, the limiting symmetric system withñ = 1 and ports' cross sections A , B on the boundaries of theñ s slab has the response shown in Fig. 5 for (a) the poles (blue circles) with σ n = ±1, (b) transmission amplitude, and (c) phase delay (blue dashed lines). The response of the actual test system (shown with black "x" and lines) obviously approaches that of the limiting case in amplitude and has an additional phase delay due to the propagation through theñ slab A-A [Figs. 5 (b) and 5(c)]. It may be tempting to think that the two systems should have the same QNMT parameters (ω n , σ n ). However, its high-Q modes now have σ n different from ±1 (e.g., ω n d/2πc = 0.165−0.039i has σ n = −0.45+0.68i), exactly because they are calculated on our suggested boundary z p = A and the structure between A-B is asymmetric. How can we then expect QNMT to give an accurate prediction? The key lies in a set of very-low-Q modes ω C n that are supported mainly by the weakly scatteringñ slab [ Fig. 5(a) below the inset] and have highly asymmetric |σ C n | 1 [e.g., ω C n d/2πc = 0.24 − 0.23i has σ C n ∼ 10 −3 (1 + 1i)]. Asñ → 1, their Γ C n → ∞ and, when taken into account, the QNMT prediction of S gives the correct result with the additional expected group delay through A-A .
To clearly understand the effect of the low-Q modes, we separate them into a background response C, shown in green lines in Fig. 5. We see that C is almost a diagonal matrix [|C 11 | 2 ≈ 1 in Fig. 5(b)] and that C 11 represents a group delay equal to 2d/c [ Fig. 5(c)]. Indeed, this spectrum of low-Q-modes matches the "one-sided equi-spaced modes" of the previous subsection: due to n s n, A acts as a highly reflective boundary, so τ C 11 models the roundtrip phase propagation through A → A → A. The group delay ofS 21 through A-B is surprisingly reduced (and not increased!) by d/c compared to the A -B S 21 , and that is because the high-Q modes have modified σ n values. In this way, S 21 ≈S 21 C 11 works out correctly to give the required additional group delay d/c. It is worth reminding thatS does not have to be symmetric. In fact, oppositely toS 21 , the A-BS 12 has group delay increased by d/c compared to the A -B transmission, so that, combined with C 22 ≈ −1 = −e iω2d/c ≈ C 11 , we get the correctS 12 C 22 ≈ S 12 = S 21 ≈S 21 C 11 .
We conclude that all structural features contribute to scattering, which may be expressed via low-Q modes that need to be included in the QNM expansion. This justifies our choice of surface for the calculation of D in Eq. (3) to be the closest port cross-section that encloses the entire scatterer z p . For a different port choice z p , additional (practically impossible to locate) infinite-Γ modes would need to be accounted for to express the extra phase shift through z p − z p .

D. Examples of Section III
We finally test this formulation on the 2-port electromagnetic structures of Figs. 2 and 4 by using the modes with lowest Qs (much lower than the other modes) to construct the C matrix (see modes in Fig. 2(c) and the Supplemental Material [35]). C 21 indeed yields a slowlyvarying transmission background, as shown in Figs. 2(a) and 4(a) (green dash-dotted lines). Sharp features in S 21 are then obtained by the high-Q modes within the SC formulation (green solid lines), which actually gives a very good approximation of the overall transmission amplitude spectrum. Due to the remaining asymmetry of SC, we find that there are some errors in the group-delay prediction close to transmission zeros, but these errors are reduced when considering realistic finite-bandwidth pulses, as shown in Appendix G.
The 4-port structure of Fig. 3 has large connected metallic sheets, so it exhibits some very-low-Q "one-sided equi-spaced modes" (Sec. IV B) due to the outermost dielectric layers that only couple ports on the same side (namely they have σ r,pq → 0 or ∞ for p, q in opposite sides, so that |C pq | → 0). This explains the sharp features observed in S 41 , as only high-Q modes contribute to it. Moreover, the layers in this example are very thin, so these low-Q modes are located very far from the frequencies of interest anyway.

V. CONCLUSIONS
We have presented an expansion of the system scattering matrix S over non-normalized QNMs, formulated to satisfy the fundamental physical conditions of realness, energy conservation and reciprocity even for a small truncated number of terms. Resonant QNMs with frequencies ω n are computed with a numerical eigensolver. Coupling coefficients D pn are evaluated as surface overlap integrals between normalized-CPM and non-normalized-QNM fields (as only ratios σ r,pn = D pn /D rn are needed). Negative-frequency modes (−ω * n , D * pn ) are included. The D matrix is then adjusted through Eq. (11) and S is finally calculated via Eq. (7). For applications where it is convenient to separate an effective background response from the high-Q resonances, C can be determined by the same procedure using only the low-Q resonances, withS from the high-Q modes, and then S ≈SC. In Sec. IV B, we discussed several limiting cases, and showed that a nearly frequency-independent C with nonzero transmission can be produced by a very-low-Q mode on the imaginary-frequency axis, while a propagation phase is modeled by a set of several low-Q modes. The agreement of our QNMT with exact simulations gives us confidence that it can be successfully employed for rapid device design. In a separate paper [7], we indeed use this formulation to design precise standard (especially elliptic) high-order filters.
Our QNMT was mainly developed for linear ports with frequency-independent transverse mode profiles, such as plane waves. However, it could also be extended to finite arbitrary-shape scatterers using a spherical CPM basis. (Note that systems with spherical symmetry studied in previous QNMT formulations [3,5] can be modeled merely as multiple 1-ports, so QNMT was not really needed, as explained in the introduction.) Difficulties arise with other types of ports, such as when QNM-to-CPM coupling coefficients D(ω) are frequency dependent, and specifically when the S matrix has branch points due to CPM cutoff frequencies. While a rigorous extension of the theory to such systems may require a different approach, our model may still provide good approximate results for slowly varying coupling coefficients (such as for dielectric waveguides with low waveguide dispersion). The Weierstrass factorization theorem [45] states that a "meromorphic" function (analytic except for poles) can be factorized into a non-zero analytic function (an exponential) and a rational function (zeros and poles). If S(ω) is meromorphic, an exponential phase factor can then be factored out of each term as S pq (ω) = e iϕpq(ω) S pq (ω), with S pq a "proper" rational function (degree of numerator polynomial not larger than degree of denominator, so finite as ω → ∞) and ϕ pq an analytic function that we assume corresponds to a real phase shift for real frequencies. Similarly to the approach in Ref. 18, the combination of Phragmén-Lindelöf theorem (giving |e iϕpq(ω) | ≤ 1 in the upper-half complex plane) andČebotarev theorem shows that ϕ pq has to be a linear (= τ pq ω with τ pq ≥ 0, since the constant term can be added to S pq ).
When the system is lossless and reciprocal, unitarity and symmetry of S combine to S * (ω)S(ω) = I. The off-diagonal pq term can be expanded as r S * pr S rq = r e −i(τpr−τrq)ω S * pr S rq = 0. Since this has to be true for all real ω, all the phase terms in the sum must be the same, namely τ pr − τ rq = θ pq for each r. Specifically, for r = p and r = q, we get τ pp − τ pq = τ pq − τ qq ⇔ τ pq = (τ pp + τ qq ) /2. Therefore, the S matrix can be written as S = e iτ ω S e iτ ω with a diagonal τ matrix with real positive elements.
We compute the coefficient (p, q) of this matrix. The second and third terms are equal to: The final term, after decomposing into simple elements through and relabelling indices (n and l), becomes: These two equations show that, for a lossless system, S in Eq. (7) can be fully computed using only the resonant frequencies Ω and the ratios σ r of modal coupling among different ports, independently of the overall scaling factors in D. For 2-port systems, we simply choose r n = 1, skip the r subscript and denote σ n = D 2n /D 1n . The S matrix then becomes (p, q = 1, 2) When the system has mirror symmetry, these ratios can only take the values σ n = ±1, so they act like eigenvalues of the symmetry operator. Reversely, we also show that S uniquely determines Ω and σ r . In particular, for two different sets {Ω, σ r }, {Ω , σ r } such that S {Ω,σr} = S {Ω ,σ r } , we see from Eq. (D1) that Ω = Ω and Combining realness S * (ω) = S(−ω), unitarity S † (ω)S(ω) = I and symmetry S t (ω) = S(ω) on the realω axis, gives S(−ω)S(ω) = I, which can then be analytically continued in the entire complex-ω plane. Expanding this equation for a 2-port, we get We can use these equations to derive some useful properties of the zeros of S coefficients. As a reminder, realness requires all poles and zeros to be symmetric across the imaginary ω axis. Let us first assume that −ω o is , a behavior known as an "all-pass" phase filter. Moreover, when S(ω o ) is singular, ω o is called an "S-matrix zero" (for example generalized in Ref. 19).
The restrictions (i) on the S 21 zeros imply that its numerator is a polynomial of ω 2 with real coefficients, optionally with multiplicative iω factors.

Appendix F: Guidelines for numerical calculation of eigenmodes
Our QNMT relies on the calculation of the system resonances, including low-Q modes, which may be spread across the complex-ω plane and may include a mode with zero real frequency (ω o = 0 − iΓ o )-required for nonzero background transmission. Accurate numerical solution of the exact field equations for such eigenmodes can be challenging due to issues arising from spatial discretization and from truncation of the simulation domain in open systems, usually implemented with a "perfectly matched layer" (PML) [46,47]. In particular, a PML introduces the so-called "PML modes" [10,24], which can pollute the complex-ω plane and thus make it difficult to find the system's low-Q QNMs. Here, we suggest some simulation-parameter choices that seem to help calculate QNMs accurately, including the low-Q ones. For subwavelength structures (like the metallic microwave metasurfaces considered here), where quasi-static resonances can be sensitive to small spatial features, finite-element methods are attractive because they enable nonuniform meshing, so we used COMSOL Multiphysics [34].
Mesh type-In the computational domain (the cladding) around the scatterer (and in the PML when present), to minimize spurious reflections from changes in discretization one should use a mesh that is periodic (or nearly so) in the direction orthogonal to the boundary (where scattered waves should escape) [46,47]. In COMSOL, we employed a "swept" mesh for this reason, and found it to give smoother convergence of the complex eigenfrequencies with mesh density compared to a standard tetrahedral mesh and to work better for modes with zero real frequency. Modes computed numerically for a metasurface (of period a = 15mm) with cross-like apertures (of slits' width 0.02a and length 0.245a) on a metallic sheet sandwiched between two dielectric layers (of permittivity = 3 and thickness 0.5a). We only calculate the normally radiating modes of even symmetry, using cladding thickness tc, PML thickness tPML, mesh element size h in the cladding/PML, and: (a) PML backed by a perfect electric conductor ("x": tc = tPML = 0.5a, h = a/10; "o": tc = 0.25a; "+": tPML = a), (b) SBC without PML ("x": t clad = a, h = a/6; "o": t clad = 0.75a; "+": h = a/10), and (c) SBC with PML ("x": tc = tPML = 0.5a, h = a/6; "o": t clad = 0.25a; "+": h = a/10). The low-Q QNM at 0.27 − 0.1i is not computed accurately when using only a PML due to the nearby radiation modes, while those are pushed down in the complex plane when using a SBC.
Boundary condition, computational-cell size and mesh density-The commonly used frequency-independent PML of finite thickness rotates the continuous radiation spectrum in the complex plane and discretizes it into socalled "PML modes" [10,24]. Only QNMs above the rotated continuum can be accurately computed. The rotation is larger for a more absorptive (stronger or thicker) PML and for a thinner cladding, while the mesh density has little effect for fixed PML absorption. We illustrate these points for a metasurface example in Fig. 6(a): the exact QNMs are indicated with black dots, while the other symbols mark all the modes obtained numerically for different sets of computational parameters. Note that the low-Q QNM at ≈ 0.27 − 0.1i cannot be calculated accurately, unless one used an even larger rotation that could be computationally expensive. Moreover, it turns out that the field profiles of the QNM and nearby PML modes are quite similar, yet another obstacle to distinguishing the former. Even worse, for structures supporting a zero-frequency QNM (−iΓ o ), this mode cannot be found with a frequency-independent PML, as the negative imaginary axis will always lie below any rotated radiation spectrum [48].
Our QNMT is mostly applicable to CPMs with a frequency-independent transverse field profile. Thus, while computing QNMs, a "scattering boundary condition" (SBC) matching each CPM's profile can be used. As we show in Fig. 6(b), using an SBC but not a PML, the discrete radiation modes are now simply pushed down in the complex plane (instead of rotated as with the PML, although note that a similar push-down occurs with a frequency-dependent PML [48]). This clears the complex plane (up to the first diffraction branch point) and allows for an accurate computation of the low-Q and zerofrequency QNMs. The modal push is deeper for a thinner cladding and a denser mesh: these SBC-based discrete radiation modes disappear in the limit of infinite resolution, in contrast to PML modes [48].
Combining a PML and a SBC gives an even better result, as the discrete radiation modes are pushed down and then rotated in the complex plane [ Fig. 6(c)]. In summary, when possible, it is advisable to use SBCs matching the CPMs, along with a densely meshed PML and a thin cladding (which can practically be kept at a minimum for a PML mesh fine enough to also adequately model the QNM near fields).
Initial guesses for QNMs-Even if the discrete radiation modes have been pushed away, locating the low-Q (and zero-frequency) QNMs initially may require good initial complex-frequency guesses and searching for a large number of modes. In some convenient cases, the initial guesses can be provided by the analytically-solveable modes of an effective average structure. We discussed such cases in paragraphs "Equi-spaced (Fabry-Perot) modes" and "One-sided equi-spaced modes" in Sec. IV B.
Number of QNMs needed-Regarding how many modes one needs in QNMT for an accurate enough calculation, one observation is that a (low-Q) mode ω n = Ω n − iΓ n will have an effect at real frequency ω, if ω lies within a bandwidth proportional to the resonance linewidth, namely if |ω − Ω n | /Γ n < some number, which depends on the level of accuracy required. Obviously, when a mode has very large Γ n , it can play a role even if it is very distant and it may be very difficult to locate. In those cases, a background response C can still be estimated via C =S −1 S, with S from a direct simulation andS from QNMT using only the high-Q modes.
Appendix G: Group delay for a finite-width pulse For finite-bandwidth signal pulses, the group delay of a system should be computed as the difference in time expectation of flux between the input and output signals [40]. In general, for a signal f (t) with Fourier transform F (ω) = A F e iφ F (amplitude A F and phase φ F ), the flux time expectation is given by: where denotes a frequency derivative d/dω and the integral involving A F (ω) cancels due to realness [A F (−ω) = A F (ω) ⇒ A F (−ω) = −A F (ω)]. Group delay of S21 and S12 for oblique incidence on photonic grating presented in Fig. 4. The asymmetric S =SC approximation (green) has errors for (a) dφ/dω but agrees much better with the exact result (black) for (b) a more realistic Gaussian input pulse of spectral width σω = 0.01 × 2πc/a.
For a system with transfer function H(ω), the output signal g(t) has Fourier transform G(ω) = H(ω)F (ω), so φ G (ω) = φ H (ω)+φ F (ω). The system group delay is then equal to τ = t g − t f and, when φ F (ω) is constant (e.g., unchirped Gaussian pulse), it simplifies to: In the limit of a very narrowband (δ-function) input frequency spectrum F , this simplifies to the usual τ → φ H . In all QNMT formulations which do not simultaneously satisfy reciprocity and energy conservation (e.g., using Eq. (7) without Eq. (11), our S =SC formula, or Refs. [3,4]), the zeros of S coefficients do not abide by the required relations in the complex plane and this leads to errors in the phase (group-delay) response. To visualize them, in Fig. 7(a) we plot the (δ-function) group delays φ H of H = S 21 and H = S 12 for our approximate non-reciprocal formula S =SC for the oblique-incidence example of Fig. 4. We see that they mostly coincide (among them and with the exact result), except for frequency regions where the exact S 21 goes to zero, and in which τ 21 , τ 12 exhibit erroneous deviating spikes, often indicating "superluminal" (< 1) or negative delays. However, Eq. (G2) shows that frequencies corresponding to small amplitude H(ω) do not contribute much to the real time delay of a finite-bandwidth pulse. We therefore plot in Fig. 7(b) again the exact and bothSC delays for an input Gaussian pulse F (ω) ∝ exp[−(ω − ω ) 2 /4σ 2 ω ] with σ ω = 0.01×2πc/a. We see that the erroneous spikes have now disappeared and we get much better agreement, which improves as σ ω increases, as we confirmed.

II. QNM PARAMETERS FOR STRUCTURES
Here, we provide all the QNMs computed via finite-element simulations. The computed QNM-to-CPM ratios are indicated by σ c , while the ones finetuned by Eq. (11) of the main text are denoted by σ. The modes used to calculate the background C matrix are marked in bold. For simplicity, Ω and Γ in the tables denote the dimensionless frequency values as defined in the corresponding figures.