Quasinormal-mode expansion of the scattering matrix

It is well known that the quasinormal modes (or resonant states) of photonic structures can be associated with the poles of the scattering matrix of the system in the complex-frequency plane. In this work, the inverse problem, i.e., the reconstruction of the scattering matrix from the knowledge of the quasinormal modes, is addressed. We develop a general and scalable quasinormal-mode expansion of the scattering matrix, requiring only the complex eigenfrequencies and the far-field properties of the eigenmodes. The theory is validated by applying it to illustrative nanophotonic systems with multiple overlapping electromagnetic modes. The examples demonstrate that our theory provides an accurate first-principle prediction of the scattering properties, without the need for postulating ad-hoc nonresonant channels.

Scattering matrices have been playing a ubiquitous role in physics since the early history of quantum field theory [1]. Nowadays, scattering-matrix techniques represent an irreplaceable tool for scientists working in nuclear physics [2], electronic transport [3], or classically chaotic systems [4], just to mention some of the several fields of application. Scattering matrices also enjoy a well deserved popularity in electromagnetic modeling, ranging from microwave devices [5] to nanophotonics applications, such as scattering and transmission from nanostructured objects [6][7][8].
Most of the systems that are usually investigated with scattering-matrix techniques display a highly structured resonant response as a function of the excitation frequency (or energy), with the resonances in the spectrum being directly related to the poles of the analytical continuation of the scattering matrix in the complexfrequency plane [9,10]. For electromagnetic systems, such poles correspond to quasinormal modes (also called resonant states), i.e., complex-frequency solutions of Maxwell's equations with outgoing-wave boundary conditions [11][12][13][14][15]. In a sense, quasinormal modes represent the bare skeleton around which the frequency-dependent response of the system is built. The interplay among different electromagnetic modes has proven to be crucial for explaining several intriguing phenomena, such as Fano resonances in optical systems [16], scattering dark states [17,18], and the optical analog of electromagnetically induced transparency and superscattering [19], and for designing new optical materials, such as optical metasurfaces for wavefront shaping [20]. For these reasons, it is desirable and extremely interesting to be able to reconstruct ab initio the entire scattering matrix of a system from the knowledge of its quasinormal modes. Not only would such quasinormal-mode expansion contribute to the understanding of complicated spectral features in terms of interference and superposition of resonant states, but it would also offer practical advantages from the numerical point of view, since a full eigenmode calculation is generally faster and more comprehensive than a large number of single-frequency simulations.
A promising theoretical platform in which to carry out this program is represented by temporal coupled-mode theory for optical resonators. Such framework has been fruitfully employed to study the transmission of layered photonic-crystal structures [16,21,22], gratings [23], coupled cavities and waveguides [24,25], and the scattering cross section of nanoparticles [17,26,27]. For the moment, however, coupled-mode theory has been typically restricted to a selection of only one or two modes of the optical system. The residual spectral response is accounted for by a slowly varying frequency-dependent background, which is typically fitted from simulation data [16,21,22,24]. Part of the difficulty in expanding coupled-mode theory by including an arbitrary number of modes lies in estimating the coupling coefficients that relate the resonant states with the input-output channels. For a small number of modes, these can be obtained from symmetry considerations [16,24] or from the temporal decay rates [22]. However, in order to address the general case of multiple modes and an arbitrary configuration of input-output channels, a direct connection between the parameters of coupled-mode theory and the far-field properties of quasinormal modes is required.
In this work, we establish such a connection and we present a general theory to expand the scattering matrix on the quasinormal modes of photonic systems, which can be directly scaled to any number of eigenmodes and incoming or outgoing channels. The theory, based on the far-field asymptotic behavior of the modes and the unitarity property of the scattering matrix, represents a fully predictive tool that does not require the fitting of an additional nonresonant background. There are formal similarities between our results and the expansion of the electromagnetic Green function on normalized quasinormal modes, which is a well known result from classical electrodynamics [11,28,29]; of course, when the expansion of the Green function is known for any point in space, then the scattering properties of the system can also be obtained [29,30]. The theory that we present is formulated in a basis of input and output channels and it differs from these approaches in requiring only the far-field behavior of the modes at the input-output ports, as opposed to the full spatial distribution of the eigenfield. Moreover, our theory is independent of the choice of the normalization of the quasinormal modes.
Modal methods offer a deeper physical insight into the properties of resonant systems, because they allow us to draw a connection between the origin of complicated spectral features and the characteristics of the underlying quasinormal modes. For these reasons, they are particularly suitable for describing, understanding, and optimizing complex photonic systems. Notably, since the formalism that we present is derived on the basis of general coupled-mode theory, its range of applicability goes beyond that of classical electrodynamics.
The work is organized as follows. In Sec. I we derive the quasinormal-mode expansion of the scattering matrix, whereas in Sec. II we numerically validate the theory in the illustrative cases of photonic crystal slabs and multilayered metallic nanoparticles.

