Nonreciprocal Cavities and the Time-Bandwidth Limit

The time-bandwidth limit inherently relates the lifetime of a resonance and its spectral bandwidth, with direct implications on the maximum storage time of a pulse versus its frequency content. It has been recently argued that nonreciprocal cavities may overcome this constraint, by breaking the strict equality of their incoupling and outcoupling coefficients. Here, we generally study the implications of nonreciprocity on resonant cavities and derive general relations, stemming from microscopic reversibility, that govern their dynamics. We show that nonreciprocal cavities do not provide specific advantages in terms of the time-bandwidth limit, but they may have other attractive properties for nanophotonic systems.

The demands of integrated photonic circuits have motivated a large body of research with the objective of scaling down and integrating crucial optical components for information processing [1][2][3][4][5][6][7][8][9][10][11][12][13]. A common approach to enable small footprints is to use resonant cavities [6][7][8][9][10][11][12][13], which can dramatically enhance light-matter interactions by storing light over the resonance lifetime. This reduced footprint, however, comes at the cost of operational bandwidth: a single resonant cavity adheres to a strict correspondence between its resonance bandwidth ∆ω and its lifetime ∆t [14]: ∆t∆ω = 2. As a relevant example in nanophotonics, resonators are frequently used to impart delays on optical pulses, given their ability to temporarily store light. The strict relation between lifetime and bandwidth of a cavity implies a trade-off between the spectral bandwidth of the pulse that can be stored and the temporal delay that can be imparted on it. Larger structures (e.g., arrangements of multiple cavities [15]) may partially relax this bound, but still follow an analogous trade-off between the delay-bandwidth product and the overall footprint [16]. Recent approaches to overcome this limitation have relied on time-varying systems [17][18][19] or on nonlinearities [20], but efficiently implementing these schemes often becomes technologically challenging, or imposes limitations in their operation.
A thought-provoking proposal to alleviate the strict correspondence between bandwidth and lifetime of linear, time-invariant structures has been recently been put forward, based on breaking reciprocity with a static magnetic bias [21]. The authors argue that the supported bandwidth ∆ω scales with the rate at which energy enters the system, while the lifetime is inversely proportional to the rate at which energy exits it. These rates are known to be equal in reciprocal cavities [22], resulting in the mentioned time-bandwidth relation, but Ref. [21] suggests that nonreciprocal cavities may support different input and output rates, and thus can surpass this limit. While this idea would have a large impact on many photonic applications relying on resonances, it raises concerns regarding its thermodynamic validity [23]: unequal input and output rates violate energy conservation in a passive system.
Inspired by this proposal, in the following we develop a general temporal coupled-mode theory broadly describing the dynamics of nonreciprocal cavities, and show that incoming and outgoing rates are always strictly related through a number of identities that generalize the known relations in reciprocal systems. For example, we prove that in nonreciprocal cavities the total rates must always be equal in magnitude. We also show that timebandwidth product is strictly limited in any linear, timeinvariant cavity, independent of whether it is reciprocal or not, and provide numerical demonstrations. We also derive bounds for the coupling coefficients of nonreciprocal cavities, and design cavities operating at these bounds.
We start our discussion by studying a general resonant cavity with complex mode amplitude a connected to n ports, with incoming complex amplitudes s + = (s +,1 , s +,2 , · · · , s +,n ) . Such a system is described by the equation of motion [22,24] where k = (k 1 , k 2 , · · · , k n ) is the vector of input coupling coefficients, ω 0 is the eigenfrequency, γ r is the radiation loss rate into the ports, and γ i is the intrinsic loss rate due to absorption inside the cavity. The amplitudes are normalized such that |a| 2 is the stored energy in the cavity and |s +,i | 2 is the incoming power at port i. The amplitude of the reflected modes is given by Here, C is the direct reflection coefficient, which is unitary in the absence of loss, and d is the output coupling coefficient from the cavity mode into the ports. Power conservation ensures that 2γ r = d † d [24]. In addition, by assuming time-harmonic excitation (s + = s +,0 exp (iωt)) we find from Eq. 1: Eqs. 1-3 and the relation 2γ r = d † d apply to both reciprocal and nonreciprocal cavities. In addition, reciprocity dictates that d = k and Cd * = −d (asterisk denotes complex conjugate) [24], which follow from time-reversal symmetry. However, if the cavity is biased with a quantity odd symmetric under time-reversal, these two identities do not apply, as the time-reversed system is no longer equivalent to the original one. In the following, we therefore denote the coefficients in the time-reversed scenario with a tilde, asC,d andk.
Here the superscripts and † denote the regular and Hermitian transpose. These relationships are markedly different from the regular expressions for reciprocal systems [24,25], particularly because they involve the original system and its time-reversed counterpart (which, as we will show, can differ strongly in its response). This is a consequence of microscopic reversibility [26]: while the system itself is not time-reversal symmetric, global time-reversal invariance still applies, and the original and time-reversed systems are necessarily related.
Eq. 4f is the fluctuation-dissipation relation, relating dissipation (d) to how the cavity responds to external fluctuations (k) [14, 23,27]. Together with Eq. 4e, this relation is arguably the most important result for the following discussions, as they put stringent restrictions on the achievable asymmetry in the coupling coefficients: while it is possible to achieve d i = k i at any individual port, the positive-definite norm of k and d has to be equal (see Eq. 4f). Finally, we note that Eqs. 4 readily reduce to the known identities for reciprocal cavities [24], if we assume that the system is time-reversal symmetric, i.e.d = d,k = k, and C = C . Eqs. (4) have important general consequences in the context of nonreciprocal cavities [28,29]. As a relevant example, consider a nonreciprocal waveguide supporting a unidirectional mode, i.e., no backward channel into which energy can be reflected, and terminated into a cavity, similar to geometries studied in [30,31] and shown schematically in Fig. 1a. For the nonreciprocal waveguide, we use the interface between a magnetized magneto-optic semiconductor and a dielectric material [28,[32][33][34], which for the combination InSb and Si is known to support a unidirectional surface wave around 1.5 THz [21,35,36] (see all geometry details in [14]). The waveguide is terminated onto a resonant lossless rectangular cavity (20 by 30 m in size), enclosed by perfect electric conducting (PEC) walls, connected through a small aperture (0.5 µm).
We use a home-built FDTD code [14] to excite the unidirectional waveguide with a broadband pulse (f 0 =1.5 THz, ∆f =0.16 THz), and record the fields inside the cavity as a function of time. Fig. 1b shows the Fourier transform of the electric field induced at the center of the cavity (blue line), showing a sharp Lorentzian resonance at 1.52 THz, with linewidth 2γ = 45.2×10 9 rad/s. When exciting the cavity from the interior, we find exactly the same decay rate γ = 22.6 × 10 9 rad/s by examining the field decay (inset of Fig. 1b). This should not be surprising, as Eq. (3) applies independent of whether the system is reciprocal or not, and ensures that the bandwidth of any linear cavity is solely determined by the decay rate γ = γ r + γ i , i.e., by the internal loss and the outcoupling coefficients d † d/2 = γ r . This is supported by the observation that the fields at the termination outside the cavity have a much larger bandwidth (Fig. 1b, red line), equal to the bandwidth of the incident pulse, but that these fields cannot be forced into the cavity. In other words, as a first important consequence of the previous theory, the time-bandwidth product of any linear cavity satisfies ∆t∆ω = 2, and it is controlled only by the output coupling coefficient d. The incoupling coefficient k, according to Eq. (3), only controls the amount of stored energy in the resonator.
More surprising is that the cavity has a finite decay rate, despite the fact that it is both lossless (γ i = 0) and connected to a waveguide that does not support backward modes (γ r = 0). While the decay process is not immediately apparent, it is consistent with Eq. 4f, which requires that if the cavity can be excited it must also be able to decay to maintain equilibrium. This paradox can be resolved by inspecting the electric field distribution during decay (without an external excitation), shown in Fig. 1c. We notice a strong hot spot at the wedge formed by the Si/InSb waveguide and the PEC termination which, given the finite material loss in InSb, dissipates all incoming energy and thus sustains the cavity decay. We also observe this hotspot when exciting from the waveguide port, both on-and off-resonance (Figs. 1d,e) [21,[36][37][38][39]. Since both the cavity and the input port can excite this wedge mode, they interfere in the process of dissipation (as evidenced in the bump through the resonance in the red curve in Fig. 1b). As a result, the wedge mode cannot be treated as a simple internal loss process, as γ i , but needs to be treated as an additional output port [14]. Then, writing k = (k r , k w ) and d = (d r , d w ) for the input and output coefficients for the radiation from/into the waveguide and wedge, respectively, we find that d r = 0 = k r is permitted, while k † k = d † d is simultaneously satisfied. In other words, this means that while the input and output coefficients may differ at each individual port, but the total input and output rates must always be equal. It follows that in the special case of a one-port system, like the one generally described in [21], the input and output coefficients must necessarily be equal in magnitude, but they may still differ in phase.
One may wonder what happens in the limit in which material loss in the unidirectional waveguide is zero, for which the wedge mode is expected to become nondissipative. This problem has been extensively discussed in the literature (see [37,[40][41][42][43] for a selection), and has been referred to as the "thermodynamic paradox", since it was believed to produce an inconsistency between Maxwell's equations and thermodynamics. The paradox was resolved by Ishimaru [36], who pointed out that the field distribution at the termination of a unidirectional waveguide is non-integrable if assumed lossless, and it sustains finite absorption even in the limit of vanishing material loss. A lossless terminated one-way waveguide is actually an ill-posed boundary-value problem, as Maxwell's equations do not guarantee a unique solution in the ideal lossless scenario [44]. The correct solution is found in the limit of vanishing loss, for which the dissipation rate remains finite. Aside from the wedge mode in this terminated unidirectional waveguide [21,[36][37][38][39], similar singularities that support finite absorption in the limit of vanishing material loss can be found in other extreme electromagnetic systems [45,46], and are a reminder that ideal lossless scenarios should be considered an artifact in electromagnetics and may lead to singularities or non-unique solutions.
While Eqs. 1-4 strictly prove that the time-bandwidth product of a linear time-invariant cavity is always equal to 2, and that therefore it is impossible to force broadband fields into a long-lived resonance, independent of reciprocity, Refs. [21,[36][37][38][39] do demonstrate broadband focusing of photons in an ultra-small volume near the termination, whose decay rate is unrelated to, and can be much slower than, the excitation. Consistently, our simulations confirm focusing of incident photons at the wedge (Figs. 1d,e) over a wide range of frequencies (Fig. 1b, red line), irrespective of the properties of the cavity resonance. We stress however that this broadband focusing is not directly a consequence of nonreciprocity: adiabatically tapered terminated plasmonic waveguides [47][48][49], which slowly focus the incoming fields towards an apex, perform the same function. The benefit of applying nonreciprocity is that the termination is automatically matched due to the absence of a backwards mode in the waveguide [30,31], relaxing the need for a carefully designed adiabatic transition that minimizes reflections.
After having demonstrated that nonreciprocity is not beneficial to break the trade-off between lifetime and bandwidth in linear, time-invariant cavities, in the following we will explore to what degree the input and output coupling may be made different in nonreciprocal cavities, and what functionalities can be enabled by such asymmetry. Achieving asymmetry in the coupling coefficients is important for e.g. circulators [50] and unidirectional heat transfer [51,52]. In the previous example, the absence of a backwards mode ensured d r = 0. However, we will now show that, even if there is a backwards mode, the input and output coefficients at a given port differ significantly in nonreciprocal cavities. To do so, we examine the same system as in Fig. 1, but now excited in the bidirectional regime at frequencies just below 1.25 THz. We therefore tune the cavity to 1.24 THz, by increasing the cavity size to 35.4 by 20 µm (see [14] for geometry details). With full-wave simulations we retrieve the complex cavity and reflection amplitudes, and through a fitting procedure we obtain k r = (2.65 + 0.308i) × 10 −4 rad/s and d r = (0.667 + 2.14i) × 10 −4 rad/s (see [14] for spectra and fits). As expected, now d r = 0 because of the presence of a propagating backward mode, but the two magnitudes are shown to be significantly different due to the presence of the wedge mode. When operating just below the unidirectional gap the system is still strongly asymmetric (in this case the forward and backward effective indices are 4.47 and 9.92 respectively), and it is to be expected that, as the asymmetry reduces, k r and d r will also approach each other. While the input and output coefficients at the input port may be made different, both in magnitude and phase, Eq. 4f requires the total rates of the cavity must be the same, and in the system under consideration this is guaranteed by the wedge mode, which balances out any asymmetry in k r and d r . Interestingly, this is, however, not the only requirement: the direct pathway places an additional bound on these coefficients through Eq. 4e. For a general two-port system with k = (k i , k j ) and d = (d i , d j ), we find using the Cauchy-Schwarz inequality: An example of this bound is shown in Fig. 2a as a function of |C ii | for |d i | 2 /d † d = 0.53 (|d i | 2 /d † d quantifies the fraction of radiated power flowing towards port i as the cavity decays). This specific example corresponds to the previously discussed bidirectional system, where the subscript i indicates radiation into/from the port (e.g. d i = d r ). The cross in the figure highlights the magnitude of k i and C ii for this geometry at 1.24 THz, falling within the bounds delimited by the solid line. The shape of the bound on k i strongly depends on d i and d j : Fig. 2b displays similar bounds for |d i | 2 /d † d = 0.1 and |d i | 2 /d † d = 0.9. In the system of Fig. 1 he absence of a backward mode dictates that both d i = 0 and C ii = 0, in which case Eq. 5 becomes |k i | = |d j |: as expected, power can only enter the cavity at the rate in which it is dissipated by the wedge mode or other output port. Similarly, if |C ii | = 1, which implies that the two ports are not directly connected, Eq. 5 requires that |k i | = |d i |, and thus input and output coefficients for each port can be nonreciprocal only in phase, as in a gyrator. This is a consequence of the fact that the ports only exchange power with the cavity itself, and hence cannot compensate for any asymmetry in its response.
Having discussed the general bounds on the input and output coupling coefficients of nonreciprocal cavities, we now investigate an extreme condition allowed by the bounds in Eq. 5, which may provide interesting functionalities for nanophotonic systems. We consider the scenario of a cavity connected to a bidirectional waveguide (|C ii | > 0) that cannot be excited from its input port (k i = 0), but that decays into it (|d i | > 0). To design such a cavity, it is sufficient to place its opening at a position with complete destructive interference between the forward and backward modes, as shown in Fig. 3a. Since nonreciprocal waveguides have different forward and backward field profiles [50], we need to ensure that the backward mode has higher fields at the top PEC plate, and that the reflection coefficient C ii compensates for differences in magnetic field amplitude. As shown in Fig. 3b, we achieve this condition by reverting the magnetic field bias with respect to the previous examples and optimizing of the reflectivity of the waveguide termination, using a dissipative surface impedance, to achieve |C ii | = 0.4. The necessity of a reduced reflection coefficient is consistent with our previous finding, in Eq. 5, that a nonreciprocal response in magnitude can only be achieved if |C ii | < 1.
We now find a locally vanishing magnetic field at the top PEC wall (Fig. 3a), implying that a cavity coupled to the waveguide at that location cannot be excited. In Fig. 3c, we place a cavity above the waveguide with its opening at the location of destructive interference, and a forward wave impinging from the channel indeed does not excite the cavity. If the waveguide were reciprocal, a zero in magnetic field at the cavity opening when excited from the port would imply that the cavity also cannot decay into the port. However, in the nonreciprocal case this is not so: if we excite the cavity from inside, we see it decay freely towards the output port (Fig. 3d). This is a result of the fact that, due to nonreciprocity and the different profiles of the forward and backward modes of the waveguide, for excitation from inside the cavity these modes are excited with different amplitudes at the aperture.
Interestingly, as indicated by Eqs. 4c,d we can swap the values of k and d by reversing the direction of the magnetic field bias (which is equivalent to a time-reversal operation). This is shown in Figs. 3e,f, which present the same structure but with a opposite magnetic bias. It can be seen that the cavity now can be efficiently excited from the port (Fig. 3e), but when the cavity is excited from inside it decays only towards the termination and get dissipated there, despite the presence of a backward propagating mode. It is interesting to point out that, in contrast to the example in Fig. 1, here d i = 0 while C ii > 0 -even though there is an available backwards mode, the cavity cannot couple to it. To conclude, in this paper we have presented a general theoretical framework describing the dynamics of nonreciprocal cavities. Our results show that, for single port systems, k and d can only differ in phase, and for k i and d i at any individual port to be different in magnitude at least one additional channel is necessary. This is consistent with the general principle that to realize a linear isolator it is necessary to have a loss channel [53,54] or an additional radiative port, as in a circulator. For systems with multiple ports we have shown that the sums of all input and output rates are necessarily equal: k † k = d † d. This implies that non-reciprocal cavities still follow strict bounds on their input and output coefficients, but nonetheless can be employed to realize highly non-trivial phenomena, such as cavities that can be pumped only one-way and release the energy into a totally different channel. We have shown that under a reversal of the magnetic field bias the functionalities of such a system strictly reverse. Finally, we have shown that the bandwidth of a linear, time-invariant cavity is always inversely proportional to its decay rate, both in reciprocal and nonreciprocal systems. The decay rate of a cavity is solely determined by the internal loss and the outcoupling coefficients d, not by the incoupling rate. Thus, it is not possible to force broadband fields into a long-lived resonance, independent of reciprocity, for which we have provided a numerical demonstration. Our results clarify claims that non-reciprocity may alleviate the strict limitations imposed by the trade-off between delay and bandwidth in photonic systems, and may help envisioning new efficient nonreciprocal components for, e.g., information processing, unidirectional transport, and thermal rectifiers. This work was supported by the Air Force Office of Scientific Research.

Details of Cavity Geometries and Materials 2
Anisotropic FDTD algorithm 3 Heuristic Derivation of CMT for Cavity and Nonreciprocal Waveguide 4 Using COMSOL for Nonreciprocal Waveguides 4 Fitting Simulations with the Coupled-Mode Theory Model 5

References 6
THE TIME-BANDWIDTH PRODUCT Starting from Eq. 3 in the main text, the stored energy of the resonance is given by: (S1) The full-width, half-maximum of this Lorentzian curve is given by 2γ. Similarly, from Eq. 1 in the main text, if follows that, if the cavity has a nonzero amplitude a 0 at t = 0 without the presence of an incoming wave, it evolves in time as Here, the lifetime of the resonance ∆t is given by the time it takes the resonance amplitude to reach e −1 of its original value: ∆t = γ −1 . Hence, taking the product of the lifetime and bandwidth, we find ∆ω∆t = 2.

PROOFS FOR NONRECIPROCAL CMT IDENTITIES
Here we provide proofs for the identities shown in the main text, which we repeat here for clarity: Before proving these relationships, it is important to point out the consequence of a time-reversal operation on general amplitudes of the mode and cavity amplitudes: performing a time-reversal operation leads to T : s + → s * − , T : s − → s * + , and T : a → a * , since it effectively conjugates the temporal exponent in each amplitude and reverses the direction of any vector (such as the propagation vector). In the following we provide the proofs not in order of Eq. 4 in the main text, but in an order that makes more sense with respect to interdependencies. These proofs are based on the assumption that the mode profile and frequency are not affected by the time-reversal operation, and for Eqs. 4b,c,e closely follow the derivations for the reciprocal identities [1].
Eq. 4a:C = C The first relationship can be proven simply by considering reflection strongly detuned from resonance, so that s − = Cs + . Under a time-reversal operation, we then find s * + =Cs * − . Taking the conjugate, left-multiplying by (C * ) −1 , and using that C is unitary if the direct pathway is lossless, we findC = C .
Eq. 4d:γr = γr To prove Eq. 4d, we again invoke the time-reversed scenario of a decaying cavity. Without input, s + = 0, there are no reflections in the time-reversed case:s − =s * + = 0. Also, as mentioned earlier, the incident signal in the timereversed case iss + =s * − . Hence, we can write for Eq. 2 in the main text:Cs * − +da * = 0. Given that s * − = d * a * , we findCd * = −d, which, when invoking unitarity of C, yields d † d =d †d , and thus Eq. 4d.
Eqs. 4b,c:d = k andk = d To prove Eqs. 4b,c, consider a lossless cavity with an initial amplitude decaying into the ports, while s + = 0. In this scenario, both a and s − decay exponentially with complex frequency ω 0 − iγ. If we then reverse time, we excite the cavity with an exponentially growing wave with amplitudes + = s * − and frequency ω 0 + iγ. Now, starting from the equation of motion (iω − iω 0 + γ)a = k s + , we find for the time-reversed scenario at frequency ω = ω 0 + iγ: Using γ = γ r and γ r =γ r (which we have shown in the previous subsection), yields 2γ r = k d * . Taking the complex conjugate of 2γ r = k d * and combining it with 2γ r = d † d, we find: As long as d is a non-zero vector, this impliesk =d. Starting from the time-reversed scenario and following the same analysis, we can also confirm thatd =k.
Finally, by using the fact that C is unitary, we obtain Eq. 4f from Eq. 4e: d † d = k † k. It is interesting to point out that there are various ways to derive this fluctuationdissipation relation: one may also prove it using balance of power, or more rigorously, using stochastic methods [2]. Furthermore, it is important to point out that in the case of internal absorption, γ i > 0, the full fluctuationdissipation relation needs to be amended to include absorptive dissipation as well.
In this context, we should stress that while these proofs rely on the unitarity of C, the relations in Eq. 4 can also be more generally applied to lossy systems by considering loss as (an) additional port(s). This is specifically demonstrated by a heuristic derivation of Eq. 4e for the system with a dissipative wedge mode in Section 5 of the Supplementary Information.

DETAILS OF CAVITY GEOMETRIES AND MATERIALS
The unidirectional waveguides we use to study nonreciprocal cavities are based on the work by Shen et al. [3,4].
The materials used in all of our simulations are the same: silicon (Si) and indium-antimonide (InSb). We use a constant permittivity for Si: ε Si = 11.68ε 0 , and for InSb we use the transversely (ẑ) magnetized permittivity tensor [3]: where ε ∞ = 15.6ε 0 and Here ω p = 4π × 10 12 rad/s is the plasma frequency, ν = 5 × 10 −3 ω p rad/s is the collision frequency, and ω c = eB/m = −0.25ω p is the cyclotron frequency (corresponding to a static magnetic field bias of 0.25 T in the −ẑ direction). In all simulations the waveguide has the same dimensions: the total height of the waveguide is 30 µm, filled with 18 µm of InSb at the bottom and 12 µm of Si at the top. The bottom and top walls of the waveguide are perfect electric conductors (see the next section). For Fig. 1 in the main text we place a cavity behind the termination of the waveguide, which is 20 × 30 µm and resonant at 1.52 THz. The cavity is connected to the waveguide through a small opening with height of 0.5 µm and a width of 0.1 µm. See Fig. S1a for a schematic drawing  Fig. 1), where the inset shows the geometry of the cavity opening. b) Schematic for the second system discussed in the main text (and in Fig. 2), where almost everything is identical as in the previous cavity except for a displacement of the cavity and a larger size. (c) The cavity is the same as in (b), except that now the cavity is on top and slightly displaced from the opening. In the waveguide, the opening is placed at a position where the forward and backward fields cancel out.
of this geometry. We operate the waveguide in the unidirectional regime, with a pulse centered at 1.5 THz and a bandwidth of 0.16 THz (see next section). For Fig. 2 in the main text, we operate the waveguide in the bidirectional regime, below 1.25 THz [3]. We increase the cavity size to 20 × 34.4 µm so that it is resonant at 1.24 THz, and while we maintain the opening at the same position in the waveguide (3 µm) from the top wall, we shift the cavity upwards so that the opening is closer to the middle of the cavity (which increases the cavity Q-factor). The displacement between the middle plane of the cavity and the middle of the opening is 2 µm. For simplicity of analysis (see Section 4) we set the collision frequency to ν = 0 but have changed the PEC termination to a lossy impedance boundary condition with a conductivity of σ = 10 4 S/m and ε = ε 0 . See Fig. S1b for schematic details on this geometry.
The final geometry we have studied, discussed in Fig. 3 in the main text, has identical features to the previous cavity, except that the cavity is now positioned on top of the waveguide and that the conductivity of the termination has been reduced to 2093 S/m. The cavity opening is centered 60 µm away from the termination, where there is a zero in the magnetic field (see Fig. S1c).

ANISOTROPIC FDTD ALGORITHM
To perform full-wave simulations of the cavity in the time-domain we employ a home-built FDTD algorithm. For the geometries of interest, it is important for the algorithm to support i) Drude-model dispersion, and ii) anisotropic materials. To incorporate anisotropic Drude dispersion, we employ the auxiliary differential equation (ADE) method, which captures the dispersion in an additional equation for the current density [5] but implemented so that it supports anisotropic materials. Starting from Ampere's law in differential form: where From the permittivity of InSb in Eqs. 2-S5, we can write for the susceptibility so we find for the current density: Inverting the matrix in this equation to bring it to the other side: Cancelling out the denominator, and using iωJ → ∂J ∂t , we finally find: This auxiliary differential equation incorporates the dispersion, and when converted into an update equation it can be added to a regular FDTD algorithm. We use the material parameters as discussed in the previous section, and a mesh of 10 by 10 nm.

HEURISTIC DERIVATION OF CMT FOR CAVITY AND NONRECIPROCAL WAVEGUIDE
In the unidirectional waveguide discussed in Fig. 1 in the main text there is no reflected power over the whole frequency range of interest, and the direct reflection coefficient C is therefore 0. We thus obtain an additional coefficient for the direct excitation of the wedge mode (which is responsible for absorbing the incident power), |C w | = 1 − |C| 2 , so that without the cavity the power absorbed by the wedge mode is P w = |C w s + | 2 . While there is no backwards mode, even in the case of a lossless cavity as in Fig. 1 in the main text, the resonance can (and must be able to) dissipate power by exciting the wedge mode, which for the cavity results in an additional loss rate γ w : It turns out that it is crucial to consider that both the incident wave and the resonance can excite the wedge mode, and that in due process they can thus interfere. The total power absorbed by the wedge mode in general is therefore given by where d w relates the excitation of the wedge mode to the cavity amplitude, with |d w | 2 = 2γ w (which can be shown from balance of power in the case that γ r = γ i = 0). If we then reconsider balance of power expanding |s + | 2 = |Cs + + d r a| 2 + 2γ i |a| 2 + |C w s + + d w a| 2 , (S16) which becomes Inserting the equation for the cavity amplitude and setting ω = ω 0 , we find: which becomes Considering that k * r k r is real, we can write: This equation is equivalent to Eq. 4d in the main text for this two-port system and it demonstrates that k and d can indeed be different, as also observed in the simulation, as long as there is an additional port. In the case that d = 0 (when the waveguide is unidirectional), we again find that |k r | = |d w |, which means that power can enter the cavity at the same rate that the cavity can dissipate it via the wedge mode. We stress again that it is crucial for the input and cavity to be able to interfere at the wedge. If this were not the case, and for example the wedge mode would dissipate incoherently as P w = |C w s + | 2 + |d w a| 2 , this would result again in the requirement |d r | = |k r |. In fact, in this scenario it is not necessary to consider the term γ w separately from γ i , and a general description can be obtained without explicit consideration of P w and γ r . This fact makes the wedge mode alike an additional channel, rather than simply an additional dissipative process.

USING COMSOL FOR NONRECIPROCAL WAVEGUIDES
For the results in Figs. 2-3 (as well as the intensity plots in Fig. 1) in the main text we use COMSOL rather than our FDTD algorithm, because we are interested in complex amplitudes of the ingoing and outgoing waves with respect to the cavity amplitude. In the following we describe how to use COMSOL for simulations with nonreciprocal media and the fitting procedure to obtain the coefficients discussed in the main text. For reciprocal media it would be most convenient to use COMSOL port boundary condition to launch and accept incoming and outgoing modes and to determine the reflection coefficient from the structure. However, COMSOL's numerical ports do not work with nonreciprocal waveguides, because the incoming and outgoing modes are different. We therefore applied the following procedure to obtain the complex mode and cavity amplitudes, assuming that we are modeling one of the structures in Fig. S1: 1. Create the geometry and add a port on the left end of the waveguide, with "wave excitation" set to "on".
2. To the left of where this port is, add a small second domain that is essentially an extension of the waveguide. Create a second Electromagnetic Waves, Frequency Domain that applies only to this domain and make sure that in the main physics domain this waveguide extension is excluded.
3. In the second physics domain, set this waveguide extension up so that there is a port on the right boundary (overlapping with the port in the main physics domain). All other boundaries can be PEC. 4. Back in the first domain, add a second port on the same boundary, set it to a user defined port and modify the expressions for the electric mode field and propagation constant so that they are obtained from the other physics domain: e.g. "emw2.tEmodex 3", etc. Make sure wave excitation is set to "off" in this port. In this main physics domain there are now separate ports for ingoing and outgoing modes.

5.
Although not strictly necessary (and slightly more complicated), it is possible to modify the weak expressions in the launching port so that it does not try to accept the returning mode. It requires adding an additional port for the electric field (following a similar procedure as before, but now in the same physics domain), and then add 0 * before all if's in emw.PortConstrx, emw.PortConstry, emw.PortConstrz, emw.PortConstrx weak, emw.PortConstry weak, emw.PortConstrz weak, and the weak expression for the domain computation. 6. Use an overlap integral on the input boundary to determine the input and output phase and amplitudes. Because these modes are nonreciprocal, the overlap integral is different [6]: Here, e T m and h T m are the transverse mode profiles, while E T and H T are the transverse field profiles of the full solution.
7. Having obtained the input and output complex amplitudes, we can obtain the complex cavity amplitude from an integral over the stored electromagnetic energy in the cavity (due to the small opening and PEC walls virtually no energy is stored outside of the cavity) for the amplitude, and a field monitor in the center of the cavity for the phase.

FITTING SIMULATIONS WITH THE COUPLED-MODE THEORY MODEL
We now describe the fitting procedure to obtain the complex coefficients reported in the main text. First, starting with the system shown in Fig. S1b but without a magnetic field bias (so that it is reciprocal), we perform a frequency sweep for three different cavity sizes: one for the resonant cavity size (20 × 36.04 µm), and for cavity sizes one micron longer and shorter. We use these nonresonant cavities to determine the direct reflection path, by averaging C = S − /S + for both simulations. Then, we obtain ω 0 and γ from a lorentzian fit of the stored energy in the cavity, which is shown in Fig. S2a: ω 0 = 1.24 THz and γ=96.7 MHz. To obtain k r we then proceed to the complex cavity amplitude, which we fit using a = k r s + i(ω − ω 0 ) + γ (S22) We use ω 0 and γ from the previous fit and s + is input from COMSOL, and hence only k r is a fitting parameter. The fit of the complex cavity amplitude is shown in Fig. S2b, for k r = (2.22 − 1.78) × 10 4 rad/s. To obtain d r we fit the reflected amplitude: Here, again, d r is the only fit parameter, and the rest we have obtained directly from COMSOL. Fig. S2c shows the resulting fit, for d r = (222 − 1.79) × 10 4 rad/s. As expected, k r = d r but for a very small difference (the ratio is d r /k r = 1.0010 − 0.0007i). The small differences most likely originates from estimating C (by changing the cavity size) and a (by assuming that all of the stored energy is inside the cavity). The fact that we find k r = d r in the reciprocal regime thus validates our method. If we now turn the magnetic field bias on again, we find the In all cases the fit is excellent, while fitting only one parameter (b,c,e,f). In (a,d) we fit three parameters: the loss rate, center frequency, and maximum stored energy (which we don't use in our analysis).
fits shown in Fig. S2d-e, with ω 0 = 1.24 THz and γ = 81.2 MHz (indicating that the decay rate has changed, which is to be expected given the change in the waveguide mode field profiles). By fitting the complex cavity and reflection amplitudes we find k r = (2.65 + 0.308) × 10 4 rad/s and d r = (0.667 + 2.14) × 10 4 rad/s. Now, clearly, the values are different: their ratio is dr kr = 0.34 + 0.77i. * aalu@gc.cuny.edu