Generalized Maxwell projections for multi-mode network Photonics

The design of optical resonant systems for controlling light at the nanoscale is an exciting field of research in nanophotonics. While describing the dynamics of few resonances is a relatively well understood problem, controlling the behavior of systems with many overlapping states is considerably more difficult. In this work, we use the theory of generalized operators to formulate an exact form of spatio-temporal coupled mode theory, which retains the simplicity of traditional coupled mode theory developed for optical waveguides. We developed a fast computational method that extracts all the characteristics of optical resonators, including the full density of states, the modes quality factors, and the mode resonances and linewidths, by employing a single first principle simulation. This approach can facilitate the analytical and numerical study of complex dynamics arising from the interactions of many overlapping resonances, defined in ensembles of resonators of any geometrical shape and in materials with arbitrary responses.

Dielectric optical nanoresonators are becoming an important platform for controlling light in nanoscale volumes of matter for many different applications [1][2][3][4][5][6][7][8][9][10] . The description of light-matter interactions in systems with two, or few, competing resonances is a relatively well understood subject [11][12][13][14][15][16][17][18][19][20][21][22][23] . Controlling systems with many overlapping resonances, conversely, is more challenging. A main difficulty lies in the fact that resonator modes are usually derived from the solution of Maxwell equations with radiating boundary conditions. These generate non orthogonal modal sets described by mathematical expressions that rapidly become difficult to manage, both analytically and in some cases also numerically, when the number of competing modes increases 24 .
In the field of optical waveguides, the study of multi-modal systems is a mature and developed area of research in both linear and nonlinear settings 25,26 . A significant contribution originates from the development in early days of exact theoretical frameworks that reduce Maxwell equations to a simplified set of equations of motion, which furnished the building block to understand complex hierarchical systems based on many interacting units 27-31 . In the field of optical resonators, an approximate form of this approach is available in time dependent coupled mode theory 32,33 , which is routinely used in many applications, such as the design of efficient broadband light energy trapping systems [34][35][36] , the study of nonlinear dynamics 37 , and the engineering of photonic crystals and metamaterials 38,39 . This theory derives dynamical equations obtained under the condition of a total energy of the system ε = ∑ | | a m m 2 expressed as the sum of independent terms | | a m 2 , each representing the energy of one resonant mode. The approximation originates from the lack of interacting contributions a a m n ⁎ , which are necessary to account for the presence of non-orthogonal modes in the electromagnetic field expansion.
In this article, we aim at unifying these two areas by studying an exact form of spatio-temporal couple mode theory (STCMT), which retains the simplicity of time dependent equations developed for photonic resonators, and the exact nature of coupled mode equations studied for multi-mode optical waveguides. The approach is inspired by the Feshback operator splitting designed to study the spectral statistics of open quantum systems 35,[40][41][42][43][44] , and here generalized to multiple projections "spaces" with the aid of different mathematics based on generalized functions 45 .
This formulation furnishes an intuitive description of Maxwell dynamics by providing an exact separation between propagating and resonant effects, within a simple set of exact equations that are particularly convenient for both analytical descriptions and numerical studies. We here illustrate fast numerical methods for calculating all the quantities of interest, ranging from the modes quality factors to the full density of states, from a single numerical simulation. This technique can also be generalized to describe wave dynamics in rigged-Hilbert spaces 46 .