A. Quasinormal modes
In order to provide a rigorous motivation for the application of the formalism of coupled-mode theory to optical systems, we begin our analysis by establishing a direct connection with the theory of quasinormal modes. We consider a system of dielectric or absorbing photonic structures, described by a spatially inhomogeneous distribution of the dielectric function ε(r, ω). We assume that in the limit r → ∞, the dielectric function ε(r, ω) tends to the constant value ε b and we define ∆ε(r, ω) = ε(r, ω) − ε b [31]. The system supports a discrete number of quasinormal modes (also called resonant states) which are defined as the transverse complexfrequency solutions (Ẽ j ,H j ) of Maxwell's equations, with outgoing radiation boundary conditions [12][13][14][15]. This linear system of equation is equivalent to a quadratic eigenproblem for the electric field: As a consequence of the complex eigenfrequencyω j , quasinormal modes are characterized by a diverging amplitude in the far field. The same considerations also apply to systems that are periodic in one or two dimensions and radiating in the remaining dimensions. In this case, the one-or twodimensional crystalline momentum k is conserved and it is possible to define a family of quasinormal modes of the form:Ẽ k,j (r) = e ik·r E k,j (r), where the function E k,j (r) has the same periodicity of the system and the field satisfies the outgoing radiation boundary conditions along the nonperiodic dimensions.
To keep the notation general, we will assume the reciprocal wavevector to be fixed and omit the index k.
For one-dimensional dielectric media and threedimensional spheres, it has been proven that the modal eigenfields form a complete basis inside the structure, i.e., in the region where ∆ε(r, ω) = 0, provided that ∆ε(r, ω) or any order of its derivative is discontinuous at the boundary of its domain [11]. In this work, we make the assumption that the completeness hypothesis holds for arbitrary resonant systems, as well [15,28].
Following the usual scattering theory, we suppose that the system is illuminated by an incident field E b , which, in turn, is a solution of the wave equation (2) with only the background dielectric constant, ε b . Splitting the total field in the incident and scattered components, E(r) = E b (r) + E s (r), the latter can be shown to satisfy the inhomogenous wave equation in the presence of a source term proportional to the incident radiation, i.e., (4) Limiting ourselves to the region where ∆ε = 0, in the assumption that quasinormal modesẼ j form a complete basis, we can expand the scattered field on them: The exact expression for the coefficients a j depends on the incident field. Eventually, the knowledge of E s in a finite region is sufficient to extract the far-field properties of the scattered field, as it is described by the same Eq. (4), which becomes

This equation has the formal solution
withG b (r, r , ω) being the dyadic Green tensor of the background electromagnetic environment with homogeneous dielectric constant ε b . Since the integral in Eq. (6) is limited to the region where ∆ε = 0, we were able to replace the field expansion of Eq. (5). At this point, we expand the input field over a set of incoming waves (or, more generally, ports), E b (r) = α s +α E (+) α , and total electric field over a corresponding set of outgoing waves, E = α s −α E (−) α , whose detailed expression depends on the specific geometry of the system. Equation (6) clearly shows that the amplitude of each outgoing wave, s −α , can be written as the sum of a direct channel, which is directly proportional to the incoming amplitudes s +α , and a resonance-mediated channel, which is proportional to the quasinormal-mode amplitudes a j . In turn, the latter amplitudes are related to the incoming field through Eq. (4). From the linearity of Maxwell equations, it follows that all these relations can be written in terms of linear operators. This is the basis of the coupled-mode formalism, which we illustrate in the following.

