Issue |
J. Eur. Opt. Society-Rapid Publ.
Volume 21, Number 1, 2025
EOSAM 2024
|
|
---|---|---|
Article Number | 5 | |
Number of page(s) | 7 | |
DOI | https://doi.org/10.1051/jeos/2024048 | |
Published online | 24 January 2025 |
Research Article
On the exact Maxwell evolution equation of resonator dynamics
1
Laboratoire Photonique, Numérique et Nanosciences (LP2N), IOGS – Université de Bordeaux-CNRS, 33400 Talence cedex, France
2
Aix-Marseille Université, Laboratoire ADEF, Campus Universitaire de Saint-Jérôme, 52 Avenue Escadrille Normandie Niémen, 13013 Marseille, France
3
CPT, Aix-Marseille Université, Université de Toulon, 13288 Marseille, France
* Corresponding author: philippe.lalanne@institutoptique.fr
Received:
29
October
2024
Accepted:
10
December
2024
In a recent publication [Opt. Express 32, 20904 (2024)], the accuracy of the main evolution equation that governs resonator dynamics in the coupled-mode theory (CMT) was questioned. The study concluded that the driving force is proportional to the temporal derivative of the excitation field rather than the excitation field itself. This conclusion was reached with a derivation of an “exact” Maxwell evolution (EME) equation obtained directly from Maxwell’s equations, which was further supported by extensive numerical tests. Hereafter, we argue that the original derivation lacks mathematical rigor. We present a direct and rigorous derivation that establishes a solid mathematical foundation for the EME equation. This new approach clarifies the origin of the temporal derivative in the excitation term of CMT and elucidates the approximations present in the classical CMT evolution equation through a straightforward argument.
Key words: Coupled-mode theory / Electromagnetic resonance / Quasinormal mode / Resonant scattering
© The Author(s), published by EDP Sciences, 2025
This is an Open Access article distributed under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.
1 Introduction
Many concepts in physics are based on normal modes of conservative systems, such as molecular orbitals and excitation energies. In these self-adjoint problems, the response of the system can be represented as a sum over the complete set of normal modes.
However, when energy dissipation occurs through processes like absorption or leakage in open systems, the system ceases to be conservative. The corresponding operator becomes non-self-adjoint and its spectrum becomes continuous. The spectrum also includes an infinite set of quasinormal modes (QNMs) characterized by complex eigenfrequencies (or eigenenergies), which produce characteristic damped oscillations in the system’s temporal response. A significant challenge across various disciplines [1–4] focuses on using QNM expansions to represent this response, akin to the methods applied to closed systems with normal modes.
In electromagnetism, this challenge was effectively addressed long ago with the temporal coupled-mode theory (CMT) [5, 6]. CMT offers a computationally effective, systematic, and intuitive framework for characterizing interactions between different modes within resonators and predicting their temporal dynamics. As a result, it has become an essential tool for designing and optimizing optical resonators, facilitating the development of innovative photonic technologies in both linear [3, 4, 7] and nonlinear optics [8].
Typically, CMT assumes the presence of m modes and n ports. In this context, a mode corresponds to a QNM, while a port represents a propagating channel where the QNMs may decay. It is perhaps worth noting that our understanding of ports may not fully reflect reality, as resonators are open systems that dissipate energy across a continuum of “ports”. A typical scenario involves a resonator on a layered substrate, illuminated by a plane wave (one of the system’s many ports). This wave scatters into a continuum of ports composed of all the plane waves in the substrate or superstate, as well as any potential guided modes of the layered substrate [9].
The CMT equations consist of two primary equations [5–8]. The first, known as the evolution equation, describes how the amplitudes of the modes evolve in response to incoming waves from the ports. This equation is crucial for understanding the dynamics of the resonator. The second equation coherently combines resonant-assisted and background pathways to predict the fields that couple out into the ports, enabling the calculation of reflection and transmission coefficients for the various ports.
Despite its extensive application and many successes, CMT remains a phenomenological theory for Fano-Feshbach resonances [10]; all coefficients describing the coupling between the resonator modes and the ports are fitted from experimental or numerical data. There are exceptions for certain geometries. For instance, for microring filters, the coupling issue may accurately reduce to calculating the interaction between guided modes in bent waveguides with a Hermitian theoretical framework [11].
In this work, we concentrate on the first equation and its theoretical derivation. In a recent study [10], by trying to establish a rigorous foundation for the CMT using electromagnetic QNM theory, an alternative evolution equation, distinct from the one proposed in the CMT, was derived directly from Maxwell equations. This new equation, termed “exact” Maxwell Evolution (EME) equation, shares some similarities with the classical CMT evolution equation but also shows notable differences. A key distinction is in the nature of the driving term: in the EME equation, the excitation is proportional to the temporal derivative of the incident field, whereas in the CMT evolution equation, it is proportional to the incident field itself. This unexpected result was validated through a numerical test, which remains unchallenged in this study. Additionally, it was supported by a mathematical demonstration.
Here, we highlight that the demonstration lacks from mathematical rigor and propose a new demonstration that is both direct and rigorous. Importantly, the new demonstration also helps clarifying the key difference between a driving force proportional to the incident field and one proportional to its temporal derivative.
The manuscript is organized as follows. Section 2 helps the reader in contextualizing the topic by reintroducing the origin of the double integral necessary for deriving the EME equation. Section 3 highlights several aspects of the demonstration in Section 2 of [10] which lack mathematical rigor. Section 4 provides a direct and explicit demonstration of the EME equation, which is further supported in Appendix A by a more general demonstration with minimum assumptions on the incident field. Section 5 clarifies the approximation carried out by assuming a driving force directly proportional to the incident field. Section 6 concludes the work.
2 Background
We start by considering recent results obtained for QNM-expansions in the spectral domain for harmonic fields. Electromagnetic QNMs, labeled by the integer m, are source-free solutions of Maxwell equations, ,
, which satisfy the outgoing-wave condition for |r|→∞ [3, 4]. Here,
and
represent the QNM electric and magnetic fields, respectively, ε denotes the possibly dispersive permittivity tensors. The QNM fields decay exponentially over time, so that
in the exp(−iωt) notation. From this point forward, we assume the QNMs are normalized, which is a crucial step in any QNM theory. There are several methods available for normalizing electromagnetic QNMs, all of which are thoroughly discussed in [12]. In this work, we utilize the so-called PML normalization approach, where PML stands for Perfectly Matched Layer.
In the scattered-field formulation, the total field at frequency ω is decomposed as a sum of the background field [Eb(r, ω), Hb(r, ω)]exp(−iωt) and the field [Es(r, ω), Hs(r, ω)]exp(−iωt) scattered by the nanoresonator. The scattered field can be expanded in an infinite set of normalized QNMs and PML modes [12, 13](1)
where αm(ω) denotes the excitation coefficient of the mth modes. Note that the normalized modes satisfy [12]. The scattered and normalized mode fields thus have different units.
The QNM framework based on a complex-mapping regularization provides a unique expression for the excitation coefficient in equation (1) [13](2)
for the case of nondispersive materials under consideration here. An equivalent EME equation for dispersive systems with Drude-Lorentz permittivities is derived in [13]. In this case, the driving force consists of two components: one proportional to the temporal derivative of the excitation field, and the other to the excitation field itself. Numerical tests conducted on a silver bowtie antenna, as shown in Figure 2b of [13], strongly validate the robustness of the analysis. However, it would be valuable to expand the study further, for example, by conducting additional numerical tests for different geometries. It would be also valuable to generalize the theory to materials with arbitrary dispersion.
Equation (2) corresponds to equation (3) in [10]. αm(ω) essentially represents an overlap integral between the QNM and the background field. Δε specifies the permittivity variation used for the scattered-field formulation, given by Δε(r) = εR(r)−εb(r), where εR(r) and εb(r) denote the permittivities of the resonator and the background, respectively. In general, Δε(r) is non-zero within a compact volume VR(r) that defines the resonator in the scattered-field formulation.
We now consider that the background field is an optical pulse, Eb(r, t), i.e. a wave packet that can be Fourier transformed . Driven by the incident pulse, the resonator scatters a time-dependent electric field, Es(r, t). Every infinitesimal frequency component of the background field Eb(r, ω)dω gives rise to an infinitesimal harmonic scattered field
, and the scattered field in the time domain is obtained by summing up all the frequency components,
. The latter is conveniently expressed with a QNM expansion [3, 13]
(3)with
(4)
Replacing the spectral component by the Fourier transform of the wavepacket Eb(r, t) (please note that Fourier transforms are represented with a horizontal bar), we obtain
(5)
Following [10], we introduce the overlap . Equation (5) becomes
(6)which coincides with equation (5) in [10].
In [10], it is shown that(7)
However, while we will see that this equality applies to a broad spectrum of complex functions Om(t′), the proof is deficient in rigor, as we delineate in the subsequent Section.
3 Critical assessment of the proof in [10]
Complex analysis and generalized functions are used to derive equation (7) in [10]. While the equality holds, we have identified that several steps lack mathematical rigor in our scrutiny of the proof.
3.1 Unjustified use of Fubini theorem
The expression in equation (6) contains two integrals that have to be calculated sequentially, first considering the integral on t′ and then on ω. In [10], the order of integration is inverted, implicitly using the Fubini theorem, whose assumptions are not satisfied. Indeed, one can first observe that an application of the Fubini-Tonelli theorem to the absolute value of the integrand shows that the latter is not integrable over (−∞, ∞)2 although Om(t) is assumed to be in , the space of functions integrable over the real line. Second, it can be checked that the integral
diverges (see below), which implies that the Fubini theorem does not apply here.
3.2 Divergent integrals
We further revisit the computation of the integral with respect to ω in [10], where it is “shown” that for t ≠ t′(8)where H stands for the standard Heaviside function.
To derive equation (7), in [10], the residue theorem is applied twice on two semicircles in the upper and lower half of the complex plane (we refer here to the calculation of terms A and C in equation (6) in [10] with the red and green arcs shown in Fig. 1 in [10]). We believe that the derivation is incorrect.
![]() |
Figure 1 Comparison of the EME and approximate equations for a 1D Fabry-Perot resonator illuminated by a short 10-fs pulse. (a) Sketch of the geometry. (b) QNMs of the resonator in the complex frequency plane, calculated using the MAN software. The circle size indicates the excitation level of the QNMs. The real part of the eigenfrequency for QNM 1 corresponds to the central frequency of the incident pulse. (c) The values of β1 and β2 obtained from both the EME and approximate equations. (d) A comparison of the total-field reconstruction using the EME and approximate equations with 250 QNMs against full-wave numerical results from COMSOL. |
Our primary concern lies in the fact that, for ω running over the green or red arcs, it is not consistently true that Im(ω)→∞ (consider instances where ω is close to or at a specific distance from the blue diameter/line in Fig. 1). Thus, the contributions of the green and red integrals over the arcs cannot be assumed to be null.
Actually, diverges for any real t, t′. To show that, we use the equality
and obtain
(9)
The second integral on the right-hand side is defined and is equal to(10)
The derivation of equation (10) can be accomplished by employing a methodology akin to that outlined in [10] utilizing the residue theorem. By employing the same integration contours depicted in Figure 1 in [10], this time, the integrand tends to 0 as the radius of each semicircle goes to ∞, simply because (which is not the case for
).
The comparison between equations (8) and (10) is eloquent. Furthermore, we observe that the first integral on the right-hand side of equation (9), , is divergent (or not defined) for any real t, t′. As the second integral in this equation is defined, see equation (10), the integral on the left-hand side,
, is also divergent.
4 Demonstration of equation (7)
In this Section, we provide a demonstration of equation (7) that does not rely on complex analysis. By introducing the spectral overlap , equation (2) becomes
(11)
We take the Fourier transform of equation (11) to move in the temporal domain(12)
To obtain equation (12), we need two assumptions on the driving pulse
-
• the Lebesgue integral of ωαm(ω) is finite, i.e.