Results exact spatio-temporal coupled mode theory via generalized Maxwell projections. The main
idea of this approach is to divide the space Ω into a set of adjoined regions (Fig. 1), and formulate the dynamics of light evolution independently in each set through the use of orthogonal eigenmodes.
In the space decomposition Ω = ∑ Ω n n , each region n Ω is composed into an interior spatial volume V n and a boundary surface S n , with a shape that is completely arbitrary. The union of all sets n Ω containing at least one optical resonator inside their volume defines the resonator space Ω = ∑ Ω r n rn , while the remaining volume of matter builds up the external space e r Ω = Ω − Ω . In each set , n Ω we formulate Maxwell's equations by resorting to the theory of generalized functions 45 and in particular by using the expression of the generalized differential operator ∇, defined as follows: S S n n δ ∇ = ∇ + ⋅ { } ∇ being the ordinary nabla operator evaluated at all the points inside the interior volume V n , n S n the unit vector normal to the surface S n , and δ δ = − x x ( ) S S n n a three-dimensional Dirac delta function centered on the surface ∈ S x n . By substituting the expression of ∇ from Eq. (1) into Maxwell equations, we obtain their generalized form: In Eq. (2), subscripts n and e indicate fields defined in the n-th resonator region Ω n and the external space e Ω , respectively, and δ δ = ∑ ⋅ − n n x x ( ) S S n S S n n n n the singular contribution terms arising on the surfaces S n separating the resonator space n Ω from the environment. We define each singular term x x ( ) (2) as arising from the limiting condition 0 ε → in which the surface S n is progressively approached from either the positive ( ) This condition leads to equations that are mathematically exact in the limit of 0 ε → , and well defined for every value of ε. The choice of which direction (±) to use to approach the surface S n is arbitrary and generates different boundary conditions for the equations defined in each region Ω n .
We here assume that the singular terms n E E ( ) are in the resonator space. The choice splits Eq. (2) into the following set, written for the n-th resonator Ω n and the external space e Ω , respectively: Figure 1. Example of space Ω = ∑Ω m partitioning with different sets m Ω characterized by elementary geometric structures. Panel b shows the equivalent dynamics inside Ω m , which contain a dielectric nanodisk resonator (orange area). Inside this space, e Ω is seen as an ideal PEC material and the dynamics are decomposed as a series of electromagnetic modes of the nanodisk terminated by PEC boundary conditions. www.nature.com/scientificreports www.nature.com/scientificreports/ represent the coupling of electromagnetic radiation at the surface of separation between the spaces n Ω and e Ω . The remaining singular terms, conversely, define an appropriate set of boundary conditions. In the case of zero electromagnetic field inside each resonator, E 0 n = and light dynamics in the external space reduce to: In order for this system to be mathematically well defined, we need to impose the absence of any singular term. This implies setting δ δ , which generates the following set of boundary conditions on S n in the limit ε → 0: S eS S eS n n n n Equation (5) show that that the whole resonator space Ω = ∑ Ω r n n is seen as a Perfect Electric Conductor (PEC) material from the external space. Analogously, by imposing the absence of any singular terms in the dynamics of the resonator space when the external field is absent, we obtain the set of boundary conditions for Ω e : Equation (6) imply that the external space is seen from within each resonator region Ω n as a Perfect Magnetic Conductor (PMC) material. Boundary conditions (5) and (6) lead to the following final set of Maxwell equations: ∆ into a linear contribution ε x E ( ) n n and a generic source term ∆ t J x ( , ) that keeps into account general types of effects, including dispersive effects, amplification and nonlinear responses. If we choose to approach the singular terms in Eq. (3) in a different way, we obtain different combinations of ideal PEC/PMC boundary conditions. The advantage of the splitting described by Eq. (7) is to decompose the dynamics of light into different spatial regions terminated by ideal PEC/PMC boundary conditions, which allow us to describe the evolution of the electromagnetic field with a complete eigenbasis of fully orthogonal modes. In the resonator space, orthogonal modes are obtained from the eigenvalue problem of Maxwell equations, written inside each space Ω n : nm nm m r nm nm 0 and terminated by PMC boundary conditions. The operator ∇ × { } is self-adjoint with PMC boundary conditions 47 . This implies that the resonator modes E nm , H nm are orthogonal, form a complete basis and possess a real frequency m ω . Mode orthogonality is calculated from the eigenvalue problem (8) using standard techniques 48 and occurs through the following relationship: V n nm n m r nm When = ′ m m , the integral expression in (9) represents the electromagnetic energy stored inside the resonator space Ω n by the m-th mode. Equation (9) is the counterpart of the orthogonality relation of guided modes in waveguides obtained via Poynting theorem (see, e.g., Eq. (2.2.52) of 27 ) and offers the same formulation advantages: when the electromagnetic field inside the resonator is expanded in terms of resonator modes (c.c. stands for complex conjugate): The time averaged electromagnetic energy t ( ) ε 〈 〉 dissipated inside the resonator space for monochromatic excitation at frequency ω becomes simply expressed as the sum of the energy of each mode: Analytic and closed form expressions of resonator modes and resonant frequencies in basic geometries are available from classical electrodynamics results of ideal metallic resonators 47 . As an example, for a single resonator space n Ω characterized by a generic cuboid volume with sides L L L , , x y z along the x y z ( , , ) axes and filled with a dielectric material with refractive index n r ( ), the frequencies of internal modes are: lpq l pq lpq x y z with l, p, q being integers. The corresponding magnetic field distributions are then expressed as follows: with A j being normalization constants. External modes, existing in the outer region Ω e , are then expanded as a series of ingoing and outgoing scattered waves. As the radiation spectrum is typically continuous, it is convenient to carry out the mode expansion in the frequency domain , PEC: www.nature.com/scientificreports www.nature.com/scientificreports/ E H , m m ± ± depend in general on ω through their wavevector k as, e.g., in the case of plane waves ± ⋅ e ik r , spherical waves ± e ikr r / or other types of traveling waves in the free-space. Following the same idea developed for the internal modes expansion in Eq. (11), we normalize radiating modes  ± E m and H m ∼ ± through an observable quantity of physical interest. We here use the optical power 49 , defined from the following integral when h h = ′: representing the union of all the surfaces of the resonator space r Ω . With the orthogonality condition (15), the effective power P flowing through S assumes the expression: S e e S n n n M 1˜˜) defining the vector of incoming ( ) + or outgoing − ( ) waves. The mode expansions carried out in Eqs. (11) and (14) reduce the time dynamics of Maxwell's equations to an exact set of spatio-temporal coupled mode equations, which relate the time evolution of the amplitudes of internal for a given set of impinging sources s t ( ) h+ . The mode expansions in Eqs. (10) and (14) express the corresponding spatial distribution of the field, providing a complete solution to the problem.
Coupled mode equations for the time varying coefficients can be found with two different approaches. One technique is to expand the electromagnetic field with Eqs. (10) and (14), substituting into Maxwell Eq. (7) and then projecting over each mode a m or ± s h by using the orthogonality relations (9) and (15). A second method is to exploit the linearity of Maxwell equations. We here employ a combination of both methods, starting from the latter.
In the external space Ω e , due to the absence of any source = with ˜ω ω D C ( ), ( ) being linear matrices. To write the equations describing the dynamics of t a( ), we first consider the case of linear materials with J 0 = ∆ . In this limit, Maxwell equations inside the resonator space Ω n are also linear, and the dynamics of t a( ) follow from the most general form of the linear time evolution equation for the modes a t ( ) m , with input sources corresponding to impinging waves s + : , and H t K t ( ), ( ) are linear matrices with Fourier pairs ω ωH K ( ), ( ) in the frequency space. The matrix H has no particular symmetry, and as such it can be decomposed as H iW = + Γ, with a Skew-Hermitian matrix iW and a Hermitian matrix Γ. Without loss of generality, we can assume that the matrix W is diagonal. If not, due to the Hermitian nature of W, we can always diagonalise W by a unitary matrix and project a m into a new orthogonal basis that will preserve the energy relation (11).
Matrices H K C D , , ,˜˜˜ are not independent, as the dynamics resulting from (17) and (18) have to satisfy energy conservation:˜ε By substituting the coupled mode Eqs. (17) and (18) into Eq. (19), we obtain the following self-consistency relations: www.nature.com/scientificreports www.nature.com/scientificreports/ with X 1/ being shorthand notation for the inverse matrix X 1 − . Equation (21) are similar to the time dependent coupled mode equations written in the frequency domain and originally introduced in 32,33 . However, there are also differences. In the traditional set 32,33 , all the linear matrices C, K, Γ are frequency independent and a t ( ) m are the amplitudes of traditional electromagnetic modes with radiating boundary conditions. Figure 2 shows a block diagram representation of Eq. (21). In the absence of any resonance, a 0  = and the system output is characterized by the open loop response +˜. This is the contribution that arises purely from propagation effects and away from any resonant light-matter interaction. When a 0 ≠  , the system response is characterized by a second term represented by the closed-loop feedback unit of Fig. 2, which forms the contribution of resonances.
Equation (21) and Fig. 2 show that the dynamics of Maxwell's equations depend only on three independent matrices: K , C, and W. The physical meaning of these matrices and their expressions is analyzed next. When K 0 = in Eq. (21), the n-th resonator n Ω and the external Ω e space are uncoupled, and the dynamics of (21) reduce to: It is always possible to obtain an input output representation of the dynamics where the equivalent scattering matrix is the identity matrix. This is accomplished by defining a new vector of outgoing scattered waves as follows: When the spaces Ω e and n Ω interact with nonzero couplings K , electromagnetic energy flows from the cavity region to Ω e , and viceversa. The coupling matrix K has in general a small, or weak dependence on the frequency ω. This condition, known as Markov approximation of open quantum systems 41,51 , is here discussed from the generalized Maxwell Eq. (7).
By substituting the field expansion (14) in (7) and by projecting over each traveling mode, we obtain the coupling coefficient elements: www.nature.com/scientificreports www.nature.com/scientificreports/ ( ) H m and E n are in phase. We discuss this condition with an illustrative example, derived in the case of a continuous external spectrum of plane waves e ik r ± ⋅ interacting with cuboid resonator structures Ω n with resonant wavevectors k lpq represented by Eq. (12). For a general frequency value ω = ck that is not resonant with any internal resonance k k lpq ≠ , the integral (24) is characterized by oscillatory terms of the type ( ) e i k k r lpq − ⋅ , which integrated through S do not furnish contribution. In all of these cases, the coupling matrix K is frequency independent.
In the situation where J 0 ≠ ∆ , the contribution of ∆ J to the dynamics is calculated by projecting Eq. (7) over the internal eigenmodes of the resonator space, thus obtaining an additional source term in the dynamics of the internal modes: being a vector of projected source terms with general contribution: L n with χ t x ( , ) being the non-instantaneous susceptibility, which can be further expressed in the frequency domain χ ω x ( , ) with a Lorentz multipole expansion 52 : where k ε ∆ is the difference between zero-frequency and infinity frequency relative permittivities, k ω is the frequency of the pole and δ k is the damping factor of each Lorentz oscillator. Equation (28) can express the response of any causal material, including losses effects via damping coefficients δ k .
For nonlinear effects, such as Kerr and other types of perturbative nonlinearities, the source term can be expanded as ε = ∆ ∂ ∂ J t P 0 NL , with 53 : is the m-th order non-instantaneous susceptibility tensor. The substitution of Eqs. (27)(28)(29) into Eq. (25) originate a system of dispersive nonlinear equations in the resonator modes a that can be used as the basis to study complex dynamics of nonlinear and/or dispersive light-matter interactions.
To study materials with gain effects, such as the ones originating in laser media, we can generalize the three dimensional quantum theory described in 54 : 4,9) are components of the atomic Bloch coherence vector = … S S S S [ , , , ] 1 2 9 , e is the electron charge, N a is the density of polarizable atoms describing the gain material, Q is the atomic length scale and x, y, and z are unit vectors along the x, y and z axes, respectively. The time evolution of the coherence S vector follows from Bloch theory: www.nature.com/scientificreports www.nature.com/scientificreports/ Quantities that can be calculated with this approach. Density of states. One of the most important quantities of a resonant system is the density of states (DOS), which is defined 55 as follows: n n and provides the number of eigenstates ω n in the frequency interval ω ω ω + d [ , ] . The most common technique for calculating the DOS of an optical structure characterized by complex geometries and dispersive effects is to extract it numerically via, e.g., finite difference time domain (FDTD) simulations, by injecting an impulsive point-dipole source S p in the electric E or magnetic H field such as: p l 0 with ê l ( = l x y z , , ) being a unit vector along one coordinate axis, x 0 the coordinates of a point inside the material whose DOS is to be computed, and δ = − p t t t ( ) ( ) 0 being a short time pulse with broadband spectrum. The local density of states (LDOS) measured at the point x 0 along the axis l x y z , , = is obtained from the power density spectrum of the electric or magnetic field measured at x 0 via the following relation 56 : The theory developed in the previous section allows for a fast calculation, which can furnish the complete DOS with just a single FDTD simulation.
By substituting the expression of the electromagnetic field inside each resonator region n Ω , given by Eq. (10), into Eq. (35) and by integrating over the volume V n defined by the resonator region, we obtain the DOS corresponding to the resonator region Ω n from the sum of the power density spectra of the internal modes: where the last step is obtained through the orthogonality relations (9). Equation (36) can be evaluated with a single FDTD computation, by injecting a broadband, three dimensional dipole source S t x ( , ) p centered at any point outside Ω n , and then measuring the time evolution of the internal modes a t ( ) m of the resonator space Ω n by projecting the electromagnetic field E n or H n over the corresponding eigenmode E nm or H nm via the orthogonality relations (9): The use of a radiative source in Eq. (33) to excite the spectrum of modes and compute the DOS via (34) and (35), or (36) implies that with this approach it is possible obtain modes that possess a finite quality factor. Dark modes, and other type of radiationless states that cannot be excited with radiating fields, are not measurable with this technique.
In the calculation of the DOS in the resonator space n Ω , it is also possible to use any orthogonal set of internal modes that result from the eigensolution of Eq. (8) with PEC/PMC boundaries.
This result is demonstrated from the orthogonality and completeness of the modes. Let Eq. (10) describe the electromagnetic field E n , H n inside a resonator structure Ω n , written in compact form as follows: we can then expand the same electromagnetic field by using a different set of eigenmodes pertaining to a material in Ω n with permittivity x ( ) www.nature.com/scientificreports www.nature.com/scientificreports/ Correspondingly, the density of states becomes: where the third equality stems from the completeness and orthogonality of the modes, as the reader can verify from (39). Equation (39) implies that the calculation of the DOS does not rely on the particular set of modes used, as long as they are computed with PEC/PMC boundary conditions. A particularly convenient choice in the case of dielectric structures are modes of completely filled cuboid resonator structures with constants ε ε ε = x ( ) r (0) 0 and (0) 0 μ μ = , which are analytically expressed by Eqs. (12) and (13) via simple trigonometric formulas. Other possible simple choices are represented by spherical or cylindrical n Ω spaces characterized by analytic combination of Bessel functions. When using these equivalent modes expansions, the resonance matrix W appearing in Eq. (22) is in general not diagonal, due to Eq. (39) that project W into a matrix of full rank. Figure 3 summarizes this procedure with an example of DOS calculation. We consider a resonator schematically illustrated in Fig. 3a (orange area). We partition the space by using a cuboid resonator region 1 Ω (Fig. 3a  green area). We illuminate the structure by a single broadband pulse source (Fig. 3 input pulse)  Photonic resonance networks. A particularly important quantity in the analysis of resonant systems is the mode quality factor = ω τ Q 2 0 , defined as the product between the mode frequency ω and the mode decay rate τ, and describing the ability of the mode to trap and release electromagnetic energy 32 . Traditionally, the evaluation of the Q factor requires to selectively excite each mode and compute ω and τ from the mode decaying rates in time, or from the mode frequency linewidth ∆ in the LDOS. This approach assumes non-interacting resonances and cannot be directly applied in the general case of overlapping resonant states in the spectrum.
The theory developed via the generalized Maxwell's equations allows to precisely evaluate the Q factor of all the resonances of the system from a single simulation in the general case of interacting modes. To illustrate the calculation, we begin by solving Eq. (18) in the Markov limit where the input frequency ω is away from a resonant frequency of the system: expressed as the sum of complex damped exponentials with mm α ′ constant coefficients arising from the matrix product of The Fourier transform of the mode a ( ) m ω is a complex rational function: By equating Eqs. (45) and (46) the quality factor Q m associated to the resonant mode at m ω is: with  and  being the real and imaginary parts of s m , respectively. The calculation of the network of mode quantities Q m , ω m and τ m can be accomplished via a single FDTD simulation, by first calculating the DOS following the procedure outlined in the previous section and by then extracting the poles via rational fitting through Eq. (45). For this task, we used the stable pole extraction algorithm recently developed and detailed in 58 , which is mathematically exact for rational models and can automatically detect the order of the rational polynomial in the DOS from its singular matrix. Figure 4 illustrates the accuracy of this technique in the example case of = .. m 1, , 7 overlapping resonances s m , with random frequencies and damping factors contained in a narrow band and generating a single apparent resonance line in the DOS (Fig. 4a). The solid markers in Fig. 4b show the position of the resonances and damping in a two dimensional ω γ ( , ) space, with the area of each marker being proportional to the Q factor of each mode. Figure 4c presents the results of the iterative algorithm for automatic detection of the polynomial order in (45). The efficiency of the algorithm increases exponentially and after a few iterations the system can correctly detect all the resonances (Fig. 4b, cross markers) composing the DOS (Fig. 4a, solid line), with differences between 10 10 − and − 10 25 (Fig. 4d). examples of applications. One dimensional structures. We begin by considering one dimensional structures, which illustrate the application of the theory via fully analytic calculations. Figure 5 shows the structure   = and ˜= =− C C 1 11 22 arise from PEC boundary condition at the resonator space n Ω , originating reflections for each incoming source when light is injected from the external space e Ω . In the geometry of Fig. 5a, the scattering matrix is already in diagonal form.
By exploiting the self consistency relations (20), we can express the diagonal elements of the damping matrix Γ from the coupling coefficients = | | θ K K e ij ij i ij : the symmetry of the resonator structure along z implies that damping factors along channels 1 and 2 are the same: . The remaining elements of the damping matrix are then: in Eq. (51) two cases are possible due to the z symmetry of the system, and each internal mode can either decay symmetrically or anti-symmetrically in the scattering channels:  (54) is sufficient to calculate all the elements of the damping matrix via Eqs. (50) and (53). Equation (54) is a direct consequence of the phase matching condition discussed for the coupling matrix coefficients expressed by Eq. (24).
We verified these results, and in particular the validity of (54), by FDTD simulations. We begin by calculating the coefficients a t ( ) n by projecting over the magnetic field eigenmodes defined by (48). To appropriately resolve narrow peaks with discretization of minimum 30 points we need at least 0.002       ω π d c 2 frequency resolution. To achieve this requirement we launch simulations for 50000 timesteps, dt = 5.87 ps, and corresponding frequency . Figure 5(b,c) shows the spectral distribution ã ( ) m 2 ω of the first seven modes in the expansion. The modes are grouped into even (Fig. 5b) and odd (Fig. 5c). In agreement with the analytic results based on (54), the amplitude of each mode vanishes at the internal resonant frequency l m ω ≠ of the other modes with the same symmetry along z.
From the amplitude intensity at the mode internal frequency ω a ( ) l 2 , we apply (54) and calculate all the elements of the coupling K and damping Γ matrix. Figure 5d illustrates the transmissivity and reflectivity of the structure calculated from the analytic solution via multilayer theory (dashed line) and from the network model based on Eq. (21). The solutions are exactly the same, with relative differences below − 10 16 . In the analysis we did not use any fitting curve, background or parameters, but calculated all coefficients from the analysis of internal modes by using Eq. (54). Figure 6 shows the representation of the network mode parameters extracted from the poles of the DOS (a, circle markers) of the resonator of Fig. 5a. The DOS is calculated from Eq. (40) by summing up the spectral mode densities a ( ) m ω | | is illustrated in Fig. 5b,c. The rational model representing the DOS via Eq. (45) is reported as a solid line in Fig. 6a. The network representation illustrated in Fig. 6b correctly predicts a mode network composed by modes with identical damping factors γ m and resonant frequencies following the exact analytic formula (red dashed line) of a cuboid resonator terminated by PMC boundary conditions. The corresponding quality factor of the modes (Fig. 6b, solid area of each marker) increases in ω due to the increasing resonant frequency ω m of each mode.
Two and three dimensional structures. We apply the STCMT approach to describe more complicated geometries for both non-periodic (Fig. 7) and periodic (Fig. 8) boundary conditions. Figure 7a illustrates the schematic representation of a two dimensional resonator ring shaped geometry with 0.15 μm and 0.5 μm sizes for the inner and outer diameters, respectively. We define a single cuboid resonator space Ω 1 (Fig. 7a green area) that includes the whole resonator, and use the orthogonal mode set of the cuboid volume without the resonator space, as employed www.nature.com/scientificreports www.nature.com/scientificreports/ in Fig. 3. The corresponding calculated DOS is fitted with the rational expression (45) on Fig. 7b. Figure 7c provides a zoom of the results on a target frequency range covering most of a visible spectra. The red dotted lines represent the boundaries of the target frequency range. The resonance network with the various Q m factors of the resonator modes is presented in Fig. 7c.  www.nature.com/scientificreports www.nature.com/scientificreports/ Figure 8a-d shows equivalent results obtained from a periodic two dimensional resonator characterized by a complex geometry consisting of the concentric superposition of a cross shaped geometry and a disk. In this case, due to stronger mode coupling via periodic neighbor resonators, the peaks in DOS spectra are wider and the overall DOS is smoother.
As discussed in the theory, the STCMT allows to obtain the full dynamics of the field, including the spatial distribution of any resonant mode existing inside the resonator. We illustrate this approach with reference both the non periodic ring resonator of Fig. 7a and the periodic resonator of Fig. 8a. We consider two resonant peaks in the DOS, labeled (c) and (d) in Fig. 9a for the non periodic resonator and Fig. 10a for the periodic one, and extract the modes spatial profile by Eq. (10) using the orthogonal eigenfunctions Eq. 12 (Fig. 9b) of the cuboid resonator space Ω 1 used in the projections of the time varying coefficients a t ( ) m . The corresponding energy distributions of the resonant modes are portrayed in Figs. 9c,d and 10b,c. They represent whispering-gallery modes of the ring-like geometry in the nonperiodic case, and more complex Floquet-Bloch modes in the periodic structure. Figure 11a-d shows the calculation results for a three dimensional resonator structure, delimited by a cubical Ω 1 resonator space. The resonator here is a compound shape consisting of a concentric superposition of a sphere and a disk (Fig. 11a). The mode network of this structure (Fig. 11d) is mainly represented by two close resonances with quality factors Q 52 1 = and Q 21 2 = .
conclusion In this work we formulate an exact spatio-temporal coupled mode theory for arbitrary resonator structures, derived with orthogonal and complete eigenmodes obtained from Maxwell equations with generalized operators. Using this theory, it is possible to provide an exact representation of the electromagnetic dynamics in both space, time and frequency via a simple set of exact equations of motion, in which all relevant quantities such as the DOS, the modes quality factors, and the modes spatial distribution can be calculated numerically from a single first principles simulation. We provide examples of this approach in one, two and three dimensional optical structures. We believe that this approach can help the design of photonics systems based on complex multi mode interactions, providing an exact formulation of coupled mode equations in general conditions of overlapping resonances and for arbitrarily defined materials and resonator geometries.