B. Coupled-mode equations
Seeking a more general formulation, we write the characteristic equation of quasinormal modes, Eq. (1), as an eigenvalue problem for the effective "Hamiltonian" Ω+iΓ, Here and in the following, we assume the convention exp(iωt) for the temporal dependence of the field. The components of the vectors a j are interpreted as the coefficients of the expansion of the electric field in terms of quasinormal modes, according to Eq. (5). Due the inherently dissipative nature of quasinormal modes, the Hamiltonian operator Ω + iΓ is non-Hermitian and it has been split in the Hermitian and skew-Hermitian parts, which are expressed in terms of the two Hermitian operators Ω and Γ. Using the same language of Eq. (7) and following our previous considerations, we relate the incoming and outgoing amplitudes of the electromagnetic field (which we express in vector form as s + and s − ) by means of a system of coupled-mode equations: The operator C represents the direct-coupling channel, whereas the operators K and D account for coupling between quasinormal modes and the incoming and outgoing ports, respectively. Although there might be in principle infinitely many quasinormal modes and ports, for practical reasons we assume that the number of modes and ports is truncated to the finite values n and m, respectively. In this way, all the operators reduce to finite-size matrices. The set of Eqs. (8) and (9) is summarized by the scheme in Fig. 1. As originally demonstrated in Refs. [16] and [24], some relations among the quantities that appear in Eqs. (8) and (9) can be directly deduced from some very general physical properties of the system. First, electromagnetic reciprocity and energy conservation imply that and respectively. In Eq. (11), we have straightforwardly extended the theory to include the (Hermitian) decay matrix Γ nr , which accounts for absorption and other potential nonradiative-dissipation channels. Moreover, by comparing the dynamics described by Eqs. (8) and (9) with the time-reversed case and employing time-reversal symmetry, it can be shown that The system in Eqs. (8) and (9) has been extensively used to model the scattering properties of various photonic structures [17,21,22,24,27], proving itself particularly valuable for investigating the physical mechanism at the basis of various phenomena, such as the formation of Fano lineshapes in the spectrum as a consequence of the interference between the resonant and the direct-coupling channels [16]. In all these cases, however, the number of modes included in the equations is limited to one or two, and the direct-coupling channel, if present, is accounted for by fitting a specific frequency-dependent background response obtained from independent numerical simulations of the spectrum (see, for instance, Refs. [16,24,27]). The need for independent frequency-by-frequency simulations restricts the suitability of coupled-mode theory as a first-principle computational tool. Moreover, accurately fitting the direct-coupling background typically requires some additional assumptions which are difficult to interpret on physical grounds (for instance, the need for a frequency-dependent effective dielectric constant). In a broader sense, the actual separation between the resonant states and the frequency-dependent background is somewhat arbitrary, since the latter is also made up of a number of broad resonances associated with additional quasinormal modes. In the light of our assumption about the completeness of quasinormal modes (Sec. I A), we expect that by enlarging the set of electromagnetic modes, so to include the resonances usually associated with the background, we could remove the need for fitting the direct-coupling background and treat all resonant states on equal grounds. In this way, in addition to getting a more transparent physical picture, we could also better elucidate the modal structure at the basis of resonant systems. Implementing this strategy represents one of the main motivations for the formalism that we present in the next section, which is easily scalable to multiple modes with varying decay rates.

