Exact Maxwell evolution equation of resonators dynamics: temporal coupled-mode theory revisited

Despite its widespread significance, the temporal coupled-mode theory (CMT) lacks a foundational validation based on electromagnetic principles and stands as a phenomenological theory relying on fitted coupling coefficients. We employ an ab initio Maxwellian approach using quasinormal-mode theory to derive an"exact"Maxwell evolution (EME) equation for resonator dynamics. While the resulting differential equation bears resemblance to the classical one, it introduces novel terms embodying distinct physics, suggesting that the CMT predictions could be faulted by dedicated experiments, for instance carried out with short and off-resonance pulses, or with resonators of sizes comparable to or greater than the wavelength. Nonetheless, our examination indicates that, despite its inherent lack of strictness, the CMT enables precise predictions for numerous experiments due to the flexibility provided by the fitted coupling coefficients. The new EME equation is anticipated to be applicable to all electromagnetic resonator geometries, and the theoretical approach we have taken can be extended to other wave physics.


Introduction
The temporal coupled-mode theory (CMT) is an acclaimed and widely used theoretical framework for modeling the continuous wave response and temporal dynamics of any integrated or free-space photonic resonant structure.Initially developed for the analysis and design of resonant systems, as well as the comprehension of energy coupling into and out of cavities and its exchange among various resonant modes [1][2][3][4], CMT has evolved over three decades to encompass advanced systems and contemporary physical phenomena.These include noninstantaneous nonlinearities, gain, PT-symmetry, 2D photonic materials, and time-reversal, underscoring its vast capabilities.
CMT holds the potential to simplify the description of resonant systems with arbitrary material compositions and geometric configurations.Even intricate systems can be treated as lumped harmonic oscillators.At its core, the dynamics of the resonator can be reduced to an ordinary differential equation with respect to time where |  | 2 represents the energy stored in the  ℎ mode of the resonator and  ̃ is the complex-valued frequency of the mode, encompassing radiative decay to excitation channels, free space and absorption.We intentionally use conventional notations in Eq. (1) and denote by | + ⟩ the vector of input amplitudes and |  ⟩ the vector of coupling constants between the mode and the channels.
The evolution equation ( 1) is typically accompanied by an additional linear equation that delineates the decay of the resonator into each individual channel.Generally, the coupling coefficients introduced by these additional equations are derived with multiple smart assumptions.We will not delve into them here, maintaining our focus on the fundamental dynamics of the resonator.
Equation (1) stands as a pivotal outcome of CMT, offering high intuitiveness, computational simplicity, and a clear delineation of the physics governing resonances.However, it should be regarded as a heuristic approach, lacking theoretical support for its validity.In standard usage, physical quantities characterizing the cavity (e.g., resonance frequency and coupling constants) are fitted, and their extraction is generally believed to be rigorously feasible through numerical or experimental means.
The primary motivation behind this work was to establish a rigorous foundation for Eq. ( 1) and derive an analytical expression for it directly from Maxwell equations.In this pursuit, we discovered that the evolution equation is only an approximation.Our objectives are thus to elucidate the approximations inherent in the CMT model, comprehend their consequences, and raise awareness when the differences become consequential.We contend that the approximations are significant enough to warrant a dedicated publication, building upon an analysis initially presented in a prior publication [5] by the same group.
Except for the study in [6], as far as we are aware, there is a lack of any Maxwellian derivation of the CMT in the existing literature.Furthermore, based on our comprehension of the analysis in [6], it appears that the Maxwellian evolution equation derived in our current work differs from that derived in [6].Also, the accuracy of the CMT evolution equation under specific conditions varies.
In Section 2, we initiate our exploration from Maxwell equations, undertaking an ab initio derivation of the differential equation governing the dynamics of resonators.This derivation leads us to an equation distinct from Eq. ( 1), which we shall refer to as the "exact" Maxwell evolution (EME) equation.The ability to derive this is facilitated by the advancements achieved over the past decade in electromagnetic quasinormal modes (QNMs) [7].
In Section 3, we highlight the two primary distinctions between the CMT evolution equation and the EME equation: the coupling coefficient is non-local, and the excitation is not directly proportional to the driving field, but to its temporal derivative.We also conduct a comprehensive assessment through a generic example, enabling us to discern the circumstances under which the predictions of the CMT and EME equations significantly differ.
In section 4, we quantitatively test the EME equation for a real optical system.With this example, we show that a term related to the temporal derivative of the driving field plays a significant role in determining the excitation of QNMs, and this term is not given by the CMT equation.
The insights gleaned from our conclusions in the realm of electromagnetism hold the potential to extend their influence on other types of waves for which QNM theory offers robust support [8].

