A coupling model for quasi-normal modes of photonic resonators

We develop a model for the coupling of quasi-normal modes in open photonic systems consisting of two resonators. By expressing the modes of the coupled system as a linear combination of the modes of the individual particles, we obtain a generalized eigenvalue problem involving small size dense matrices. We apply this technique to dielectric rod dimmer of rectangular cross section for Transverse Electric (TE) polarization in a two-dimensional (2D) setup. The results of our model show excellent agreement with full-wave finite element simulations. We provide a convergence analysis, and a simplified model with a few modes to study the influence of the relative position of the two resonators. This model provides interesting physical insights on the coupling scheme at stake in such systems and pave the way for systematic and efficient design and optimization of resonances in more complicated systems, for applications including sensing, antennae and spectral filtering.


Introduction
Metamaterials are a class of material engineered to produce properties that do not occur naturally [35]. Their fundamental building blocks, usually called metamolecules or meta-atoms, allow to mould the flow of light in an unprecedented way [20]. The study of these individual building blocks is of prime importance since their shape and electromagnetic properties governs the effective behaviour of the metamaterial. In particular, it is of great interest to distinguish which effects arise from the meta-atoms themselves or from their arrangement (periodic or random). In addition, there is an ever increasing demand for controlling optical fields at the nanoscale for applications ranging from medical diagnostics and sensing to optical devices and optoelectronic circuitry [33,25,19,16]. In particular, resonant local field enhancement is of paramount importance in phenomena such as surface enhanced Raman scattering (SERS) [34,26], improved non linear effects [5,8,13], optical antennae and the control of the local density of states [11,1] or spectral filtering [29,28]. Resonant phenomena are largely at the core of all those applications, calling for a systematic study of modal properties of complex nanostructures. Recently, a substantial amount of work have been devoted to the study of resonant interaction of light with nanophotonic systems, with particular focus on the computation and analysis of their quasi-normal modes (QNMs [24,23]), also known as leaky modes [22], resonant states [4,17], quasimodes [14,30] or quasi-guided modes [27] in the literature). These modes are an intrinsic properties of open resonators, and the associated eigenfrequencies are complex-valued (even for lossless materials), the real part giving the resonance frequency and the imaginary part being related to the linewidth of the corresponding resonance. A method based on quasi-normal modes expansion has been developed to compute the electromagnetic response of open periodic and non-periodic systems to an arbitrary excitation, allowing a simple and intuitive description of their resonant features [30]. A similar approach has been employed to define mode volumes and revisit the Purcell factor in dispersive nanophotonic resonators [23], and to compute QNMs of perturbed nanoparticles for sensing applications [31]. Another technique called the Resonant State Expansion (RSE [17,2]) consists in treating the system as a perturbation of a canonical problem which spectral elements are known in closed form. The idea is to compute these perturbed modes and to use them in the modal decomposition. However, those approaches are inherently perturbative and often rely on the fact that the considered perturbation must be much smaller or with low refractive index compared the unperturbed system. To study the interaction of electromagnetic systems, Coupled mode theories (CMT) [12] have been developed using various formulations and extensively applied in particular in the context of optical waveguides (see Ref. [9] and references therein). More recently, it has been extended to leaky eigenmodes with application in free space resonant scattering [6], absorption in semiconductors [32] and Fano resonances in photonic crystal slabs [3]. After expanding the modes of the two uncoupled structures onto the modes of isolated simpler systems, one obtains a system of ordinary differential equations, either in time or space. However, the coupled modes equations are of a rather heuristic nature and rely upon a slowly varying approximation, thus limiting the domain of application of the method. More recently, hybridization models have been introduced in the context of plasmonics [21], providing a simple and intuitive picture (an electromagnetic analogue of molecular orbital theory) to describe the plasmon response of complex nanostructures of arbitrary shape as the interaction of plasmons of simpler sub-structures. However, the formalism is limited to the quasi-static approximation. The method we develop here is general and relies only on the assumption that the coupled eigenmodes can be expressed as a superposition of the modes of the two uncoupled systems, which can be computed by a suitable numerical method or analytically in some simple cases. It can be applied to arbitrary shaped open resonators with both non trivial (possibly anisotropic) permittivity and permeability. We stress here that it handles quasinormal modes associated with complex eigenvalues. The imaginary part of the eigenfrequency giving the leakage rate of the mode is fully taken into account in our method. In contrast with the CMT where the coupled modes amplitudes are obtained through solving a system of ordinary differential equations, the proposed model reads a generalized eigenvalue problem for the complex coupling coefficients and eigenfrequencies. We illustrate the validity of the proposed model on a 2D problem of a high index dimmer of rectangular cross section rods for Transverse Electric (TE) polarization. The agreement between the model and full-wave finite element simulations is excellent and the convergence is analysed as function of the number of basis modes used. Using our model, we then study the coupling of the two nanorods lowest order modes for symmetric and asymmetric dimers by using the first two modes of the isolated resonator. This reduced order model depicts a simple hybridization picture, allowing a better understanding of the coupled modes in terms of their fundamental spectral building blocks.