C. Expansion of the scattering matrix
The scattering matrix of the system connects the amplitude of the outgoing waves with the amplitude of the incoming waves: where we used identity (10). Here, we derive an expression for the expansion of the scattering matrix on quasinormal modes, on the basis of the system of Eqs. (8)- (9). To this purpose, in addition to the complex eigenfrequencies of the quasinormal modes,ω j (j = 1, . . . , n), we also assume the knowledge of the asymptotic behavior of the quasinormal-mode eigenfield in the output ports, which is equivalent to the knowledge of the relative complex amplitudes of the vectors For simplicity, we will refer to the vectors b j as the "scattering eigenvectors" of the system. As it is the case for all eigenproblems, the (complex) normalization constant of the eigenvectors can be set arbitrarily; however, as proven in App. B, the final expression for the scattering matrix does not depend on the choice of such constant. As a consequence, our approach is inherently normalization-free, at variance with other works dealing with the expansion of the dyadic Green function, which require the quasinormal modes to be normalized in a specific fashion [13][14][15].
In practice, the complex eigenvalues and the scattering eigenvectors need to be computed by numerical eigensolvers. The specific method depends on the definition of the input-output ports, but, in general, it involves calculating the electric field at a point or on a surface in the far-field region of the system, and, possibly, computing the projection integral of the field with the modal profile of the port. Some examples are provided in Sec. II. We stress that, since the scattering eigenvectors depend only on the far-field behavior of the resonant states, they can be obtained without computing the full distribution of the electromagnetic field over all space. This characteristic is particularly helpful, for instance, when quasinormal modes are calculated with numerical techniques such as the boundary-element method or the multipole expansion method, which typically benefit from a faster rate of convergence for far-field calculations.
Since the matrix Ω + iΓ is not Hermitian, the right eigenvectors alone are not orthogonal. However, as it is known from the theory of complex Hamiltonians [32], right eigenvectors (a j in our case) form a biorthogonal basis together with left eigenvectors, which are defined by the equation To simplify the notation, we introduce the n × n matrix A whose columns are the right eigenvectors a j and the corresponding matrix L of the column left eigenvectors l j . With this new notation, Eq. (7) becomes: withΩ being the diagonal matrix of the complex eigenvaluesω j . Moreover, we define the m×n matrix B whose columns are the vectors b j . The complex Hamiltonian of Eq. (7) can then be expanded on the biorthogonal basis as follows [32]: Even if the right eigenvectors are not orthogonal, they are however linearly independent [32]; thus, we can formally write L = (A † ) −1 . Replacing Eq. (17) into Eq. (13), we obtain the quasinormal-mode expansion of the scattering matrix, where we define Λ . = A T A and we use the relation B = DA, which comes directly from Eq. (14). For the moment, Eq. (18) represents only a formal result, which can be also seen as a special case of Mittag-Leffler's theorem on the pole-expansion of meromorphic functions [33]. For all practical purposes, it is crucial to derive an expression for the matrix Λ. This latter matrix plays a fundamental physical role, because the amplitude and phase of its terms determine the oscillator strength of each resonance and affect the degree of interference among the modes, which, in turn, has been found responsible for the appearance of interesting spectral features, such as Fano lineshapes [16] or the optical analogue of electromagnetically-induced transparency [19].
First of all, it can be shown that Λ is diagonal. This result follows from the symmetry of the complex Hamiltonian Ω + iΓ, which can be proven by combining Eqs. (11) and (12). The same result can also be derived from the requirement that the resulting scattering matrix must be symmetric [24]. Next, by multiplying each side of Eq. (12) by A * and after some algebraic manipulations, we obtain CB * = −BΛ −1 (A † A) * , which we can recast in the more compact form which defines the matrix Q = A † A. By multiplying Eq. (16) by A † on the left, taking the difference with its Hermitian conjugate, and employing Eqs. (11) and (14), we arrive at In general, the solution for Q cannot be written explicitly in terms of matrix products; however, it is straightforward to express it componentwise. First, in the case of no absorption (Γ nr = 0), we can write: This latter equation allows us to clarify the physical meaning of Eq. (19). With the aid of Eqs. (18) and (21), it can be shown that Eq. (19) is equivalent to the condition From the inversion of the scattering matrix, on the other hand, we obtain that S −1 (ω j )b j = 0, since quasinormal modes are defined as the self-sustaining solution of Maxwell's equations in the absence of any input radiation. Comparing the two results, it is clear that Eq. (19) guarantees that the scattering matrix is unitary at the modal eigenfrequencies, as required by energy conservation.
In the presence of absorption (Γ nr = 0), the energy balance must account also for the additional dissipation. In a broad sense, nonradiative processes represent a number of input-output channels that it is impractical to take into account directly. It is possible, however, to quantify their total effect on the decay rate of each quasinormal mode, for instance by calculating the shift of the imaginary part of the eigenfrequency with respect to the case when all losses are turned off. An example of this approach is discussed in Sec. IID. When the nonradiative decay rate is small compared to the frequency of the mode, the nonradiative term Γ nr can be treated as a first-order perturbation of the total Hamiltonian, i.e., we can assume A † Γ nr A A † AΓ nr =Γ nr A † A, whereΓ nr is the diagonal matrix of the first-order nonradiative decay rates, γ nr,j (j = 1, . . . , n). In this way, we can write the following generalized expression for Q: Equations (19) and (21) [or (23) for absorbing systems] allow us to fully determine the matrix Λ, and, hence, the quasinormal-mode expansion of Eq. (18). However, a closer inspection of Eq. (19) reveals that the system has m×n equations (the dimension of B) and only n unknows (the diagonal of Λ). Thus, for a given direct coupling matrix C, the system is generally overdetermined and a solution is not always guaranteed to exist. From a different perspective, the direct matrix C cannot be chosen freely, but it must satisfy some constraints that depend on the properties of the resonant states. In practice, it might be difficult to choose a direct matrix with a simple analytical form and, at the same time, consistent with Eq. (19), especially when a large number of quasinormal modes is involved.
For all these reasons, it is essential to develop a general theory that encompasses also the case when the matrix C is an approximation of the exact direct-coupling matrix. To this end, instead of looking for an exact solution of Eq. (19), we search for an approximate solution in the least-square sense. To be more precise, having defined the vectors x j (j = 1, . . . , n) as the columns of the matrix we look for the diagonal matrix Λ in Eq. (19) whose diagonal terms, λ j , minimize the objective function This reformulation of the problem does not affect the generality of the theory, because, if Eq. (19) has an exact solution, then such solution must coincide with the leastsquare one [34]. A simple calculation of the stationary points of the objective function leads to the result λ j = −x † j b j /(x † j x j ), which, once replaced into Eq. (18), provides us with the final expression Equation (26), together with Eqs. (24) and (21), is the desired expansion of the scattering matrix and it represents the main result of the present work. Using Eq. (24), the expansion coefficients can be explicitly written as The denominator of the coefficient can be regarded as a modified inner product that renormalizes the scattering eigenvectors in order to guarantee the total scattering matrix to be unitary. In the limiting case when the offdiagonal elements of Q are negligible, the expression in Eq. (26)