The exact Maxwell evolution equation
The EME equation is derived by employing electromagnetic quasinormal mode (QNM) theory, which has experienced substantial advancements over the past decade and has now reached a level of maturity that allows its secure application to numerous contemporary electromagnetic challenges in photonics and plasmonics.Electromagnetic QNMs are source-free solutions of Maxwell equations, ∇ ×  ̃ = − ̃ 0  ̃ , ∇ ×  ̃ =  ̃( ̃)  ̃ , which satisfy the outgoing-wave condition for || → ∞.  ̃ and  ̃ respectively denote the normalized electric and magnetic fields [7],  denotes the possibly-dispersive permittivity tensors.The QNM field exponentially decays in time and Im( ̃) < 0 with the exp(−) notation.
Hereafter, we consider non-magnetic materials, with  0 denoting the permeability of vacuum.Extending the analysis to magnetic cases is straightforward [7].We also assume frequency-independent permittivity for simplicity, emphasizing that the distinctions we uncover between the evolution equations are unrelated to dispersion.The EME equation for resonators composed of Drude-Lorentz materials will be presented after the derivation of the EME equation.We will not derive it since it was previously established in [5].
Subsequently, we adopt the QNM theory with a complex-mapping regularization.This approach has produced highly compelling results, particularly in addressing intricate electromagnetic problems such as those involving branch cuts due to the presence of substrates [5].The concept of QNM regularization with complex mapping originated in quantum mechanics in the late 1970s.In the realm of electromagnetism, the favored mapping employs a complex scaling transformation known as perfectly matched layers (PMLs).This transformation maps the inherently divergent QNMs of open space onto square-integrable regularized QNMs [9].For an in-depth examination of QNM regularization, including the use of complex mapping, refer to Section 5.3 in [10].
The regularized QNMs exhibit numerous similarities with the normal modes of Hermitian systems.Specifically, they are square-integrable, can undergo normalization, and can be employed to expand the electromagnetic field   (, ) exp(−) scattered by a resonator illuminated by a monochromatic field oscillating at frequency    (, ) = ∑   () ̃()  , (2) where   is the complex-valued excitation coefficient of the  ℎ QNM.Throughout the manuscript,  denotes a bivector formed by electric and magnetic fields, for instance  ̃ = [ ̃,  ̃].
The extension of Eq. ( 2) is complete for  in the interior of compact resonators placed in free space [10,11].We will consider this special case (a simple 1D slab in air) in the final numerical test.For non-uniform backgrounds, e.g.layered substrate, there are branch cuts.To account for this possibility, we assume that the extension of Eq. ( 2) encompasses not only QNMs that are unaffected by the PML mapping but also an infinite set of "numerical" eigenvectors.These eigenvectors are computed as source-free solutions of Maxwell equations within the mapped space.Utilizing this expanded basis, highly accurate reconstructions of the field scattered by the resonator have been achieved, even for complex geometries such as resonators with dispersive materials positioned on substrates.Accurate reconstructions have been reported in both spectral [10,[12][13][14][15][16][17] or temporal [5,12] domains.
The QNM framework based a complex-mapping regularization provide a unique expression for the excitation coefficient [5,13] for the case of nondispersive materials under consideration here.This expression holds significant importance as the subsequent mathematical analysis depends on it.It's worth noting that all existing QNM frameworks, not just the one employed here, suggest an identical expression for non-dispersive materials, as indicated in Table 1 of reference [7].For dispersive materials  a case that will be addressed at the end of the Section  several expressions exist, see Table 1 in reference [7] and the comprehensive discussion in reference [16].
Except for the pole prefactor,   essentially represents an overlap integral between the QNM (or numerical modes) and the background field   = [  ,   ] in the absence of resonator.Δ() specifies the permittivity variation used for the scattered-field formulation, given by Δ() =   () −   (), where   () and   () denote the permittivities of the resonator and the background, respectively.Δ() ≠ 0 within a compact volume   that defines the resonator in the scattered-field formulation [7].
Unlike other rigorous methods featuring a finite number of QNMs and a background term [18,19], the intrinsic strength of the regularization approach for theoretical works  like the present one  lies in providing a unified mathematical framework.Within this framework, all eigenstates, QNMs or numerical modes, are treated identically.Owing to the analytical dependence of the excitation coefficient, the resonator response to any incident driving wavepacket can be analytically calculated using a spectral decomposition, summing over all distinct frequencies.Subsequently, the field   (, ) scattered by the resonator illuminated by a wavepacket can be expressed as a sum of contributions from QNMs and numerical modes [7]   (, ) = Re(∑   () ̃()  ), Here the overlap is now defined as   () = ∫ Δ()  (, ) •  ̃ 3    .Hereafter, the EME equation is derived from Eq. ( 5) utilizing complex analysis and generalized functionsdistributions  with a particular emphasis on Cauchy's integral theorem.Another demonstration, which is straightforward and explicit, is currently under elaboration [20].
For dispersive PMLs characterized by a complex scaling factor inversely proportional to the frequency, causality is preserved, and all eigenfrequencies of the regularized problem are situated in the lower half of the complex plane.As a result, all temporal excitation coefficients   () can be analytically computed utilizing contour integration and the Cauchy integral formula for the pole  ̃.
We recast the temporal integral in Eq. ( 5) in three terms, where the three integrands are strictly identical to that in Eq. ( 5) and are represented by '…' hereafter. is a small and positive quantity, and we will consider the limit  → 0 later on.As  → ∞, the exponential factor exp(( ′ − )) has different behaviors for these three terms.We will show later that , , and  can be computed analytically by wisely selecting the contour integration paths shown in Fig. 1.  we show that the modulus of the arc integral is bounded by a finite positive number.Thus, taking the limit  → 0 in the temporal integral, we obtain  1 ′ = 0.An identical same reasoning holds for  1 ′′ by considering the green arc.There is no pole contribution and the arc integral being bounded for the same reason, we have  1 ′′ = 0.
Altogether, we have  ̃() 3  with  an arbitrary 3 × 3 matrix.We intend to publish another demonstration of Eq. ( 8) at a later time, which will not involve the combination of Schwartz distributions and complex analysis.
We recall that Eq. ( 8) holds for non-dispersive materials.A similar derivation can be conducted for dispersive materials [5] using the auxiliary-field method [21,22] and an analytical expression for the permittivity.For instance, for resonators with a standard multi-pole Lorentz-Drude permittivity and for a dispersionless background in the volume   (  () may depend on the frequency for  ∉   ), it was previously shown [5] that where  ∞ denotes the high-frequency permittivity constant of the resonator.Equation ( 9) has been quantitatively tested for a gold bowtie antenna illuminated by a Gaussian pulse with a central frequency in the visible range, see Fig. 2 in [5].It will not be further tested hereafter.Equations ( 8) or (9) share similarities with the CMT evolution equation, yet differences are apparent, underscoring the incorporation of heuristics in the CMT.These heuristics are discussed in the next section.
Further deeper insights into Maxwellian derivation of the CMT regarding outcoupling with ports can be found in reference [6].However, it appears that Eqs. ( 8) or ( 9) deviate from Eq. ( 45) in reference [6], notably in terms of the final term, which incorporates the temporal derivative of the incident electric field   (, ).
Additionally, the following qualitative discussion regarding instances where the CMT and EME evolution equations diverge differs from the one presented in reference [6].

Qualitative discussion
We observe three primary distinctions between the EME equation and the classical CMT evolution equation:  In the CMT, the squared absolute value of the modal excitation coefficient, |  | 2 , signifies the energy stored in the mode  of the resonator.The energy stored in QNMs lacks meaning in a non-Hermitian context: the mode energy is not square integrable since the QNM fields exponentially diverge outside the resonator, see related discussions on the quality factor or mode volume in [7]. The second difference, unrelated to hermiticity, arises from the CMT assumption that the mode excitation by the wavepacket occurs locally [3,4].The local excitation approximation assumes that  + is proportional to the wavepacket electric field   (  , ) at an arbitrary point   typically chosen inside the resonator volume.In contrast, the QNM approach predicts a nonlocal excitation that spans the entire resonator volume through a spatial convolution of the incident wavepacket and QNM fields.Only in scenarios where the resonator is deep subwavelength, e.g.lump element circuits, does   () become directly proportional to the incident wavepacket. The third difference has not been predicted in prior works (except in [5]) and is unexpected: it highlights that the driving force is proportional to the derivative of the incident field.We do not precisely know the physical origin of this term.Nevertheless, we note that it implies that the excitation coefficient  responds instantaneously to the driving field, as observed with the last term in Eq. ( 7).We further note that this term is absent in the advanced formulation of Eq. ( 9) when dispersive materials are considered, specifically for  ∞ =   .For most materials, the high-frequency dielectric constant is anticipated to be approximately 1, indicating that materials behaves as if they were transparent or non-polarizable at high frequency.We have verified that, with the inclusion of background permittivity dispersion, this term vanishes for   ( → ∞) =   ( → ∞) ≈ 1 .Nonetheless, it is essential to highlight that the derivative term should be considered in usual simplified models assuming non-dispersive dielectrics.
The temporal CMT has been extensively utilized for several decades, consistently demonstrating its relevance in modeling the temporal dynamics of both integrated and freespace photonic resonant structures.The whole literature does not report any substantial deviations between CMT predictions and experimental or full-numerical data.Hereafter, we try to understand why.Specifically, we elucidate why the CMT evolution Eq. ( 1) can yield predictions in very good agreement with those provided by the EME equation, despite the inherent mathematical disparities between the two formulations.
In this section, the analysis is carried out qualitatively with the aim of providing a broad discussion applicable to various systems, such as small plasmonic cavities or ring resonators.The next Section will provide a quantitative examination for a Fabry-Perot cavity.

Slowly-varying envelop approximation.
The first approximation arises from the fact that incident wavepackets can often be described as the product of a slowly-varying envelope function, denoted as (), and a fast carrier signal, exp(−).Thus, the electric field of the incident wavepacket can be expressed as  8) approximately becomes −⟨  (, )|Δ()| ̃⟩.In the next Section, we will show that this approximation is only valid for on-resonance pulses whose central frequency (or carrier frequency) matches the QNM frequency ( ≈  ̃).