-
• the Lebesgue integral of ωOm(ω) is finite, i.e.

Let us consider the left-hand term: . It is equal to
. We can do exactly the same for the last term of equation (12),
. Then, directly from equation (12), we have
(14)which corresponds to the main result obtained for non-dispersive materials in [10], specifically the EME evolution. It is straightforward to verify that the expression for βm(t) provided in equation (7) satisfies equation (14). This completes the proof of equation (7).
In Appendix A, we present an alternative demonstration based on less stringent (possibly minimal) assumptions regarding the driving pulse: and
.
5 Clarification of the origin of the temporal derivative in the driving force
This section emphasizes the approximation made in classical CMT, where the driving force is assumed to be proportional to the incident field rather than its time derivative.
To clarify the approximation introduced in the CMT, we remember that a derivative in the real temporal domain amounts to multiply by the conjugated variable in the frequency domain in Fourier (or Laplace) analysis. We thus intuitively multiply the expression of the excitation coefficient by to remove the ω-dependence in the numerator of equation (2). We obtain an approximate expression of the excitation coefficient denoted by
(15)
We expect the approximate expression to be most accurate when the frequency of the monochromatic incident field is close to the QNM resonance frequency. As demonstrated in Section 4, the temporal excitation coefficient is given by . This leads by differentiation to
(16)in contrast to the EME equation
(17)
Equation (16) aligns with the usual approach in CMT, here the driving force is proportional to an overlap integral between the incident field (rather than its derivative) and the QNM field, with a prefactor of as expected from dimensional analysis.
Next, we evaluate the accuracy of equation (16) for a 1D Fabry-Perot resonator in air. The resonator is illuminated at normal incidence by a 10-fs plane-wave Gaussian pulse (Fig. 1a). The central frequency of the pulse corresponds to the real part of the eigenfrequency of one of the QNMs (QNM 1 in Fig. 1b).
Figure 1c compares |β1| and |β2|, computed using the approximate and EME equations. For each QNM, Om(t) is obtained by simply computing the overlap between the QNM electric field and the Gaussian pulse Eb(r, t):
. The values of |β1| and |β2| are then determined by solving the differential equations (16) and (17). Significant differences are observed primarily for the off-resonant QNM 2, for which the resonance frequency differs from the pulse central frequency – the approximation made in equation (15) by multiplying by the excitation coefficient by
becomes inaccurate. The reverse occurs when we adjust the central frequency of the incident pulse to match
, implying the QNM field distribution does not play any role in the difference reported. Additionally, note that in classical CMT, the coupling coefficient is typically fitted, and more precise predictions than those derived from equation (17) are possible, as discussed in Section 3 of [10].
Temporal excitation coefficients, being abstract quantities, cannot be directly measured. Figure 1d offers a comparison using measurable quantities, such as the total field. The two panels in Figure 1d show the reconstructed total field: first, inside the resonator for a fixed time t = 65 fs, and second, at a fixed coordinate x = 150 nm as a function of time. For the reconstruction, 250 QNMs are used in the expansion of equation (1). The EME equation (17) yields a reconstruction that perfectly matches reference data from COMSOL time-domain simulations, while the approximate equation (16) exhibits significant deviations.
It is important to note that the driving pulse is ultra-short and has a broad spectral range. Furthermore, all the QNMs have low quality factors, making them easily excited even when the central frequency of the driving pulse deviates significantly from the QNM resonance frequency. This is why we used 250 QNMs in the reconstruction shown in Figure 1d. With such a large number of QNMs, the maximum difference between the black curves of equation (17) and the COMSOL data is less than 0.005 for both curves in Figure 1d, providing a strong quantitative validation of the approach.
6 Conclusion
The derivations of the EME equation in Section 4 and Appendix A are straightforward and do not depend on complex analysis or generalized functions, providing a solid foundation for the equation across a broad range of incident pulses.
The numerical tests in Section 5 are intentionally carried out on a simple 1D geometry, commonly used with CMT, which emphasizes the validity of the EME equation. This is achieved by using a driving term instead of
when Δε(r) is non-dispersive. As noted in [10], under the slowly-varying envelope approximation,
can often be replaced by −iωEb(r, t), where ω is the carrier frequency of the driving pulse. This suggests that fitting a coupling constant with a driving term proportional to Eb(r, t) will yield accurate predictions when
, as demonstrated in Section 5.
Nonetheless, many modern electromagnetic software tools, including general solvers like COMSOL and specialized options [14–16] now compute QNMs. Additionally, QNM normalization has been mastered and is implemented in specialized software [16]. As a result, the time-derivative-based analytical expression, , can be efficiently calculated as a spatial overlap integral. Therefore, we anticipate that CMT models without fitting parameters will be readily used to interpret both experimental and numerical results. Notably, the temporal toolbox in [16] effectively addresses the exact Maxwell evolution equation, and efforts are already underway to model the second CMT equation for port coupling [13, 16–18].
Funding
PL acknowledges financial supports from the WHEEL (ANR-22CE24-0012-03) Project. Rachid Zarouf is partially supported by the pilot centre for research in education Ampiric, funded by the France 2030 Investment Program, as part of the action “Territories of Educational Innovation” operated by the Caisse des Dépôts.
Conflicts of interest
The authors declare no conflict of interest.
Data availability statement
Data underlying the results presented in this paper are not publicly available at this time but may be obtained from the authors upon reasonable request.
Author contribution statement
All authors contributed equally. TW provided the computational results of the figure.
References
- Zworski M, Mathematical study of scattering resonances, Bull. Math. Sci. 7, 1–85 (2017). https://doi.org/10.1007/s13373-017-0099-4. [CrossRef] [Google Scholar]
- Nollert H-P, Quasinormal modes: the characteristic sound of black holes and neutron stars, Class. Quantum Grav. 16, R159–R216 (1999). https://doi.org/10.1088/0264-9381/16/12/201. [NASA ADS] [CrossRef] [Google Scholar]
- Lalanne P, Yan W, Kevin V, Sauvan C, Hugonin J-P, Light interaction with photonic and plasmonic resonances, Laser Photon. Rev. 12, 1700113 (2018). https://doi.org/10.1002/lpor.201700113. [Google Scholar]
- Benisty H, Greffet JJ, Lalanne P, Introduction to Nanophotonics. Chapter 6 (Oxford University Press, New York, 2022). [Google Scholar]
- Haus HA, Waves and Fields in Optoelectronics, 1st ed. (Prentice-Hall, Englewood Cliffs, 1984). [Google Scholar]
- Shu W, Wang Z, Fan S, Temporal coupled-mode theory and the presence of non-orthogonal modes in lossless multimode cavities, IEEE J. Quantum Electron. 40, 1511–1518 (2004). https://doi.org/10.1109/JQE.2004.834773. [NASA ADS] [CrossRef] [Google Scholar]
- Zhang H, Miller OD, Quasinormal Coupled Mode Theory. arXiv preprint. arXiv:2010.08650. https://doi.org/10.48550/arXiv.2010.08650. [Google Scholar]
- Christopoulos T, Tsilipakos O, Kriezis E, Temporal coupled-mode theory in nonlinear resonant photonics: from basic principles to contemporary systems with 2D materials, dispersion, loss, and gain, J. Appl. Phys. 136, 011101 (2024). https://doi.org/10.1063/5.0190631. [Google Scholar]
- Yang J, Hugonin JP, Lalanne P, Near-to-far field transformations for radiative and guided waves, ACS Photon. 3, 395–402 (2016). https://doi.org/10.1021/acsphotonics.5b00559. [Google Scholar]
- Wu T, Lalanne P, Exact Maxwell evolution equation of resonators dynamics: temporal coupled-mode theory revisited, Opt. Express 32, 20904–20914 (2024). https://doi.org/10.1364/OE.517237. [CrossRef] [Google Scholar]
- Popović MA, Manolatou C, Watts MR, Coupling-induced resonance frequency shifts in coupled dielectric multi-cavity filters, Opt. Express 14, 1208–1222 (2006). https://doi.org/10.1364/OE.14.001208. [NASA ADS] [CrossRef] [Google Scholar]
- Sauvan C, Wu T, Zarouf R, Muljarov EA, Lalanne P, Normalization, orthogonality, and completeness of quasinormal modes of open systems: the case of electromagnetism, Opt. Express 30, 6846–6885 (2022). https://doi.org/10.1364/OE.443656. [NASA ADS] [CrossRef] [Google Scholar]
- Yan W, Faggiani R, Lalanne P, Rigorous modal analysis of plasmonic nanoresonators, Phys. Rev. B 97, 205422 (2018). https://doi.org/10.1103/PhysRevB.97.205422. [CrossRef] [Google Scholar]
- Lalanne P, Yan W, Gras A, Sauvan C, Hugonin J-P, Besbes M, Demesy G, Truong MD, Gralak B, Zolla F, Nicolet A, Binkowski F, Zschiedrich L, Burger S, Zimmerling J, Remis R, Urbach P, Liu HT, Weiss T, Quasinormal mode solvers for resonators with dispersive materials, J. Opt. Soc. Am. A 36, 686–704 (2019). https://doi.org/10.1364/JOSAA.36.000686. [Google Scholar]
- Betz F, Binkowski F, Burger S, RPExpand: software for Riesz projection expansion of resonance phenomena, SoftwareX 15, 100763 (2021). https://doi.org/10.1016/j.softx.2021.100763. [CrossRef] [Google Scholar]
- Wu T, Arrivault D, Yan W, Lalanne P, Modal analysis of electromagnetic resonators: user guide for the MAN program, Comput. Phys. Commun. 284, 108627 (2023). https://doi.org/10.1016/j.cpc.2022.108627. [CrossRef] [Google Scholar]
- Betz F, Binkowski F, Hammerschmidt M, Zschiedrich L, Burger S, Resonance expansion of quadratic quantities with regularized quasinormal modes, Phys. Status Solidi A 220, 2200892 (2023). https://doi.org/10.1002/pssa.202200892. [NASA ADS] [CrossRef] [Google Scholar]
- Zolla F, Nicolet A, Demésy G, Photonics in highly dispersive media: the exact modal expansion, Opt. Lett. 43, 5813–5816 (2018). https://doi.org/10.1364/OL.43.005813. [NASA ADS] [CrossRef] [Google Scholar]
Appendix A
Demonstration of equation (7) via Fourier inversion
In this Appendix, we again aim at calculating (see Eq. 7)(A1)with minimum assumption on the incident field, i.e. on Om(t′). The approach just assumes that Om is in
, the space of functions integrable over the real line. It is based on the well-known Fourier inversion formula.
For fixed values of ω, we first calculate the integral , which is equal to
. Further integrating over
with respect to ω, we aim at calculating
. Replacing
by
, we first obtain
and thus find
(A2)
To evaluate the integral on the right-hand side, we demonstrate that the function G(t) = , where
satisfies a simple differential equation that we can solve. It can be verified that G is a differentiable function of the variable
, because the conditions of the well-known differentiation theorem for integrals are met for functions in
. Differentiating G with respect to t yields
. Considering that
(the limit and the integral can be interchanged due to the integrand being suitably dominated), we obtain
.
Thus and coming back to equation (A2), we find
(A3)
This concludes a demonstration of equation (7).
All Figures
![]() |
Figure 1 Comparison of the EME and approximate equations for a 1D Fabry-Perot resonator illuminated by a short 10-fs pulse. (a) Sketch of the geometry. (b) QNMs of the resonator in the complex frequency plane, calculated using the MAN software. The circle size indicates the excitation level of the QNMs. The real part of the eigenfrequency for QNM 1 corresponds to the central frequency of the incident pulse. (c) The values of β1 and β2 obtained from both the EME and approximate equations. (d) A comparison of the total-field reconstruction using the EME and approximate equations with 250 QNMs against full-wave numerical results from COMSOL. |
In the text |
Current usage metrics show cumulative count of Article Views (full-text article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.
Data correspond to usage on the plateform after 2015. The current usage metrics is available 48-96 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.