A. Photonic crystal slab
As an illustrative example, we consider a photonic crystal slab composed of a square lattice of circular holes etched in a silicon membrane (ε = 12.1). Indicating with a the lattice constant, we assume the slab thickness and the hole radius t = 0.4a and r = 0.2a, respectively. For normally incident light polarized along one of the lattice axes, we can limit ourselves to a single polarization of light; moreover, in the range of frequency ω < 2π/a, only the zeroth order of diffraction is available. As a consequence, the system can be effectively described with two ports, corresponding to the plane waves E 1,+ = s 1,+ e −iωz/c and E 2,+ = s 2,+ e iωz/c , propagating along the normal direction to the slab, which we indicate as the z axis.
In Fig. 2(a) we show the complex eigenfrequencies of the quasinormal modes of the system for normally incident light. Although all the modes represent equally valid solutions of the same characteristic equation (2) and they are treated on equal grounds in the expansion of the scattering matrix (26), it is useful from a physical point of view to distinguish between two categories of quasinormal modes: weakly dissipating quasi-guided modes and leaky modes with much larger radiation rates. As it appears from Fig. 2(a), the threshold between the two families can be set around Imω ≈ 10 −2 (2πc/a), with a difference of more than one order of magnitude between the corresponding imaginary parts of the eigenfrequencies. The leaky modes (Imω > 10 −2 2πc/a) have strong similarities with the Fabry-Pérot resonances of a homogeneous dielectric slab with an average refractive index n av , displaying a roughly constant frequency spacing of the order of the free spectral range δω = πc/(n av t). The deviation from the equal spacing behavior grows when the frequency increases, due to the wavelength becoming more sensitive to the dielectric-function inhomogeneity in the system [21].
Quasi-guided modes can be easily computed in various ways, including, e.g., frequency-domain [14] or timedomain [13,28] methods, or by determining the poles of the scattering or transmission coefficient in the complex frequency plane [9]. These techniques can also be combined, in order to exploit specific advantages. For instance, in the present example, the modes with Reω j > 0 have been computed by solving a linearized version of the eigenproblem in Eq. (2) with a commercial finiteelement package [36], whereas, for better numerical accuracy, leaky modes have been obtained separately by looking for the complex-frequency poles of the transmission amplitude computed with the Fourier modal method using a freely available solver [35]. Since the wave equation (2) is second order in the frequency, for each quasinormal mode with Reω j > 0 there exists a corresponding state withω j = −ω * j andẼ j (r) =Ẽ * j (r) [11], which has been included in the calculations, raising the total number of quasinormal modes under consideration in this example to n = 33. Due to numerical difficulties in performing the calculations near the imaginary axis of the complexfrequency plane, the decay rate of the Reω j = 0 mode has been estimated using the analytical formula for a homogeneous dielectric slab with an averaged refractive index [12].
The asymptotic behavior of the eigenfield is entirely determined by the inversion symmetry of the system with respect to the middle plane of the slab. Since the electric field amplitude is either even or odd with respect to the inversion, as indicated in Fig. 2(a), we can directly assume the scattering eigenvectors for even ("+") and odd ("−") modes. As we already remarked, since the scattering matrix expansion is indepen-dent of the eigenfield normalization, any other choice of the normalization in Eq. (28) would have been equally suitable. Finally, in agreement with our assumption about the completeness of quasinormal modes for photonic systems, we take the 2 × 2 identity matrix as the direct-coupling matrix In this way, we can derive the expression of the scattering matrix of the photonic crystal slab by applying Eq. (26) with the complex eigenfrequencies of Fig. 2(a) and the scattering eigenvectors of Eq. (28). The transmission intensity obtained from the resulting scattering matrix is shown by the solid curve in Fig. 2(b), and it is compared with an independent calculation by the Fourier modal method (dashed line) [35]. The agreement between the curves is excellent, highlighting the validity of the quasinormal-mode expansion of the scattering matrix. The comparison confirms that the first-principle description of the optical properties of the system provided by the theory is complete and accurate; moreover, we stress that such a description does not require any ad-hoc assumptions on the direct coupling channel and is based only on the complex eigenfrequencies of the quasinormal modes.

