Abstract
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 spatiotemporal 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.
Introduction
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 lightmatter 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 multimodal 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,28,29,30,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 \(\varepsilon ={\sum }_{m}{a}_{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}_{m}{a}_{n}^{\ast }\), which are necessary to account for the presence of nonorthogonal modes in the electromagnetic field expansion.
In this article, we aim at unifying these two areas by studying an exact form of spatiotemporal 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 multimode 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 riggedHilbert spaces^{46}.
Results
Exact spatiotemporal coupled mode theory via generalized Maxwell projections
The main idea of this approach is to divide the space \(\Omega \) 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 \(\Omega ={\sum }_{n}{\Omega }_{n}\), each region \({\Omega }_{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 \({\Omega }_{n}\) containing at least one optical resonator inside their volume defines the resonator space \({\Omega }_{r}={\sum }_{n}{\Omega }_{rn}\), while the remaining volume of matter builds up the external space \({\Omega }_{e}=\Omega {\Omega }_{r}\).
In each set \({\Omega }_{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 \(\nabla \), defined as follows:
\(\{\nabla \}\) being the ordinary nabla operator evaluated at all the points inside the interior volume \({V}_{n}\), \({{\bf{n}}}_{{S}_{n}}\) the unit vector normal to the surface \({S}_{n}\), and \({\delta }_{{S}_{n}}=\delta ({\bf{x}}{{\bf{x}}}_{{S}_{n}})\) a threedimensional Dirac delta function centered on the surface \({\bf{x}}\in {S}_{n}\). By substituting the expression of \(\nabla \) from Eq. (1) into Maxwell equations, we obtain their generalized form:
In Eq. (2), subscripts n and e indicate fields defined in the nth resonator region \({\Omega }_{n}\) and the external space \({\Omega }_{e}\), respectively, and \({\delta }_{{S}_{n}}{{\bf{n}}}_{{S}_{n}}={\sum }_{n}{{\bf{n}}}_{{S}_{n}}\cdot \delta ({\bf{x}}{{\bf{x}}}_{{S}_{n}})\) the singular contribution terms arising on the surfaces \({S}_{n}\) separating the resonator space \({\varOmega }_{n}\) from the environment.
We define each singular term \(\delta ({\bf{x}}{{\bf{x}}}_{{S}_{n}})\) in Eq. (2) as arising from the limiting condition \(\varepsilon \to 0\) in which the surface \({S}_{n}\) is progressively approached from either the positive \((+)\) or negative \(()\) side, indicated as \(\delta ({\bf{x}}{{\bf{x}}}_{{S}_{n}}\pm \varepsilon )\). This condition leads to equations that are mathematically exact in the limit of \(\varepsilon \to 0\), and well defined for every value of \(\varepsilon \). 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 \({\varOmega }_{n}\).
We here assume that the singular terms \({\delta }_{{S}_{n+}}{{\bf{n}}}_{{S}_{n}}\times ({{\bf{E}}}_{n}{{\bf{E}}}_{e})\) and \({\delta }_{{S}_{n+}}{{\bf{n}}}_{{S}_{n}}\cdot ({{\bf{B}}}_{n}{{\bf{B}}}_{e})\) are approached from the external space, with \({\delta }_{{S}_{n+}}=\delta ({\bf{x}}{{\bf{x}}}_{S}+\varepsilon )\), while the remaining singularities \({\delta }_{{S}_{n}}{{\bf{n}}}_{{S}_{n}}\times ({{\bf{H}}}_{n}{{\bf{H}}}_{e})\) and \({\delta }_{{S}_{n}}{{\bf{n}}}_{{S}_{n}}\cdot ({{\bf{D}}}_{n}{{\bf{D}}}_{e})\) are approached from within the resonator space, with \({\delta }_{{S}_{n}}=\delta ({\bf{x}}{{\bf{x}}}_{S}\varepsilon )\). This choice implies that, at every value of \(\varepsilon \), the terms \({\delta }_{{S}_{n+}}{{\bf{n}}}_{{S}_{n}}\times ({{\bf{E}}}_{n}{{\bf{E}}}_{e})\) and \({\delta }_{{S}_{n+}}{{\bf{n}}}_{{S}_{n}}\cdot ({{\bf{B}}}_{n}{{\bf{B}}}_{e})\) are contained in the environment, while \({\delta }_{{S}_{n}}{{\bf{n}}}_{{S}_{n}}\times ({{\bf{H}}}_{n}{{\bf{H}}}_{e})\) and \({\delta }_{{S}_{n}}{{\bf{n}}}_{{S}_{n}}\cdot ({{\bf{D}}}_{n}{{\bf{D}}}_{e})\) are in the resonator space.
The choice splits Eq. (2) into the following set, written for the nth resonator \({\Omega }_{n}\) and the external space \({\Omega }_{e}\), respectively:
In Eq. (3), singular terms \({\delta }_{{S}_{n+}}{{\bf{n}}}_{{S}_{n}}\times {{\bf{E}}}_{n},{\delta }_{{S}_{n}}{{\bf{n}}}_{{S}_{n}}\times {{\bf{H}}}_{n},{\delta }_{{S}_{n+}}{{\bf{n}}}_{{S}_{n}}\cdot {{\bf{B}}}_{e}\), and \({\delta }_{{S}_{n}}{{\bf{n}}}_{{S}_{n}}\cdot {{\bf{D}}}_{n}\) represent the coupling of electromagnetic radiation at the surface of separation between the spaces \({\Omega }_{n}\) and \({\Omega }_{e}\). The remaining singular terms, conversely, define an appropriate set of boundary conditions. In the case of zero electromagnetic field inside each resonator, \({{\bf{E}}}_{n}=0\) 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 \({\delta }_{{S}_{n+}}{{\bf{n}}}_{{S}_{n}}\times {{\bf{E}}}_{e}={\delta }_{{S}_{n+}}{{\bf{n}}}_{{S}_{n}}\cdot {{\bf{B}}}_{e}\mathrm{=0}\), which generates the following set of boundary conditions on \({S}_{n}\) in the limit \(\varepsilon \to 0\):
Equation (5) show that that the whole resonator space \({\Omega }_{r}={\sum }_{n}{\Omega }_{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 \({\varOmega }_{e}\):
Equation (6) imply that the external space is seen from within each resonator region \({\Omega }_{n}\) as a Perfect Magnetic Conductor (PMC) material. Boundary conditions (5) and (6) lead to the following final set of Maxwell equations:
In Eq. (7) we have expanded the electric displacement \(\varepsilon ({\bf{x}})\frac{\partial }{\partial t}{{\bf{E}}}_{n}={\varepsilon }_{0}{\varepsilon }_{n}({\bf{x}})\frac{\partial }{\partial t}{{\bf{E}}}_{n}+{{\bf{J}}}_{\Delta }\) into a linear contribution \({\varepsilon }_{n}({\bf{x}}){{\bf{E}}}_{n}\) and a generic source term \({{\bf{J}}}_{\Delta }({\bf{x}},t)\) 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 \({\Omega }_{n}\):
and terminated by PMC boundary conditions. The operator \(\{\nabla \}\times \) is selfadjoint with PMC boundary conditions^{47}. This implies that the resonator modes \({{\bf{E}}}_{nm}\), \({{\bf{H}}}_{nm}\) are orthogonal, form a complete basis and possess a real frequency \({\omega }_{m}\). Mode orthogonality is calculated from the eigenvalue problem (8) using standard techniques^{48} and occurs through the following relationship:
When \(m=m{\prime} \), the integral expression in (9) represents the electromagnetic energy stored inside the resonator space \({\Omega }_{n}\) by the mth 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 \(\langle \varepsilon (t)\rangle \) dissipated inside the resonator space for monochromatic excitation at frequency \(\omega \) becomes simply expressed as the sum of the energy of each mode:
with \({a}_{m}(t)=\tilde{a}(\omega ){e}^{i\omega t}\) and \(\tilde{{\bf{a}}}=[{\tilde{a}}_{1}(\omega ),\ldots ,{\tilde{a}}_{M}(\omega )])]\) defining the vector of amplitudes of the internal modes in the frequency domain. Equation (11) is the counterpart of the expression of the power in multimode waveguides, furnished by the squared sum of noninteracting terms (see, e.g., Eq. (2.2.56) of^{27}). The expression (11) is only exact in the formulation of Eq. (8) with PMC boundary conditions.
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 \({\Omega }_{n}\) characterized by a generic cuboid volume with sides \({L}_{x},{L}_{y},{L}_{z}\) along the \((x,y,z)\) axes and filled with a dielectric material with refractive index \(n({\bf{r}})\), the frequencies of internal modes are:
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 \({\Omega }_{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 \({{\bf{E}}}_{e}({\bf{x}},t)={\tilde{{\bf{E}}}}_{e}({\bf{x}},\omega ){e}^{i\omega t},{{\bf{H}}}_{e}({\bf{x}},t)={\tilde{{\bf{H}}}}_{e}({\bf{x}},\omega ){e}^{i\omega t}\):
with time varying amplitude coefficients \({s}_{m\pm }(t)={\tilde{s}}_{m\pm }(\omega ){e}^{i\omega t}\), which describe the time evolution of incoming \({{\bf{E}}}_{m+},{{\bf{H}}}_{m+}\) and outgoing \({{\bf{E}}}_{m},{{\bf{H}}}_{m}\) waves through \(m=\mathrm{1,}\ldots ,M\) different scattering channels. Traveling waves \({{\bf{E}}}_{m\pm },{{\bf{H}}}_{m\pm }\) depend in general on \(\omega \) through their wavevector \({\bf{k}}\) as, e.g., in the case of plane waves \({e}^{\pm i{\bf{k}}\cdot {\bf{r}}}\), spherical waves \({e}^{\pm ikr/r}\) or other types of traveling waves in the freespace.
Following the same idea developed for the internal modes expansion in Eq. (11), we normalize radiating modes \({\tilde{{\bf{E}}}}_{m\pm }\) and \({\tilde{{\bf{H}}}}_{m\pm }\) through an observable quantity of physical interest. We here use the optical power^{49}, defined from the following integral when \(h=h{\prime} \):
with \(S={\sum }_{n}{S}_{n}\) representing the union of all the surfaces of the resonator space \({\Omega }_{r}\). With the orthogonality condition (15), the effective power P flowing through S assumes the expression:
with \({\tilde{{\bf{s}}}}_{\pm }=[{\tilde{s}}_{1\pm }(\omega ),\ldots ,{\tilde{s}}_{M\pm }(\omega )]\)) 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 spatiotemporal coupled mode equations, which relate the time evolution of the amplitudes of internal modes \({a}_{m}(t)\), with outgoing scattered waves \({s}_{h}(t)\) for a given set of impinging sources \({s}_{h+}(t)\). 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\pm }\) 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 \({\varOmega }_{e}\), due to the absence of any source \({{\bf{J}}}_{\Delta }=0\), Maxwell’s equations are linear and the scattered modes \({s}_{h}(t)\) follow a linear evolution as a function of modes \({a}_{m}(t)\) and impinging fields \({s}_{h+}(t)\). The time dynamics of the scattered field \({{\bf{s}}}_{}(\omega )\) in the frequency domain is then expressed as a linear superposition of \(\tilde{{\bf{a}}}(\omega )\) and \({\tilde{{\bf{s}}}}_{+}(\omega )\):
with \(\tilde{D}(\omega ),\tilde{C}(\omega )\) being linear matrices. To write the equations describing the dynamics of \({\bf{a}}(t)\), we first consider the case of linear materials with \({{\bf{J}}}_{\Delta }=0\). In this limit, Maxwell equations inside the resonator space \({\varOmega }_{n}\) are also linear, and the dynamics of \({\bf{a}}(t)\) follow from the most general form of the linear time evolution equation for the modes \({a}_{m}(t)\), with input sources corresponding to impinging waves \({{\bf{s}}}_{+}\):
In Eq. (18), \(\dot{{\bf{a}}}={\rm{d}}{\bf{a}}/{\rm{dt}}\), and \(H(t),K(t)\) are linear matrices with Fourier pairs \(\tilde{H}(\omega ),\tilde{K}(\omega )\) in the frequency space. The matrix \(H\) has no particular symmetry, and as such it can be decomposed as \(H=iW+\Gamma \), with a SkewHermitian matrix \({\rm{iW}}\) and a Hermitian matrix \(\Gamma \). 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 \(\tilde{H},\tilde{K},\tilde{C},\tilde{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 selfconsistency relations:
where \(\tilde{\Gamma }(\omega )\) is the Fourier transform of \(\Gamma \). Equation (19) are a particular form of the fluctuation dissipation theorem^{50} applied for Maxwell equations. The conditions imposed by Eq. (20) solve Eqs. (17) and (18) in the frequency domain:
with \(\mathrm{1/}\tilde{X}\) being shorthand notation for the inverse matrix \({\tilde{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\), \(\varGamma \) are frequency independent and \({a}_{m}(t)\) 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, \(\tilde{{\bf{a}}}=0\) and the system output is characterized by the open loop response \({\tilde{{\bf{s}}}}_{}=\tilde{C}(\omega )\cdot {\tilde{{\bf{s}}}}_{+}\). This is the contribution that arises purely from propagation effects and away from any resonant lightmatter interaction. When \(\tilde{{\bf{a}}}\ne 0\), the system response is characterized by a second term represented by the closedloop 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: \(\tilde{K}\), \(\tilde{C}\), and \(W\). The physical meaning of these matrices and their expressions is analyzed next. When \(\tilde{K}=0\) in Eq. (21), the nth resonator \({\Omega }_{n}\) and the external \({\Omega }_{e}\) space are uncoupled, and the dynamics of (21) reduce to:
The first equation describes the undamped motion of the internal modes \({a}_{m}(t)\) at frequencies \({\omega }_{m}\) arising from the diagonal elements of the (resonances) matrix \({W}_{mm{\prime} }={\omega }_{m}{\delta }_{mm{\prime} }\) (Fig. 1b). This dynamics represent free oscillations of noninteracting, orthogonal modes of the resonator space solution of Eq. (17).
In the external space \({\Omega }_{e}\), as obtained from the second of Eq. (22), light dynamics reduce to a scattering process of input sources \({{\bf{s}}}_{+}\) impinging on the PMC material defined in each \({\Omega }_{n}\) representing a resonator. The matrix \(\tilde{C}(\omega )\) is the unitary scattering matrix describing this process. The offdiagonal terms \({\tilde{C}}_{mn}=\frac{{\tilde{s}}_{m}}{{\tilde{s}}_{n+}}\) of the matrix \(\tilde{C}(\omega )\) represent the scattering of energy from incoming waves \({\tilde{s}}_{m+}\) into outgoing radiation on different modes \({\tilde{s}}_{m}\), while the diagonal terms \({\tilde{C}}_{nn}=\frac{{\tilde{s}}_{n}}{{\tilde{s}}_{n+}}\) are the reflection coefficients of incoming waves \({\tilde{s}}_{n+}\) into contributions propagating in the same mode but in opposite directions \({\tilde{s}}_{n}\).
The scattering matrix can be expressed in exponential form \(\tilde{C}(\omega )={e}^{i\varphi (\omega )}\), with \(\phi \) being \(M\times M\) matrix. 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:
Equation (23) does not alter the space partitioning \({\Omega }_{n},{\Omega }_{e}\), nor the mode \({a}_{n}(t)\) evolution. The transformation (23) defines a new set of scattering modes via (14) that diagonalise the scattering matrix and, as such, provide only reflections for each input channel excited in the dynamics. An example of this representation is furnished in the next section.
When the spaces \({\Omega }_{e}\) and \({\Omega }_{n}\) interact with nonzero couplings \(\tilde{K}\), electromagnetic energy flows from the cavity region to \({\Omega }_{e}\), and viceversa. The coupling matrix \(\tilde{K}\) has in general a small, or weak dependence on the frequency \(\omega \). 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:
with \({{\bf{E}}}_{n}\) being the field inside the resonator space and \({{\bf{H}}}_{m\pm }^{\ast }\) the scattering modes in the mth channel. The time average of the integral gives a nonzero contribution only when the fields \({{\bf{H}}}_{m\pm }^{\ast }\) and \({{\bf{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}^{\pm i{\bf{k}}\cdot {\bf{r}}}\) interacting with cuboid resonator structures \({\Omega }_{n}\) with resonant wavevectors \({{\bf{k}}}_{lpq}\) represented by Eq. (12). For a general frequency value \(\omega =c{\bf{k}}\) that is not resonant with any internal resonance \({\bf{k}}\ne {{\bf{k}}}_{lpq}\), the integral (24) is characterized by oscillatory terms of the type \({e}^{i({\bf{k}}{{\bf{k}}}_{lpq})\cdot {\bf{r}}}\), which integrated through S do not furnish contribution. In all of these cases, the coupling matrix \(\tilde{K}\) is frequency independent.
In the situation where \({{\bf{J}}}_{\Delta }\ne 0\), the contribution of \({{\bf{J}}}_{\Delta }\) 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:
with \({{\bf{j}}}_{\Delta }=[{j}_{\varDelta 1}(t),\ldots ,{j}_{\varDelta m}(t)]\) being a vector of projected source terms with general contribution:
equations (25) and (26) model the exact dynamics of light matter interactions in multimodal material structures with arbitrarily defined linear and nonlinear responses. All nonlinear responses and dispersive effects are taken into account by specializing the form of the generic source term \({{\bf{J}}}_{\Delta }({\bf{x}},t)\).
For dispersive material responses including complex distribution of absorption in the frequency domain \(\omega \), the source term can be expressed as \({{\bf{J}}}_{\Delta }={\varepsilon }_{0}\frac{\partial {{\bf{P}}}_{L}}{\partial t}\), with:
with \(\chi ({\bf{x}},t)\) being the noninstantaneous susceptibility, which can be further expressed in the frequency domain \(\tilde{\chi }({\bf{x}},\omega )\) with a Lorentz multipole expansion^{52}:
where \(\Delta {\varepsilon }_{k}\) is the difference between zerofrequency and infinity frequency relative permittivities, \({\omega }_{k}\) is the frequency of the pole and \({\delta }_{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 \({\delta }_{k}\).
For nonlinear effects, such as Kerr and other types of perturbative nonlinearities, the source term can be expanded as \({{\bf{J}}}_{\Delta }={\varepsilon }_{0}\frac{\partial {{\bf{P}}}_{NL}}{\partial t}\), with^{53}:
where \({\chi }^{(m)}({\bf{x}},t)\) is the mth order noninstantaneous susceptibility tensor. The substitution of Eqs. (27–29) into Eq. (25) originate a system of dispersive nonlinear equations in the resonator modes \({\bf{a}}\) that can be used as the basis to study complex dynamics of nonlinear and/or dispersive lightmatter 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}:
where \({S}_{i}\) (\(i=\mathrm{1,}\,\mathrm{4,}\,9\)) are components of the atomic Bloch coherence vector \({\bf{S}}=[{S}_{1},{S}_{2},\ldots ,{S}_{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 \({\bf{x}}\), \({\bf{y}}\), and \({\bf{z}}\) are unit vectors along the \(x\), \(y\) and \(z\) axes, respectively. The time evolution of the coherence \({\bf{S}}\) vector follows from Bloch theory:
with \(\Gamma ({{\bf{E}}}_{n})\) being the nonlinear coupling matrix function of the electric field \({{\bf{E}}}_{n}({\bf{x}},t)\) inside the \(n\)th resonator, \(\gamma \) the damping matrix and \({{\bf{S}}}^{\mathrm{(0)}}\) the initial state of \({\bf{S}}\). The derivation of this equation with detailed expressions for \(\Gamma \) and \(\gamma \) is found in^{54}. Equations (25)–(26) and (30)–(31) form a nonlinear system of coupled mode equations capable of describing generic processes of gain, such as stimulated emission and lasing effects, in complex optical resonators.
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:
and provides the number of eigenstates \({\omega }_{n}\) in the frequency interval \([\omega ,\omega +d\omega ]\). 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 pointdipole source \({{\bf{S}}}_{p}\) in the electric \({\bf{E}}\) or magnetic \({\bf{H}}\) field such as:
with \({\hat{{\bf{e}}}}_{l}\) (\(l=x,y,z\)) being a unit vector along one coordinate axis, \({{\bf{x}}}_{0}\) the coordinates of a point inside the material whose DOS is to be computed, and \(p(t)=\delta (t{t}_{0})\) being a short time pulse with broadband spectrum.
The local density of states (LDOS) measured at the point \({{\bf{x}}}_{0}\) along the axis \(l=x,y,z\) is obtained from the power density spectrum of the electric or magnetic field measured at \({{\bf{x}}}_{0}\) via the following relation^{56}:
once the LDOS is calculated, the DOS is obtained by integrating over the volume of the material V and by summing up the contributions arising from different axis:
Equation (35) is rarely employed in practice due to the requirement of performing a large number of simulations, in principle one for each point \({{\bf{x}}}_{0}\) and polarization considered.
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 \({\Omega }_{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 \({\Omega }_{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}_{p}({\bf{x}},t)\) centered at any point outside \({\Omega }_{n}\), and then measuring the time evolution of the internal modes \({a}_{m}(t)\) of the resonator space \({\Omega }_{n}\) by projecting the electromagnetic field \({{\bf{E}}}_{n}\) or \({{\bf{H}}}_{n}\) over the corresponding eigenmode \({{\bf{E}}}_{nm}\) or \({{\bf{H}}}_{nm}\) via the orthogonality relations (9):
Equation (37) can be evaluated during the FDTD simulation, as the electric field \({{\bf{E}}}_{n}({\bf{x}},t)\) is available at each time \(T\), and modal distributions \({{\bf{E}}}_{nm}({\bf{x}})\) are easily calculated for resonators terminated by ideal PEC/PMC boundary conditions by eigenvalue solvers^{57}. Once the distribution of \({a}_{m}(t)\) is known, the DOS is directly computed with Eq. (36).
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 \({\Omega }_{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 \({{\bf{E}}}_{n}\), \({{\bf{H}}}_{n}\) inside a resonator structure \({\Omega }_{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 \({\Omega }_{n}\) with permittivity \({\varepsilon }^{\mathrm{(0)}}({\bf{x}})\) and permeability \({\mu }^{\mathrm{(0)}}({\bf{x}})\) calculated by using the same set of PEC/PMC boundary conditions.
By projecting Eq. (38) on the new set of modes \({{\bf{E}}}_{nm}^{\mathrm{(0)}}({\bf{x}})\), \({{\bf{H}}}_{nm}^{\mathrm{(0)}}({\bf{x}})\), we obtain a new set of time varying amplitudes \({b}_{m}(t)\) related to \({a}_{m}(t)\) as follows:
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 \({\varepsilon }^{\mathrm{(0)}}({\bf{x}})={\varepsilon }_{0}{\varepsilon }_{r}\) and \({\mu }^{\mathrm{(0)}}={\mu }_{0}\), which are analytically expressed by Eqs. (12) and (13) via simple trigonometric formulas. Other possible simple choices are represented by spherical or cylindrical \({\Omega }_{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 \({\Omega }_{1}\) (Fig. 3a green area). We illuminate the structure by a single broadband pulse source (Fig. 3 input pulse) and calculate the amplitudes \({a}_{m}(t)\) of the internal modes by Eq. (37). Figure 3b,c shows the time evolution of \({a}_{m}(t)\) and \(\tilde{a}(\omega {)}^{2}\) for the first modes \(m=\mathrm{1,}\,\mathrm{2,}\,3\). The resulting DOS, calculated from Eq. (36), is reported in Fig. 3d. The energy density distribution from FDTD simulations at the time \(t\cdot \frac{c}{2\pi d}\mathrm{=2}\) shown on Fig. 3e with red dotted rectangle representing a cuboid resonator region \({\Omega }_{1}\).
Photonic resonance networks
A particularly important quantity in the analysis of resonant systems is the mode quality factor \(Q=\frac{{\omega }_{0}\tau }{2}\), defined as the product between the mode frequency \(\omega \) and the mode decay rate \(\tau \), 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 \(\omega \) and \(\tau \) from the mode decaying rates in time, or from the mode frequency linewidth \(\Delta \) in the LDOS. This approach assumes noninteracting 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 \(\omega \) is away from a resonant frequency of the system:
with \({{\bf{a}}}_{0}={\bf{a}}(t=\mathrm{0)}\) being the initial (excited) state and \({e}^{H(t)}\) the matrix exponential. The matrix exponential can be expressed in closed form by diagonalizing \(H=Q\Lambda {Q}^{1}\), with \({\Lambda }_{mm{\prime} }={\sigma }_{m}{\delta }_{mm{\prime} }\) being the diagonal matrix of complex eigenvalues \({\sigma }_{m}=i{\omega }_{m}{\gamma }_{m}\):
with \({[{e}^{\Lambda {\rm{t}}}]}_{mm{\prime} }={e}^{{\sigma }_{m}t}{\delta }_{mm{\prime} }\). If we launch an impulsive source \({{\bf{s}}}_{h+}(t)=\delta (t){\delta }_{hh{\prime} }\) on a single scattering channel \(h{\prime} \), by substituting Eq. (42) into (41) we obtain the mode solution \({a}_{m}(t)\) at \(t\, > \,0\):
expressed as the sum of complex damped exponentials with \({\alpha }_{mm{\prime} }\) constant coefficients arising from the matrix product of \({[Q{e}^{\varLambda t}{Q}^{1}]}_{mm{\prime} }{a}_{m{\prime} }\mathrm{(0)}\). The Fourier transform of the mode \({a}_{m}(\omega )\) is a complex rational function:
with poles \({\omega }_{m{\prime} }+i{\gamma }_{m{\prime} }\). The corresponding DOS is also a rational function:
with poles \({s}_{m}={\omega }_{m}+i\frac{1}{{\tau }_{m}}\). To extract the Q factor, we proceed as follows. The time varying amplitude \({a}_{m}(t)\) of the electromagnetic field oscillating at frequency \({\omega }_{m}\) and decaying constant \({\tau }_{m}=\frac{2Q}{{\omega }_{m}}\) of an internal mode is \({a}_{m}(t)={a}_{0}{e}^{i({\omega }_{m}t\frac{t}{{\tau }_{m}})}\), and generates a contribution to the DOS equal to:
with \({s}_{m}={\omega }_{m}+i\frac{1}{{\tau }_{m}}\). By equating Eqs. (45) and (46) the quality factor \({Q}_{m}\) associated to the resonant mode at \({\omega }_{m}\) is:
with \( {\mathcal R} \) and \( {\mathcal I} \) being the real and imaginary parts of \({s}_{m}\), respectively. The calculation of the network of mode quantities \({Q}_{m}\), \({\omega }_{m}\) and \({\tau }_{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 \((\omega ,\gamma )\) 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).
Complete representation of resonant modes
Once the mode evolutions are obtained and stored in the \({a}_{m}(t)\) coefficients, it is possible to obtain the expression of each resonant mode from Eq. (10) in space (\({\bf{x}}\)), time (\(t\)) and frequency (\(\omega \)) after transforming the mode amplitudes \({a}_{m}(t)\) in the spectral domain. The main advantage of this approach lies in the fact that the quantities \({a}_{m}(t)\) are computed from a single FDTD simulation with the same setting used for the calculation of the DOS, and the complete spectrum of modes is directly available from the DOS via Eq. (36).
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 setup. The resonator region \({\Omega }_{1}\) (Fig. 5a, blue region) is composed of a cuboid with thickness d along the propagation axis \(z\), and with infinite sides along \(x\) and \(y\). The resonator space is filled with a uniform dielectric material of refractive index \(n=3.5\). Despite its simplicity, this structure is sufficiently general to allow a detailed discussion of many important properties of photonics networks.
The system of Fig. 5a has two scattering channels: when only source \({s}_{1+}\) is active, the reflection \(R\) is measured in \({s}_{1}\) and the transmission \(T\) in \({s}_{2}\). Conversely, when source \({s}_{2+}\) is launched on the structure, its reflection \(R\) is retrieved in \({s}_{2}\) and the transmission \(T\) in \({s}_{1}\). As the dielectric slab is symmetric along \(z\), only one case (\({s}_{1+}\) or \({s}_{1}\) active) is sufficient to calculate the material response.
Following Eq. (13), the frequencies of the internal modes are \({\omega }_{m}=\frac{c\pi m}{nd}\), and the spatial distribution of the magnetic modes reduce to \({H}_{rl}={A}_{0}sin\left(\frac{l\pi z}{d}\right)\) polarized either along \(x\) or \(y\). External modes, conversely, are represented by incoming and outgoing plane waves \({e}^{\pm ikz}\):
We then express the scattering matrix for the 2 × 2 system of Fig. 5:
in which \({\tilde{C}}_{12}={\tilde{C}}_{21}=0\) and \({\tilde{C}}_{11}={\tilde{C}}_{22}=1\) arise from PEC boundary condition at the resonator space \({\Omega }_{n}\), originating reflections for each incoming source when light is injected from the external space \({\Omega }_{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 \(\Gamma \) from the coupling coefficients \({K}_{ij}={K}_{ij}{e}^{i{\theta }_{ij}}\):
the symmetry of the resonator structure along \(z\) implies that damping factors along channels 1 and 2 are the same: \({K}_{i1}={K}_{i2}=\sqrt{{\Gamma }_{ii}}\). 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 antisymmetrically in the scattering channels:
these relations imply that when internal modes \(i\), \(j\) possess opposite symmetry along \(z\), we have \({\theta }_{i1}{\theta }_{j1}=({\theta }_{i2}{\theta }_{j2})\) and \({\Gamma }_{ij}=0\). Conversely, when modes \(i\), \(j\) have the same symmetry \({\theta }_{i1}{\theta }_{j1}=({\theta }_{i2}{\theta }_{j2})\) and \({\Gamma }_{ij}=\pm \sqrt{{\Gamma }_{ii}{\Gamma }_{jj}}\), with the plus sign if both modes are even and with minus sign if modes are odd. These cases are summarized as follows:
by substituting the expression of \({\Gamma }_{ij}\) in the coupled mode Eq. (41), we find the following main result, which holds for all modes having the same parity:
this equation states that the amplitude of each internal mode \({\tilde{a}}_{m}({\omega }_{m})\), evaluated at the resonant frequency \({\omega }_{m}\) of the mode, furnishes the amplitude of the damping factor \({\varGamma }_{ll}\), and goes to zero at all the resonance frequencies of the other modes. Equation (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}_{n}(t)\) 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 \([\frac{\omega d}{2\pi c}]\) frequency resolution. To achieve this requirement we launch simulations for 50000 timesteps, dt = 5.87 ps, and corresponding frequency resolution of 0.0015 \([\frac{\omega d}{2\pi c}]\). Figure 5(b,c) shows the spectral distribution \({{\tilde{a}}_{m}(\omega )}^{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 \({\omega }_{l\ne m}\) of the other modes with the same symmetry along \(z\).
From the amplitude intensity at the mode internal frequency \({{\tilde{a}}_{l}(\omega )}^{2}\), we apply (54) and calculate all the elements of the coupling \(K\) and damping \(\varGamma \) 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 \({\tilde{a}}_{m}(\omega )\) 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 \({\gamma }_{m}\) and resonant frequencies following the exact analytic formula \(\frac{{\omega }_{m}d}{2\pi c}=\frac{m}{2\sqrt{{\varepsilon }_{r}}}\) (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 \(\omega \) due to the increasing resonant frequency \({\omega }_{m}\) of each mode.
Two and three dimensional structures
We apply the STCMT approach to describe more complicated geometries for both nonperiodic (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 \({\Omega }_{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 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.
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 \({\Omega }_{1}\) used in the projections of the time varying coefficients \({a}_{m}(t)\).
The corresponding energy distributions of the resonant modes are portrayed in Figs. 9c,d and 10b,c. They represent whisperinggallery modes of the ringlike geometry in the nonperiodic case, and more complex FloquetBloch modes in the periodic structure.
Figure 11a–d shows the calculation results for a three dimensional resonator structure, delimited by a cubical \({\Omega }_{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}_{1}=52\) and \({Q}_{2}=21\).
Conclusion
In this work we formulate an exact spatiotemporal 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.
Data availability
Data and a source code used in this work are available upon a reasonable request.
References
 1.
Zheludev, N. I. & Kivshar, Y. S. From metamaterials to metadevices. Nature Materials 11. Review Article (Oct. 2012).
 2.
Jahani, S. & Jacob, Z. Alldielectric metamaterials. Nature Nanotechnology 11. Review Article (Jan. 2016).
 3.
Kivshar, Y. Alldielectric metaoptics and nonlinear nanophotonics. National Science Review 5,144–158 (Jan. 2018).
 4.
Mayer, B. et al. Monolithically Integrated HighβNanowire Lasers on Silicon. Nano Letters 16, 152–156 (Jan. 2016).
 5.
Miao, P. et al. Orbital angular momentum microlaser. Science 353, 464–467 (2016).
 6.
Ma, Z. et al. Terahertz AllDielectric Magnetic Mirror Metasurfaces. ACS Photonics 3, 1010–1018 (June 2016).
 7.
Yang, Z.J. et al. Dielectric nanoresonators for light manipulation. Physics Reports 701, 1–50 (July 2017).
 8.
Pilozzi, L. & Conti, C. Topological lasing in resonant photonic structures. Phys. Rev. B 93, 195317 (19 May 2016).
 9.
Xu, Z., Song, W. & Crozier, K. B. Direct observation of optical trapping of a single quantum dot with an allsilicon nanoantennain Frontiers in Optics 2017 (Optical Society of America, 2017), FM3B.2.
 10.
Shibanuma, T., Grinblat, G., Albella, P. & Maier, S. A. Efficient Third Harmonic Generation from MetalDielectric Hybrid Nanoantennas. Nano Letters 17, 2647–2651 (Apr. 2017).
 11.
Koshelev, K., Favraud, G., Bogdanov, A., Kivshar, Y. & Fratalocchi, A. Nonradiating photonics with resonant dielectric nanostructures. Nanophotonics 8, 725–745 (Mar. 2019).
 12.
Papasimakis, N., Fedotov, V. A., Savinov, V., Raybould, T. A. & Zheludev, N. I. Electromagnetic toroidal excitations in matter and free space. Nature Materials 15, 263–271 (Feb. 2016).
 13.
Miroshnichenko, A. E. et al. Nonradiating anapole modes in dielectric nanoparticles. Nature Communications 6 (Aug. 2015).
 14.
Plotnik, Y. et al. Experimental Observation of Optical Bound States in the Continuum. Physical Review Letters 107 (Oct. 2011).
 15.
Kodigala, A. et al. Lasing action from photonic bound states in continuum. Nature 541, 196–199 (Jan. 2017).
 16.
Miroshnichenko, A. E., Malomed, B. A. & Kivshar, Y. S. Nonlinearly PTsymmetric systems: Spontaneous symmetry breaking and transmission resonances. Phys. Rev. A 84, 012123 (1 July 2011).
 17.
Hodaei, H., Miri, M.A., Heinrich, M., Christodoulides, D. N. & Khajavikhan, M. Paritytime–symmetric microring lasers. Science 346, 975–978 (2014).
 18.
Hsu, C. W., Zhen, B., Stone, A. D., Joannopoulos, J. D. & Soljačić, M. Bound states in the continuum. Nature Reviews Materials 1 (July 2016).
 19.
Kupriianov, A. S. et al. Metasurface engineering through bound states in the continuum 2019. eprint:arXiv:1904.04688.
 20.
Totero Gongora, J. S., Miroshnichenko, A. E., Kivshar, Y. S. & Fratalocchi, A. Anapole nanolasers for modelocking and ultrafast pulse generation. Nature Communications 8 (May 2017).
 21.
Huang, Y.W. et al. Toroidal Lasing Spaser. Scientific Reports 3 (Feb. 2013).
 22.
Rybin, M. V. et al. HighQ Supercavity Modes in Subwavelength Dielectric Resonators. Phys. Rev. Lett. 119, 243901 (24 Dec. 2017).
 23.
Bohn, J. et al. Active Tuning of Spontaneous Emission by MieResonant Dielectric Metasurfaces. Nano Letters 18, 3461–3465 (June 2018).
 24.
Bohren, C. & Huffman, D. R. Absorption and Scattering of Light by Small Particles (Wiley Science Paperback Series, 1998).
 25.
Yeh, C. & Shimabukuro, F. I. The Essence of Dielectric Waveguides 179–219 (Springer US, Boston, MA, 2008).
 26.
Agrawal, G. Nonlinear Fiber Optics (Elsevier Science, 2007).
 27.
Kogelnik, H. Theory of Dielectric Waveguides 13–81 (Springer Berlin Heidelberg, Berlin, Heidelberg, 1975).
 28.
Huang, W.P. Coupledmode theory for optical waveguides: an overview. J. Opt. Soc. Am. A 11, 963–983 (Mar. 1994).
 29.
Yariv, A. Coupledmode theory for guidedwave optics. IEEE Journal of Quantum Electronics 9, 919–933 (Sept. 1973).
 30.
Marcuse, D. Theory of dielectric optical waveguides English, xi, 257 p. (Academic Press New York, 1974).
 31.
Tamir, T. & Garmire, E. Integrated optics (Springer, 1979).
 32.
Haus, H. A. Waves and Fields in Optoelectronics 402 (Prentice Hall, London, 1983).
 33.
Suh, W., Wang, Z & Fan, S. Temporal coupledmode theory and the presence of nonorthogonal modes in lossless multimode cavities. IEEE Journal of Quantum Electronics 40, 1511–1518 (Oct. 2004).
 34.
Liu, C. et al. Enhanced energy storage in chaotic optical resonators. Nature Photonics 7, 473–478 (Jan. 1, 2013). published.
 35.
Liu, C., Falco, A. D. & Fratalocchi, A. Dicke phase transition with multiple superradiant states in quantum chaotic resonators. Physical Review X 4 (Jan. 1, 2014). published.
 36.
Gomard, G., Peretti, R., Drouard, E., Meng, X. & Seassal, C. Photonic crystals and optical mode engineering for thin film photovoltaics. Opt. Express 21, A515–A527 (May 2013).
 37.
Shcherbakov, M. R. et al. Photon acceleration and tunable broadband harmonics generation in nonlinear timedependent metasurfaces. Nature Communications 10, 1345 (2019).
 38.
Joannopoulos, J. D., Johnson, S. G., Winn, J. N. & Meade, R. D. Photonic Crystals: Molding the Flow of Light (Second Edition) 2nd ed. (Princeton University Press, 2008).
 39.
Galinski, H., Fratalocchi, A., Döbeli, M. & Capasso, F. Light Manipulation in Metallic Nanowire Networks with Functional Connectivity. Advanced Optical Materials 5.cited By 4 (Jan. 1, 2017). published.
 40.
Hackenbroich, G., Viviescas, C. & Haake, F. Field Quantization for Chaotic Resonators with Overlapping Modes. Phys. Rev. Lett. 89, 083902 (8 July 2002).
 41.
Viviescas, C. & Hackenbroich, G. Field quantization for open optical cavities. Phys. Rev. A 67, 013805 (1 Jan. 2003).
 42.
Antenucci, F., Conti, C., Crisanti, A. & Leuzzi, L. General Phase Diagram of Multimodal Ordered and Disordered Lasers in Closed and Open Cavities. Phys. Rev. Lett. 114, 043901 (4 Jan. 2015).
 43.
Gongora, J. S. T., Favraud, G. & Fratalocchi, A. Fundamental and highorder anapoles in alldielectric metamaterials via Fano–Feshbach modes competition. Nanotechnology 28, 104001 (Feb. 2017).
 44.
Cao, H. Review on latest developments in random lasers with coherent feedback. Journal of Physics A: Mathematical and General 38, 10497–10535 (Nov. 2005).
 45.
Kanwal, R. P. Generalized Functions Theory and Applications (Springer, London, 2004).
 46.
Marcucci, G. & Conti, C. Irreversible evolution of a wave packet in the riggedHilbertspace quantum mechanics. Phys. Rev. A 94, 052136 (5 Nov. 2016).
 47.
Jackson, J. D. Classical electrodynamics 3rd ed. (Wiley, New York, NY, 1999).
 48.
Sakurai, J. J. & Napolitano, J. Modern Quantum Mechanics 2nd ed. (Cambridge University Press, 2017).
 49.
Tamir, T. Integrated Optics (Springer, Berlin, 1975).
 50.
Kubo, R. The fluctuationdissipation theorem. Reports on Progress in Physics 29, 255–284 (1966).
 51.
Fyodorov, Y. V. & Sommers, H.J. Statistics of resonance poles, phase shifts and time delays in quantum chaotic scattering: Random matrix approach for systems with broken timereversal invariance. Journal of Mathematical Physics 38, 1918–1981 (1997).
 52.
Landau, L. & Lifshitz, E. In Electrodynamics of Continuous Media (Second Edition) (eds. Landau, L. & Lifshitz, E.) Second Edition, 290–330 (Pergamon, Amsterdam, 1984).
 53.
Landau, L. & Lifshitz, E. In Electrodynamics of Continuous Media (Second Edition) (eds. Landau, L. & Lifshitz, E.) Second Edition, 372–393 (Pergamon, Amsterdam, 1984).
 54.
Fratalocchi, A., Conti, C. & Ruocco, G. Threedimensional ab initio investigation of lightmatter interaction in Mie lasers. Phys. Rev. A 78, 013806 (1 June 2008).
 55.
Economou, E. N. Green’s function in quantum physics 7 (Springer, 2006).
 56.
Taflove, A., Oskooi, A. & G. Johnson, S. Advances in FDTD Computational Electrodynamics: Photonics and Nanotechnology 74–78 (Artech House, Jan. 2013).
 57.
Logg, A. et al. Automated Solution of Differential Equations by the Finite Element Method (Springer, 2012).
 58.
Ito, S. & Nakatsukasa, Y. Stable polefinding and rational leastsquares fitting via eigenvalues. Numerische Mathematik 139, 633–682 (July 2018).
Acknowledgements
The authors acknowledge support from KAUST (OSR2016CRG52995) and Shaheen supercomputer from the Kaust Supercomputing Laboratory (KSL).
Author information
Affiliations
Contributions
M.M. and A.F. developed the theory and performed the simulations. F.G. and A.B.L. provided conceptual comments. M.M. and A.F. wrote the manuscript, with input from all the authors.
Corresponding author
Ethics declarations
Competing interests
The authors declare no competing interests.
Additional information
Publisher’s note Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Rights and permissions
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The images or other third party material in this article are included in the article’s Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the article’s Creative Commons license and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this license, visit http://creativecommons.org/licenses/by/4.0/.
About this article
Cite this article
Makarenko, M., BurgueteLopez, A., Getman, F. et al. Generalized Maxwell projections for multimode network Photonics. Sci Rep 10, 9038 (2020). https://doi.org/10.1038/s41598020652936
Received:
Accepted:
Published:
Further reading

Broadband vectorial ultrathin optics with experimental efficiency up to 99% in the visible region via universal approximators
Light: Science & Applications (2021)
Comments
By submitting a comment you agree to abide by our Terms and Community Guidelines. If you find something abusive or that does not comply with our terms or guidelines please flag it as inappropriate.