Coupled mode model
Consider a system A (material properties µ a and ε a ) with eigenvalues Λ a n = (ω a n /c) 2 and non zero eigenvectors E a n and a system B (material properties µ b and ε b ) with eigenvalues Λ b n = (ω b n /c) 2 and eigenvectors E b n . We study the situation where the two resonators A and B form a coupled system denoted C (material properties µ c and ε c ), with eigenvalues Λ c n = (ω c n /c) 2 and eigenvectors E c n . The spectral equation for the three systems u ∈ {a , b , c} is: In the three cases, the eigenmodes are normalized according to [7]: where Ω denotes the computational domain. The only assumption we make here is to consider that the coupled eigenmodes can be written as a linear combination of the modes of the two isolated systems: Inserting Eq. (3) in the eigenvalue equation (1), taking into account the linearity of operator M µ c and projecting onto eigenmodes E a j and E b j , we obtain a generalized eigenvalue problem: Here we defined the eigenvectors C = [A i , B i ] T , i ∈ N containing the expansion coefficients and the block matrices N and M : where the coupling coefficients N uv ij and M uv ij , u ∈ {a , b}, v ∈ {a , b}, are defined as: Thus finding the coupled eigenmodes E c (i. e.the expansion coefficients A i and B i ) and the coupled eigenvalues Λ c boils down to solving the generalized eigenvalue problem (4). In practice we retain M a and M b modes of systems A and B respectively in the expansion (3), so that N and M are dense matrices of size M a + M b . On the other hand, the matrices involved in FEM based methods are sparse and of generally larger size. We continue by rewriting the permittivity and inverse permeability of the coupled system as ε c = ε a + ∆ε B = ε b + ∆ε A and µ −1 Here we introduced the permittivity and inverse permeability contrasts ∆ε U and ∆µ −1 U which are bounded by the resonator Ω U , U ∈ {A , B}: where ½ ΩA (resp. ½ ΩB ) is the indicator function of domain Ω A (resp. Ω B ). Using the linearity of the involved differential operators and the eigenvalue equations we obtain Exploiting the orthonormality of the modes, we have L aa ij = L bb ij = δ ij . We define the block matrices is the identity matrix of size M a (resp. M b ), and Λ a (resp. Λ b ) is a diagonal matrix containing the eigenvalues Λ a i (resp. Λ b i ). We can finally recast the eigenproblem (4) as Writing the eigenvalue problem under this form has two main advantages. First, it shows explicitly the relation between the spectrum of the coupled system and the spectra of the two uncoupled resonators through the matrix Λ. Secondly, the elements of matrices K and P are computed as overlap integrals over a smaller domain, either resonator A or B, thus the numerical calculations are faster and more accurate.
3 Numerical example