B. Asymmetric photonic crystal structure
A specific advantage of the scattering-matrix expansion is the straightforward applicability to generic systems lacking any particular symmetry. In order to illustrate this point, we consider a square lattice of L-shaped void structures partially patterned in a silicon slab. The shape and size of the structures is schematized in the inset of Fig. 3(a). The height of the patterned region (h = 0.2a) is one half of the total thickness of the slab (t = 0.4a), resulting in a configuration which is not symmetric by inversion along z. Moreover, for incident light polarized along one of the lattice main axes, the transmitted and reflected radiation will include a cross-polarized fraction. Therefore, we can model the system by defining four ports, corresponding to plane waves propagating above and below the slab and polarized along the two inplane crystal axis (which we indicate as x and y). In agreement with the assumption that quasinormal modes form a complete basis, we also assume the identity matrix as the direct-coupling matrix, i.e., C = I 4×4 .
The complex eigenfrequencies of the quasinormal modes, computed with the finite-element method [36], are presented in Fig. 3(a). Even in this case we can distinguish between a set of quasi-guided modes and a set of roughly equispaced leaky modes with a larger dissipation rate. Similarly to the previous example, the decay rate of the pair of modes with Reω j = 0 has been estimated using the analytical results for a homogeneous dielectric slab, and, moreover, we have also explicitly included the modes withω j = −ω * j . However, in this case the scattering eigenvectors b j must be obtained from the asymptotic behavior of the calculated quasinormal-mode eigenfield [37]. To this purpose, we consider the x and y electric-field components of each quasinormal mode in two planes located above and below the silicon slab at a sufficiently large distance to make the near-field contributions negligible. The specific choice of the distance does not affect the results, since only the relative amplitudes among the field components are relevant for the theory. It is interesting to note that the scattering eigenvector can also be computed with a near-to-far-field transformation of the quasinormal modes [38]. From the expansion of the scattering matrix in Eq. (26), we computed the total transmission intensity, T , and the cross-polarized transmission intensity, T xy (i.e., intensity of x-polarized transmitted light for ypolarized incident radiation). These quantities are shown (solid curves) in Figs. 3(b)-(c) and they are compared with the exact results (dashed curves) obtained from the Fourier modal method [35]. There is good agreement between the curves, especially in the vicinity of multiple narrow resonances, further confirming the validity of our approach as a predictive tool for computing the scattering matrix of electromagnetic systems. The small deviation from the exact result in the high-frequency region of Fig. 3(b) is likely due to the lower number of leaky modes included in this example with respect to the case of Sec. II A. The large radiative width of leaky modes (with a quality factor of the order of 10) implies that additional states beyond the frequency range under consideration may still have a small effect on the transmission in Fig. 3(b). To corroborate this hypothesis, we verified that the agreement with simulation data can be further improved when an additional pair of leaky modes at Reω j 1.1(2πc/a) is included in the scattering matrix expansion [37].
Computing the eigenvalues of Fig. 3(a) with the finiteelement method takes about one hour on a multiprocessor workstation. By comparison, on the same workstation the time required by a single frequency-point calculation of the transmission using the same finiteelement solver and the same mesh is about three minutes, implying that computing the transmission spectrum of Fig. 3(b) (roughly 1000 points) would require about 50 hours with the finite-element method. This 50fold reduction of computational time highlights the computational advantage of modal methods over frequencydomain full-wave simulations using the same electromagnetic solver.

C. Hybrid plasmonic system
Hybrid nanophotonic devices combining different photonic elements hold great promise for enhancing the functionality and the performance of various optical elements [39]. For instance, hybrid photonic-plasmonic systems made of plasmonic nanostructures coupled to optical resonators have been demonstrated to combine strong localization of light with precise control of the emission properties, enhancing the interaction with quantum emitters and the optical biosensing capabilities [40,41].
In order to exemplify the applicability of our theory to this class of systems, we consider the example of a square array of 60-nm-diameter metallic particles embedded in 200-nm-thick dielectric slab. For simplicity, we assume the metal dielectric function to follow the dissipationless Drude's model, ε(ω) = 1 − ω 2 p /ω 2 , with the plasma frequency ω p = 6eV. The complex eigenfrequencies of the quasinormal modes have been computed with the finiteelement method and are shown in Fig. 4(a). In addition to the Fabry-Pérot resonances of the slab, the calculation reveals the presence of a number of narrower modes. Such modes originate from the hybridization of the multipolar modes of the metallic particle due to the interaction with the polarizable dielectric. Similarly to the case of Sec. IIA, the modes can be classified as even or odd with respect to inversion symmetry along the direction perpendicular to the slab.
The transmission spectrum of the system is derived from the quasinormal-mode expansion of the scattering matrix [Eq. (26)] and it is shown by the solid curve in Fig. 4(b). Like the previous examples, we include the Fabry-Pérot zero-frequency mode with Reω j = 0 and we assume the unitary direct coupling matrix C = I 2×2 . The presence of a large number of modes results in a highly-structured spectrum with several closely-spaced minima and maxima of transmission, resulting from the reciprocal interference of light scattered by the polarization currents in the metal and the dielectric.
The transmission computed from the modal expansion of the scattering matrix is in very good agreement with the results of a frequency-by-frequency calculation with the finite-element method [dashed curve in Fig. 4]. The calculation of the complex eigenvalues in Fig. 4(a) takes a few hours on a multiprocessor workstation, comparing very favourably with the frequency-by-frequency computation, which requires about 30 hours using the same mesh. Furthermore, the modal expansion of the scattering matrix allows us to accurately resolve even the narrowest resonances, as demonstrated by the inset of Fig. 4(b), displaying a close-up of the spectrum in a small frequency range. This characteristic emphasizes an advantage of modal methods over the direct frequencyby-frequency computation, where a reduction of the frequency resolution over the whole extent would be highly impractical on grounds of the increased computational cost.
All these considerations can be directly extended to more realistic devices, such as metallic nanoparticles in interaction with large optical resonators and photonic cavities. In these cases, the systems are expected to benefit even further from the advantages of the modal expansion method, due to the increased size and complexity. Notably, the interest of determining the quasinormal modes is not limited to accessing the scattering properties of the system. For instance, it has been demonstrated that the quasinormal modes of an array of metallic particles can interact with molecular excitons, giving rise to plasmon-exciton-polaritons [42]. Thus, in addition to providing access to the scattering matrix, the modal information is also essential for describing and understanding the polaritonic effects.