Local excitation approximation.
Another commonly encountered scenario involves resonators characterized by subwavelength dimensions, such as RLC lump-element circuits.In these instances, it is reasonable to assume that the incident wavepacket does not vary at the scale of the resonator, and the overlap integrals ⟨  |Δ()| ̃⟩ reduces to ⟨  (  , )|Δ(  )| ̃⟩,   being the 'center of coupling' of the resonator.This approximation holds true even for significantly larger systems, including ring resonators or photonic crystal cavities, which are typically much larger than the wavelengths of the confined radiation.In such cases, the input port typically consists of a waveguide, and the overlap integral ⟨  (  , )|Δ(  )| ̃⟩ remains predominantly localized.This is primarily due to the fact that the evanescent tail of the waveguide mode generally only intersects a small portion of the cavity.Typically, the overlap extends over a volume no greater than a fraction of  3 .
By utilizing these two approximations, we effectively achieve a primary goal of this study: deriving an analytical expression for the coupling constant outlined in Eq. ( 1) where  ̃ represents the electric field of the normalized QNM and  is the carrier frequency.
The prevalence of slowly-varying envelopes and local excitations in many applications and the absence of a QNM framework in the recent past may explain why earlier experiments did not reveal any discrepancy in CMT predictions.
We now provide a qualitative assessment of the impact of these two approximations.For the simulation, we consider a parallelepiped resonator characterized by an isotropic spatially uniform Δ (see Fig. 2).The resonator is assumed to be illuminated by a pulse with an electric field  0 () exp(−), where  =  − ( −   )/,  = 1,  = 1 − 0.05 is a complex frequency and () = (1 +  − ) −1 is the sigmoid function, approaching a Heaviside function for large values of .Note that the carrier frequency  nearly matches the resonance frequency  ̃ = 1.1 − 0.0025 .The pulse is assumed to hint the 'center of coupling'   of the resonator at time  = 0.
To simplify our qualitative analysis, we also assume that the electric field of the QNM is spatially uniform in the volume   = .Thus, the overlap integral in Eq. ( 8) simplifies to  Δ  ̃ •  0 〈() exp(−)〉  .Here, the parameters  , Δ ,  ̃ and  0 are scaling parameters and their values are inconsequential for the comparison.What is crucial for assessing finite size effects is the resonator length, , along the direction of propagation of the driving pulse (see Fig. 2).Complex frequency illuminations have recently attracted significant interest for a range of intriguing applications, such as loss compensation [23].In our investigation, we examine them to conveniently adjust the sharpness of the pulse envelope.Two values of  are considered:  = 0.1 (gentle envelop, Fig. 3a) in and  = 1 (steep envelop, Figs.3b and 3c).
Figure 3 summarizes the results of a comparison of the predictions obtained with the EME equation (cyan curve) and Eq. ( 10) (red curve), for increasing resonator sizes () and two pulse shapes () shown in grey.
We start with a subwavelength resonator ( = 0.1) and a gentle wavefront ( = 0.1) that resembles a Gaussian pulse.Consistently with other studies (not reported for compactness) with Gaussian pulses, we find a perfect agreement, see Fig. 3a.
To further test the predictive force of Eq. ( 10), we keep  = 0.1 and consider a steep envelop ( = 1).Again, we found a remarkable agreement with numerical data obtained directly by solving Eq. ( 8).This result is not reported in Fig. 3. Then we increase the resonator size ( = 0.5).We observe a discrepancy.However, if we multiply the prediction of Eq. ( 10) by 0.4, the rescaled data aligns perfectly with the prediction of the EME equation, see Fig. 3b.This is a noteworthy result.The CMT with a fitted coupling constant remains very accurate even when the resonator size is half a wavelength and fast-varying envelopes.
We finally consider a resonator size of one wavelength ( = 1).In this case, even with a large rescaling parameter, the CMT cannot predict the double peak response obtained with the EME equation (Fig. 3c).This implies that, even with fitted coupling coefficients, the CMT equation cannot be accurate and the EME formula emerges as especially important for investigating the fast dynamics of large-sized optical cavities involving linear or nonlinear effects [24].
Illumination by pure harmonic waves.Before we draw our conclusions, let us mention the special case where the resonator is illuminated by a harmonic oscillation wave, This expression remains accurate regardless of the size of the resonator.It explains why the CMT has been successfully employed with monochromatic illumination for large resonators, such as racetrack type resonators, for which the coupling sections may span several tenths of wavelengths.

Fig. 2 .
Fig. 2. Parameters used for the simulations reported in Fig. 3 for various pulse shapes () and resonator dimensions () characterized by the parameter  along the direction  of the incident plane wave  0 () exp(−).