Symmetric dimer
To illustrate the method we consider a 2D example for the TE polarization case, i. e. E = E z e z . The structure is dimer consisting of two identical dielectric rods (ε = 16, µ = 1) of rectangular cross section (L x = 50 nm, L y = 3 L x = 150 nm) embedded in air and separated by a gap d = 1.4L x = 70 nm (cf. inset in Fig. 1 (a)). Results are plotted on Fig. 1 (a) for M = 200, where we can see that each eigenfrequency ω of an isolated rod produce a pair of coupled frequencies ω c ± corresponding to a bright and dark mode. Our model can predict the position of these coupled eigenfrequencies in the complex plane with very good accuracy. To quantify this we computed the relative error for the eigenfrequencies calculated with our model (denoted ω c n ) and the reference values for the dimer computed directly with the FEM (denoted ω c, ref We can see on Fig. 1 (b) that the relative error is between 10 −6 and 10 −10 for M = 200. Taking more basis modes into account improves the overall convergence, particularly for higher order modes, as can be seen from the relative error curves for M = 50, 100 and 200. This proves the convergence of the method and its ability to compute the spectrum of the coupled system with high accuracy. Note that here we included proper modes (i. e.quasi-normal modes associated with point spectrum) as well as improper modes (radiation modes associated with continuous spectrum, dicretized by using PMLs [18]). From a practical point of view, the computational time increases as M 2 with t c = 77s, 297s and 1233s for M = 50, 100 and 200 respectively.
To gain insight into the coupling mechanism between different eigenstates, we restrain the number of basis elements included in our model. We study the first four coupled modes (i. e.with the lowest real parts of their eigenfrequency) as a function of the gap of the dimer d. The traditional hybridization picture can be seen as a particular case of our model where one retains only one QNM in the expansion. Indeed, one mode usually dominates the coupling process for the low frequency modes, but this is not true for higher order modes. In the particular case studied here, we use only the two lowest order modes of the isolated rod : the magnetic Figure 2: Coupled modes of a symmetric dimer. Real (top) and imaginary (bottom) parts of the eigenfrequencies as a function of the normalized gap d/L x for the first four coupled modes: the magnetic dipole (MD, blue), the electric dipole along x (EDx, green), the electric dipole along y (EDy, orange) and the electric quadrupole (EQ, red). The coupled eigenfrequencies computed using our model with M = 1 (circles) and M = 4 (crosses) agree well with the coupled eigenfrequencies computed directly (solid lines). The first two modes of the isolated rod, the magnetic dipole (md) and the electric dipole along y (edy), are indicated with horizontal dashed lines. dipole (denoted md) and the electric dipole along y (denoted edy) (cf. Fig. 3). Symmetry considerations prove that these two modes do not couple to each other, so we can use a one mode model to find the coupled resonant frequencies of the first four modes of the dimer. In the case of non magnetic materials (µ = 1) we have K = 0 and because the system is symmetric with respect to the x = 0 plane, we have L ab = L ba = L, P aa = P bb = P ′ , P ab = P ba = P ′′ , and Λ a = Λ b = Λ. Under these conditions we obtain a simple analytical formula for the coupled eigenfrequencies: The associated eigenmodes are C + = (1, 1) which represent an in phase coupling (the so-called bright mode), and C − = (1, −1) accounting for an out of phase coupling (dark mode). In our case, the magnetic dipole (md) modes of the two isolated particles give rise to the magnetic dipole (MD) and the electric dipole along x (EDx) of the dimer thanks to a symmetric and antisymmetric coupling respectively (cf. Fig. 3). Similarly, the in-phase (resp. out-of phase) coupling of the two electric dipole modes along y creates the coupled electric dipole mode along y EDy (resp. the coupled electric quadrupole mode EQ). This single mode approximation is clearly valid for the particular setup used here as can be seen from Fig. 2, where the results from the one mode model (circles) are really close to the reference coupled eigenfrequencies computed directly (solid lines). However Eq. (7) is less accurate for the imaginary part of the coupled eigenvalues in the case of the coupling of two magnetic dipoles (cf. blue and green curves for MD and EDx). We thus included three modes of the continuous spectrum in the model (denoted cs1, cs2 and cs3 with eigenvalues ω cs1 = (0.0481−0.0726i) ω 0 , ω cs2 = (0.0593−0.1080i) ω 0 and ω cs3 = (0.0594 − 0.1090i) ω 0 ), which improved greatly the accuracy on the imaginary part (cf. crosses on Fig. 2). This result suggest that one need to include radiation modes in the model to account properly for the leakage of the coupled modes. We have checked that the accuracy does not improve if we add higher order QNMs for both real and imaginary parts. On the contrary, the accuracy increases when we add modes from the continuous spectrum. It is not clear from a mathematical point of view whether the QNMs form a complete basis outside the resonator in the general case as to the best of our knowledge it has been proved to be valid only inside the resonator and only for one-dimensional systems [15]. However, there is numerical evidence that QNMs and radiation modes are both needed in modal expansion for computing accurately electromagnetic fields [30]. From a physical point of view, our interpretation is that including modes from the continuous spectrum helps reconstructing the eigenfield outside of the resonator, describing more accurately its radiative nature and thus increasing the precision on the imaginary part of the eigenfrequency, which reflects the leakage rate of the mode.

Breaking the symmetry
Finally we study the case where one of the resonator is shifted of a distance d y along y, so that the symmetry of the original dimer is broken. The two fundamental modes of the dielectric rods (md and edy) can now couple to each other, so we include them both in our model. Results are reported on Fig. 4 (circles) and show good agreement with full wave simulations (full lines). When the shift increases, the coupled magnetic dipole resonance barely changes, EDx and EQ frequencies are redshifted while EDy is blueshifted. The prediction of the model on the imaginary part is correct but less precise. Adding the radiation modes cs1 and cs2 in the expansion reduces those discrepancies significantly (see crosses on Fig. 4 for M = 4).
The study of the coupling coefficients allow a quantitative evaluation of the coupling between modes induced by the breaking of symmetry. As the two resonators are identical, we have for a given coupled mode |A n | = |B n | (n being the index of the isolated mode, either md or edy), with a relative phase difference of 0 for the sym- metric coupling and π/2 for the antisymmetric coupling. We plot those coefficients as a function of the shift for the four coupled modes on Fig. 5. For coupled MD and EDx, the md still dominates the coupling, either symmetrically or antisymmetrically, which is illustrated by the fact that the norm of the related coupling coefficient |A| associated with md (red solid line) is greater than the one for edy. The situation is more complicated for the other two modes: for EDy, increasing the shift creates more symmetric coupling between the two isolated modes, resulting in a hybrid mode where the contribution of md and edy is almost the same for large shifts (cf. the bottom left panel on Fig. 5, where |A md | ≃ |A edx | for d y /L y = 0.5 and the associated field map on Fig. 4 (b)). The situation is similar for the EQ mode with an antisymmetric coupling of both nanorod modes, with an increasing contribution of the magnetic dipole as d y increases (cf. the bottom right panel on Fig. 5 and field map on Fig. 4 (b)).

Conclusion
We have developed a coupling model for computing the eigenmodes of two coupled electromagnetic resonators. Its validity is illustrated in the 2D TE case of high dielectric object through its good agreement in comparison with full wave simulations. The vectorial formulation is general and is suitable for arbitrarily shaped finite resonators with complex permittivity and permeability (possibly spatially varying and anisotropic). Future research directions will focus on applying the model to the Transverse Magnetic (TM) case and three dimensional problems, as well as extending it to periodic systems. Note that our approach can be straightforwardly extended to the coupling of an arbitrary number of objects by simply adding more terms corresponding to the modes of the different objects in Eq. (3), allowing the study of photonic oligomers [10]. This model provides physical insights on the coupling of individual modes and may ease the design and understanding of resonant processes required for many applications in photonics.

Acknowledgments
This work has been funded by the Engineering and Physical Sciences Research Council (EPSRC), UK under a Programme Grant (EP/I034548/1) "The Quest for Ultimate Electromagnetics using Spatial Transformations (QUEST)."