D. Layered metallic particle
In order to highlight the generality of the theory, we consider a very different example. As demonstrated in Refs. [17,26,27], coupled-mode theory can be used to model the scattering and absorption cross-sections of spatially confined scatterers, such as metallic nanoparticles. Although in these works only one or two quasinormal modes are included in the application of the theory, our formalism allows the extension of the number of modes and channels in a straightforward way. Moreover, we also use this example to illustrate the application of the theory to absorbing materials.
For three-dimensional scattering objects, the ports correspond to incoming and outgoing spherical waves of degree l, order m, and both transverse electric (TE) and transverse magnetic (TM) polarization [6]. For simplicity, we consider a spherically symmetric system, where we can limit ourselves only to multipole terms with m = 1 and l > 0 [6,17]. The scattering and absorption cross section can be expressed as a function of the reflection coefficients, i.e., the diagonal terms of the scattering matrix, as follows [43]: (the index σ indicates polarization: σ = TE, TM). These expressions can be generalized to nonspherical scatterers by including the additional dependence on the order m of the modes [27]. For the sake of illustration, we consider a multilayered spherical particle with alternating layers of dielectric (ε = 2.1) and metallic materials, according to the structure sketched in the inset of Fig. 5. Core-shell metallic nanoparticles are a viable and well established platform for obtaining a significant local field enhancement together with a broad frequency tunability in the spectral response [44]. Here, we are mainly interested in the presence of multiple modes in each scattering channel, which underlines the advantages of our theoretical treatment in dealing with complex electromagnetic systems. Since we are considering a subwavelength particle sustaining plasmonic resonances, we limit ourselves to the lowest order TM-polarized modes (l = 1 and l = 2). We assume the metal dielectric function to follow Drude's model with plasma frequency ω p and nonradiative damping rate κ nr . In all calculations, we assume κ nr = 0.01ω p . This value is in agreement with those obtained from the fitting of the dielectric function of noble metals (e.g., gold) at frequencies lower than the onset of interband transitions [14]. The complex eigenfrequencies of the modes have been extracted from the position of the poles of the exact reflection coefficient in the complex-frequency plane [43] and they are presented in Fig. 5(b).
In the presence of absorbing materials, the theory requires the knowledge of the nonradiative decay rate of the modes, which is not directly available from our calculations, since the imaginary part of the complex eigenfrequency includes both the radiative and nonradiative components. In the case of low absorption, it is possible to distinguish the two contributions in an approximate way, by computing the complex eigenfrequencies twice, the second time upon setting Drude's damping rate to zero, and by taking the difference between the imaginary parts of the frequency in both calculations: The absorption cross section of the multilayered particles as calculated with our theory [Eqs. (26) and (23)] is depicted in Fig. 5(a) and it is compared with the exact result of generalized Mie theory [6,43]. The agreement of the curves is excellent, especially considering the additional level of approximation involved in estimating the nonradiative decay rates. Notable spectral features, such as the dip around ω = 0.33ω p , which is due to the interference between partially overlapping l = 1 modes, or the significantly different oscillator strengths of l = 1 and l = 2 modes, are well reproduced by the scattering matrix expansion. These results demonstrate that the theory can be easily extended to non-unitary systems, when an estimate of the radiative efficiency of each quasinormal mode is available [37].

III. DISCUSSION AND CONCLUSIONS
In this work, we derived a general approach to expand the scattering matrix of optical systems on the basis of quasinormal modes and we validated it with illustrative examples. The theory is directly scalable to any number of modes and input-output channels. This particular feature allows us to treat all resonant modes on equal grounds, going beyond the traditional partition of a system in a small set of narrow modes and a frequencydependent background fitted from simulation data. In this way, we achieve a more transparent picture of the modal structure of the system and, at the same time, we solve the ambiguity that could arise in defining the background channel in complex optical structures with a wide distribution of resonance widths. Eliminating the need for fitting a frequency-dependent background allows us to turn the qusinormal-mode expansion into a first-principle and self-consistent computational tool, which only requires the knowledge of the complex eigenfrequencies and the far-field behavior of the electromagnetic modes.
Creating artificial optical materials is an important goal in current nanophotonic research [45]. Such materials allow us to precisely control the intensity, phase, and polarization of scattered and transmitted light and to enhance light-matter interaction at the nanoscale. Spatial arrangements of optical resonators have been used, for instance, to realize high-contrast gratings [46], photonic metasurfaces [20,47,48], and zero-refractive-index metamaterials [49]. When the constituting optical resonators are chiral, several intriguing effects can be observed, such as the asymmetric transmission of circularly and linearly polarized light [50,51]. Even for a single optical resonator, like a multilayered particle, the interference of different resonant states give rise to interesting phenomena, such as, for instance, the optical analog of electromagnetically induced transparency and superscattering [19] and the formation of scattering dark states [17,18]. Multiple-resonance effects can also be exploited to tailor the scattering cross section of a scatterer, making it transparent to an outside observer [52]. Furthermore, hybrid photonic-plasmonic systems allow us to tailor the interaction with quantum emitters [40,41] and evidence polaritonic effects [42]. All these optical systems are typically characterized by a complex spectral structure, due to the presence of multiple electromagnetic modes coupled to the environment via various incoming and outgoing channels.
Our theory establishes a direct connection between the electromagnetic modes and the spectral properties of photonic resonant systems. The expression for the quasinormal-mode expansion that we derive is reminescent of the Breit-Wigner formula [2, 4], albeit with the some notable distinctions. A crucial difference is that the coefficient of each resonant term in the expansion depends on the frequencies and the amplitudes of all the other modes via a specifically introduced coupling matrix Q. This additional dependence reflects the fact that, whereas the application of the Breit-Wigner formula is restricted to non-overlapping resonances, no such limitation applies to the present theory, which accounts in a natural way for the effective interaction among different states originating from the coupling to a common external environment.
Typically, as an alternative to modal expansion, the scattering matrix and the derived quantities (such as transmission or scattering intensities) can also be computed with a full-wave solver on a frequency-by-frequency basis. The expansion on quasinormal modes, however, offers several advantages over direct frequency-domain computations on several aspects. In the first place, modal methods allow for a significant reduction of computational times [15,29], especially when the presence of narrow resonances dictates a very fine frequency resolution. The most computationally demanding phase of the modal expansion is the calculation of the quasinormal modes. After that, the method allows us to arbitrarily reduce the frequency resolution at no further computational cost.
More importantly, the scattering matrix expansion provides a more complete amount of information and offers a deeper physical insight with respect to a frequencyby-frequency calculation. This aspect is especially helpful, for instance, in the process of designing and optimizing optical materials. Building upon the connection between quasinormal modes and scattering properties established by the theory, instead of looking for a specific spectral feature among a large number of simulated spectra with varying parameters, one could equivalently search for a quasinormal mode with specific attributes. This strategy is generally faster, more transparent, and more suggestive of the relation among the physical parameters. For all these reasons, the quasinormal-mode expansion of the scattering matrix is particularly suitable for investigating the physical mechanisms at the heart of highly structured spectra, such as those arising from the interference of several closely spaced modes. Indeed, as we noted above, this is the case for many photonic systems which are currently the subject of intense research efforts. At the same time, the theory also represents a powerful and predictive tool for the first-principle calculation of the scattering behavior of general physical systems.

ACKNOWLEDGMENTS
This work is part of the research program of the Netherlands Organisation for Scientific Research (NWO). The authors acknowledge support from the European Research Council (ERC Advanced Grant 340438-CONSTANS) and from an industrial partnership program between Philips and NWO.

Appendix A: Case of orthogonal modes
If the scattering amplitudes of the quasinormal modes are orthogonal (i.e., b † i b j = 0 for i = j), or the spec-tral overlap between the modes can be neglected, the coupling matrix Q of Eq. (21) becomes diagonal. The least-square solutions of Eq. (19) can, then, be written as λ j = −b T j C † b j /(2Imω j ). As a result, the scatteringmatrix expansion of Eq. (26) assumes the simpler expression: This equation can be understood as a modified version of the Breit-Wigner formula [2, 4], in which the interaction between overlapping modes is neglected, but where the relation between the phase of each resonant term and the direct-coupling matrix C is retained.
Appendix B: Free choice of the normalization of the scattering amplitudes Here, we show that the result in Eq. (26) is independent of the normalization of the scattering amplitudes of the quasinormal modes. To this end, we consider a different set of amplitudes b j , which differ from the original b j by some complex multiplicative constants φ j (which can be different for different modes): Introducing the diagonal matrix Φ = diag(φ j ) and the matrix Q , defined by the expression in Eq. (21) with the modified eigenvectors, it is straightforward to verify that Q = Φ * −1 Q Φ −1 . In a similar fashion, we observe that Eq. (19) retains exactly the same form provided that Q, B, and Λ are replaced by Q , the column matrix of the new eigenvectors, B , and Λ = ΦΛΦ, respectively. Then, substituting these replacements in Eq. (18), we obtain that i.e., the expansion of the scattering maintains exactly the same formal expression independently of the choice of the eigenvector normalization constants.