Synchronization-based lossless non-reciprocal scattering]Synchronization-based lossless non-reciprocal scattering

1]Tiemo Pedergnana

1]Abel Faure-Beaulieu

[2]Romain<EMAIL_ADDRESS>
[1]Nicolas<EMAIL_ADDRESS>
[1]Department of Mechanical and Process Engineering, ETH Zürich, Zürich, Switzerland

[2]Institute of Electrical and Micro Engineering, EPFL, Vaud, Switzerland

Breaking the reciprocity of wave propagation is a problem of fundamental interest, and a mucht-sought functionality in  practical applications, both in photonics and phononics. Although it has been achieved using resonant linear scattering from cavities with broken time-reversal symmetry, such realizations have remained inescapably plagued by inherent passivity constraints, which make absorption losses unavoidable, leading to stringent limitations in transmitted power. In this work, we solve this problem by converting the cavity resonance into a limit cycle, exploiting the uncharted interplay between non-linearity, gain, and non-reciprocity.  Remarkably, strong enough incident waves can synchronize with these self-sustained oscillations and use their energy for amplification. We theoretically and experimentally demonstrate that this mechanism can simultaneously enhance non-reciprocity and compensate absorption. Real-world acoustic scattering experiments allow us to observe perfect non-reciprocal transmission of audible sound in a synchronisation-based 3-port circulator with full immunity against losses.

§ INTRODUCTION

Wave propagation is said to be reciprocal if it is symmetric with respect to exchanging the locations of any sender-receiver pair. Due to its close links to time-reversal symmetry, reciprocity is ubiquitous in physics [1, 2, 3]. Breaking reciprocity is necessary in applications which aim to guide waves in one direction while isolating another, for example to protect a source from unwanted reflections, or to induce strong topological protection against obstacles or imperfections [4, 5]. Such non-reciprocal transmission is of high practical interest and has recently been reviewed for electromagnetic [6, 7] and mechanical [8, 9] waves.

Analogous to how light can be non-reciprocally routed using biased magneto-optical media [10], unidirectional acoustic waves can be induced in a resonant cavity with a biased air flow [11]. Both aforementioned effects are based on scattering waves at a pair of isolated, degenerate cavity modes whose symmetry is broken by a frequency splitting similar to the electronic Zeeman effect [12], where an external magnetic field lifts the degeneracy of the atom's eigenstates. The exploitation of this effect to construct a circulator is well-known, and illustrated in Figs. <ref>a and <ref>b, which represents the state of the art in linear, resonance-based non-reciprocal scattering. Of course, any macroscopic realizations, albeit inspired from atomic-scale quantum phenomena, will suffer from dissipative effects [13, 14]. The present study proposes a solution to remedy this fundamental problem by unveiling a synchronization-based mechanism for achieving lossless non-reciprocal transmission. Our method requires no internal active elements to control the field [15, 16].

Instead, we show here that a lossless circulator can be designed based on the non-linear interference between a wave and a cavity limit cycle, which emits self-sustained radiation in the absence of an external source (Fig. <ref>c). The concept proposed in this work is sketched in Fig. <ref>d. Not only can a wave impinging on a cavity limit cycle extract energy from its radiated signal by synchronization, but the non-linear interference simultaneously enhances non-reciprocal isolation between the ports. In a three-port circulator with a properly tuned bias, the above scenario can enable lossless non-reciprocal transmission, overcoming the dissipation induced by low $Q$ factors in resonant cavities. After extending classical coupled-mode theory [17, 18] to the non-linear case, we can explain and predict the synchronization-based transmission using only a few parameters. To validate our theory of synchronization-based lossless non-reciprocity, we finally present an experimental realization of a loss-immune acoustic circulator, which exploits flow-induced whistling [19, 20, 21] to create spinning acoustic limit cycles [22, 23].

no external sources

(symmetry intact)

non-reciprocal scattering

(symmetry broken)

limit cycle

Iphv 0
Jphv time
radiation losses

from all ports

lossy non-reci-

procal scattering

self-sustained radia-

tion from all ports

lossless non-reci-

procal scattering

1phv $1$
0phv $0$
2phv $T_{1\shortrightarrow 2}$ and
Lphv $T_{1\shortrightarrow 3}$
Mphv $T_{1\shortrightarrow 3}$
Nphv $T_{1\shortrightarrow 2}$
wphv $f$
nphv 1
ophv 2
pphv 3
qphv 4
3phv min.
4phv 0
5phv max.
6phv Amplitude
Fphv Port 1
Gphv Port 2
Hphv Port 3
Cphv Port 1
Dphv Port 2
Ephv Port 3
9phv output in sync
Zphv or
>phv decaying disturbance

Concept of synchronization-based lossless non-reciprocal scattering in a circulator. a, b State of the art. Inducing a Zeeman-type bias in a resonator cavity (a) leads to non-reciprocal transmission of incident waves (b). Due to the passivity of the resonance, non-reciprocal transmission will be lossy, even if the isolation of the blocked direction is perfect. c A limit cycle in a round cavity with multiple ports continuously radiates a signal through the ports. d A harmonic wave incident on the cavity can synchronize with and gain energy from the limit cycle's emissions. In the presence of a bias, this enables lossless non-reciprocal transmission.  e In absence of external forcing, the acoustic field is either spinning right or left.  f Under external forcing, the spinning acoustic field becomes a standing wave.

§ NON-LINEAR WAVE-MODE COUPLING

In this section, we describe the non-linear interference between the cavity modes and the incident wave in a circulator with three ports, highlighting the key ingredients required to achieve lossless non-reciprocal transmission. The concept and the theoretical foundations are not limited to acoustic circulators. As in [10], cross-coupling between the modes is neglected under the assumption that the dominant mechanisms are the coupling between the incident and the outgoing wave (W$\shortrightarrow$W), coupling from the incident wave to the cavity modes (W$\shortrightarrow$M) and coupling from the cavity modes to the outgoing wave (M$\shortrightarrow$W). As discussed in the next section, only one mode ($a_-$) is undergoing a limit cycle. The non-linear coupling between the cavity modes $a=(a_+,a_-)^T\in\mathbb{C}^2$ and the incident wave $\vert s_\mathrm{in}\rangle\in\mathbb{C}^3$ is governed by the following system:

\begin{eqnarray}
\dot{a}_\pm&=&\underbrace{i\omega_{\pm} a_\pm }_{\text{oscillation}} -\underbrace{\gamma a_\pm}_{\text{dissipation}} +\underbrace{\frac{\beta_\pm a_\pm}{1+\kappa \vert a_\pm\vert^2}}_{\text{amplification}}  +\underbrace{D^\dag_\pm \vert s_{\mathrm{in}}\rangle}_{\text{W}\shortrightarrow \text{M}},\label{modal dynamics 2m}\\
\vert s_{\mathrm{out}}\rangle&=&\underbrace{C \vert s_{\mathrm{in}}\rangle}_{\text{W}\shortrightarrow \text{W}} + \underbrace{D a}_{\text{M}\shortrightarrow \text{W}}, \label{input-output 2m}
\end{eqnarray}

where $\omega_{\pm}$ are the eigenfrequencies, $\gamma$ is the decay rate, $\beta_\pm\in\mathbb{R}^+$ are the linear gain coefficients, $\kappa \in \mathbb{R}^+$ is the non-linearity constant, $C\in \mathbb{R}^{3\times 3}$ is the background scattering matrix, $\vert s_\mathrm{out}\rangle\in\mathbb{C}^3$ is the outgoing wave, $D\in\mathbb{C}^{3\times 2}$ is the coupling matrix. A dot denotes the time derivative and $\dag$ is the Hermitian transpose. Explicit expressions of $D$ and $C$ are given in the Supplementary Notes [24]. Note that Eq. (<ref>) corresponds to two decoupled equations, one for $a_+$ and one for $a_-$. The eigenfrequencies $\omega_\pm$ are generally assumed to be quasi-degenerate ($\omega_+\approx \omega_-$) and well-isolated in the cavity spectrum, meaning that there are no nearby resonances. Analogously, the linear gain coefficients are assumed to be similar: $\beta_+\approx\beta_-$. Due to symmetry, the decay rate $\gamma$ and the non-linearity coefficient $\kappa$ are assumed to be the same for both modes. To compute the steady-state forced response of $a_\pm$, we assume that the modes are perfectly synchronized with the incident wave. The functional form of the amplification term in Eq. (<ref>) is motivated by experimental results [21] and this type of saturable gain was also used in  [25]. For a wave incident from the first port, $\vert s_\mathrm{in}\rangle=s e^{i\omega t} \vert 1\rangle$, where $s$ is the incident wave amplitude, $\omega=2\pi f$ is the angular forcing frequency and $\vert 1\rangle=(1,0,0)^T$ is a unit vector. We define the nondimensionalized amplitude as $\tilde{s}=s \sqrt{\kappa/\beta}$, where $\beta=(\beta_++\beta_-)/2$ is the arithmetic mean of the linear gain coefficients. Using the definition of the scattering matrix $S$ given by $\vert s_\mathrm{out}\rangle=S\vert s_\mathrm{in}\rangle$, a general expression for $S$ can be found using the Moore-Penrose pseudoinverse [25]:

\begin{eqnarray}
S=\dfrac{\vert s_\mathrm{out}\rangle\langle s_\mathrm{in}\vert}{\langle s_\mathrm{in}\vert s_\mathrm{in}\rangle}. \label{general scattering matrixm}
\end{eqnarray}

In general, $S$ is a non-linear function of $\vert s_\mathrm{in}\rangle$ because the outgoing wave $\vert s_\mathrm{out}\rangle$ is a linear function of the cavity modes $a_\pm$, which, in turn, depend non-linearly on $\vert s_\mathrm{in}\rangle$.

In the following, we provide analytical results which numerically illustrate the mechanism visualized in Fig. <ref>. We recall that, in order to synchronize with the limit cycle so as to profit from the energy of its self oscillations, the incident amplitude must be sufficiently large. In the limit of infinitely large amplitudes $\vert a_\pm\vert$, the modal equations (<ref>) reduce to the classical temporal coupled mode theory which describes linear scattering at a resonant cavity [17, 18, 10]. This limit plays an important role in the non-linear transmission scenario we propose. One can imagine the limit of infinitely large amplitude as a dissipative ground state which we seek to overcome by synchronizing with the cavity limit cycle. In practice, the goal is to increase the incident wave amplitude to a large enough value such that the resulting scattering matrix is a weakly perturbed version of its infinite amplitude limit, and then to tune the amplitude and the modal parameters until the desired scattering state is reached. By expanding the modal dynamics (<ref>) in a perturbation series in $s^{-1}$, one can compute the leading-order correction to the scattering matrix in the limit of large incident wave amplitudes. The derivation of the non-linear scattering matrix and its large amplitude approximation is provided in the Supplementary Notes [24]. To model dissipation in our model, we introduce the unitarity factor $\sigma\in\mathbb{R}^+$, which is equal to $1$ for a perfectly reversible process and takes values between 0 and 1 if the scattering is lossy. Thus, $\sigma$ is a measure of the dissipation present in a real-world system that must be overcome to achieve lossless scattering. Next, we assume that the Zeeman splitting and incident wave frequency satisfy $\omega_\pm=\omega\pm\sqrt{3}\gamma\sigma/(\sigma+2)$, which corresponds to perfect matching (see the Supplementary Notes [24]) Computing the synchronized forced response of the modal amplitudes $a_\pm = \vert a_\pm \vert e^{i\omega t+\varphi_\pm}$ for an incident wave from port  “1”, $\vert s_\mathrm{in}\rangle=s e^{i \omega t}\vert 1 \rangle$, yields the amplitude transmission coefficients for large normalized amplitudes $\tilde{s}$ [24]:

\begin{eqnarray}
T_{1\shortrightarrow 2}&\approx& \bigg \vert \dfrac{1-\sigma}{3} -\dfrac{(2\sigma^2+2\sigma-1)}{  \tilde{s}^2 (\sigma^2 + \sigma + 1)}\bigg \vert+O(\tilde{s}^{-4}), \label{T12lim4m} \\
T_{1\shortrightarrow 3}&\approx&\bigg \vert \dfrac{ (1+2\sigma)  }{3}+\dfrac{(1+4\sigma+\sigma^2)}{ \tilde{s}^2  (\sigma^2 + \sigma + 1)}\bigg \vert+O(\tilde{s}^{-4}).  \label{T13lim4m}
\end{eqnarray}

For mildly dissipative systems with $0.4<\sigma < 1$, and large enough $\tilde{s}$, Eq. (<ref>) shows that the non-linear interference of the incident wave with the cavity limit cycle decreases transmission to the blocked port “2”, eventually reaching an identically zero value. Under these conditions, Eq. (<ref>) indicates that transmission to port “3” is enhanced and reaches values close to or even exceeding unity, depending on level of dissipation $\sigma$. The above results show that qualitatively, this synchronization-based mechanism enables an increase in the non-reciprocal behavior of the system, while at the same time compensating transmission losses. Quantitative comparisons between theory and experiments will be provided in Sec. <ref>.

§ EXPERIMENTAL VALIDATION

For our scattering experiments, we built and characterized a circulator which consists of a deep annular cavity subject to a swirling flow [23] and is equipped with three equispaced ports along its circumference. The bulk flow speed, the swirl number of the flow and the amplitude of the incident wave through one of the ports can be set individually to tune the non-linear scattering matrix. Three waveguides are connected to the cavity ports, thus forming an acoustic circulator. The experimental setup is shown in Fig. <ref>a. Technical information is given in Appendix <ref>. Let us first consider Fig. <ref>a in the absence of acoustic forcing. A swirling flow is imposed in the blue pipe, and passes through the center of the axisymmetric cavity. At suitable conditions, the swirling jet in the cavity exhibits self-sustained oscillations thanks to the intrinsic coupling between one of its hydrodynamic modes and a spinning acoustic mode of the cavity. The trapped aeroacoustic mode producing this cavity whistling is shown in Fig. <ref>b, where one can see the acoustic pressure field that interacts with the coherent vorticity fluctuations. The latter fluctuations develop along the shear layer in the cavity aperture region, and were measured with particle image velocimetry (see Fig. <ref>c). A bifurcation diagram of the acoustic pressure inside the cavity, depicting the transition between the quiet and the whistling phase, is presented in Appendix <ref>.

Once an aeroacoustic limit cycle has been established and the cavity is whistling, acoustic forcing with a specified amplitude is imposed from a compression driver mounted on one of the ports. The incident and outgoing signals are reconstructed using microphones mounted on the ports to obtain the transmission coefficients. For sufficiently strong and in-tune forcing, the cavity limit cycle can synchronize with the forcing, as demonstrated in Appendix <ref>. Therefore, this experimental set-up is a realistic representation of the model assumptions made in Sec. <ref>.

aphv Port 1
bphv Port 2
cphv Port 3
ephvSwirling air flow [$\shortrightarrow$]

-phv$\tilde{\Omega}_\theta$ [s-1]
nphv Source
Bphv 70 cm
Zphv cylindrical cavity
?phvflow direction

Acoustic circulator based on a spinning cavity limit cycle. a Experimental setup. A swirling flow with fixed axial and azimuthal components is imposed in the blue pipe which passes through a cylindrical cavity, leading to a limit cycle (whistling) at around $800$ Hz. In contrast with classic wind instruments whose sound radiation originates from self-oscillating standing modes, our special whistle involves spinning modes. For the scattering experiments, three 70 cm long square waveguides are connected to the upstream wall of the round cavity at equispaced angular positions ($120^\circ$ apart). These waveguides are equipped with sound absorbing foam at their left termination in order to provide quasi-anechoic conditions. For the scattering experiments, incident waves are sent to the cavity using a compression driver (source) mounted on one of the waveguides. Microphones in the waveguides are used to reconstruct the incident and outgoing waves. b Acoustic mode of the circulator cavity, also showing isosurfaces of the azimuthal component of the coherent vorticity fluctuations $\tilde{\Omega}_\theta$ (see Appendix <ref> and <ref>). c Zoomed-in view of this instantaneous snapshot of the vorticity fluctuations, which was obtained from stereoscopic particle image velocimetry.

§ RESULTS

In order to facilitate the comparison of theoretical and experimental results, we define the incident wave amplitude relative to the maximum range of amplitudes $[s_\mathrm{min},s_\mathrm{max}]$, imposed in the scattering process: $\Delta s=(s-s_\mathrm{min})/(s_\mathrm{max}-s_\mathrm{min})$, and analogously for $\Delta \tilde {s}$. We show in Fig. <ref> that our theoretical model reproduces the non-linear scattering experiments with remarkable accuracy using only few parameters. The ranges of excitation amplitudes used in the model and in the experiments were $\tilde{s}\in[0.04,2.1]$ and  $s\in [10,120]$ Pa, respectively. No fit was performed to adjust the model parameters, but they were manually defined to match the experiments. The parameter values used in the model were $\sigma=0.61$, $\gamma=27.5$ s-1, $f_+=\omega_+/2\pi=804.3$ Hz, $f_-=\omega_-/2\pi=800.8$ Hz, $\beta_+=21.2$ s-1, $\beta_-=30.3$ s-1, $\kappa=1$ (in units of $[a^{-2}]$). For these values, perfect matching in the model is achieved at 802.5 Hz, similar to what is observed in the experiments (see also Appendix <ref>).

nphvfrequency (Hz)
mphv$ s_\mathrm{in} $ (Pa)
gphv10 120
aphv$s_\mathrm{in}$ (Pa)
yphvincident wave amp. $ {s}_\mathrm{in} $ (Pa)
vphvtransmission 1$\shortrightarrow$3 (dB)
tphv-20  -11  -1.4   8.1  18
Ephv-41  -29  -16   -3.2  10
Gphvisolation (dB)
Zphv$s_\mathrm{in}$=100 Pa
Yphv$s_\mathrm{in}$=50 Pa
Xphv$s_\mathrm{in}$=10 Pa
Aphv780  790 800
810 820
Cphvfrequency (Hz)
1phvExperiments  Theory
2phvincreasing excitation amplitude
Xphv$\Delta s =0$
Yphv$\Delta s =0.36$
Zphv$\Delta s =0.82$
3phv$\Delta \tilde{s} =0$
4phv$\Delta \tilde{s} =0.36$
5phv$\Delta \tilde{s} =0.82$
7phv$T_{1\shortrightarrow 2}$
8phv$T_{1\shortrightarrow 3}$

Experimental validation of a perfect lossless circulator. Amplitude transmission coefficients for increasing incident wave amplitudes $\Delta \tilde{s}$ and $\Delta s$ (expressed relative the intervals $\tilde{s}\in[0.04,2.1]$ and $s\in [10,120]$ Pa, respectively; see Sec. <ref>). Top row: Theoretical results obtained from the non-linear coupled-mode theory described in Sec. <ref>. Bottom row: experimental measurements performed on the setup described in Sec. <ref>. In both cases, a harmonic incident wave is imposed in port “1”. The blue and black curves refer to the transmitted wave amplitudes at ports “2” and “3”, respectively, relative to the incident wave amplitude and frequency.

For small incident wave amplitudes, sharp peaks are visible because the self-sustained radiation of the limit cycle dominates transmission through the ports. As the incident wave amplitude is increased, a non-reciprocal, lossless scattering state emerges. To verify the robustness of the experimental results shown in Fig. <ref>, we investigated the sensitivity of the lossless, non-reciprocal state (see Appendix <ref>). For this, we repeated the scattering experiments for a fixed value of $\Delta s=9$ while varying the axial and azimuthal components of the swirling flow. The results show that there exists a continuum of low-loss, non-reciprocal scattering states in the vicinity of the operating point treated in this study. Therefore, our measurements evidence the relevance of loss-compensated non-reciprocity as a robust synchronisation phenomenon which persists under small variations of the operating conditions. As limit cycles are ubiquitous in dynamical systems with saturable gain, we surmise that synchronization-based lossless non-reciprocity may be transposed to other physical platforms, from electronics to non-Hermitian photonics. Semiconductor-based electronic amplifiers and continuous optical gain media may allow for a direct implementation of our concept in high-frequency circuits or silicon nanophotonics. More fundamentally, our loss-compensation mechanism may open new experimental opportunities in topological physics, where time-reversal symmetry plays a pivotal role in the robustness of chiral edge states, which can be induced in circulator networks [4, 5]. The complex interplay between non-linearity and gain to produce limit cycles capable of breaking reciprocity remains a promising avenue, and we envision that extensions to periodic lattices, or systems coupled to a radiation continuum may lead to a wealth of unforeseen progress in a field that used to be thought as fundamentally hindered by absorption losses.

§.§ Data availability

The data that support the findings of this study are available at <https://polybox.ethz.ch/index.php/s/CLxLxoDsWkSuIgq>.

§.§ Code availability

The codes that support the findings of this study are available at <https://polybox.ethz.ch/index.php/s/CLxLxoDsWkSuIgq>.

§.§ Contributions

N.N. conceptualized and supervised the research. T.P. and A.F. conducted the acoustic and flow measurements. T.P. and N.N. developed the theory. T.P. performed data analysis and validated the model. T.P., A.F., R.F., and N.N. discussed the results. T.P. wrote the first draft of the article. N.N. and R.F. reviewed and edited it. All authors approved the final version of the article.

[1]
authorOnsager, L.
titleReciprocal relations in irreversible processes. I.
journalPhys. Rev.
volume37 (4), pages405 – 426
10.1103/PhysRev.37.405 .

[2]
authorOnsager, L.
titleReciprocal relations in irreversible processes.
journalPhysical Review
volume38 (12), pages2265 – 2279
10.1103/PhysRev.38.2265 .

[3]
authorCasimir, H.
titleOn Onsager's principle of microscopic
journalRev. Mod. Phys.
volume17 (2-3), pages343 – 350
10.1103/RevModPhys.17.343 .

[4]
authorZhang, Z., authorDelplace, P. &
authorFleury, R.
titleSuperior robustness of anomalous non-reciprocal
topological edge states.
volume598 (7880), pages293 – 297
10.1038/s41586-021-03868-7 .

[5]
authorZhang, Z., authorDelplace, P. &
authorFleury, R.
titleAnomalous topological waves in strongly amorphous
scattering networks.
journalSci. Adv. volume9 (12)
10.1126/sciadv.adg3186 .

[6]
authorVerhagen, E. & authorAlù, A.
titleOptomechanical nonreciprocity.
journalNat. Phys.
volume13 (10), pages922 – 924
10.1038/nphys4283 .

[7]
authorCaloz, C. et al.
titleElectromagnetic nonreciprocity.
journalPhys. Rev. Appl.
volume10 (4) (year2018).
10.1103/PhysRevApplied.10.047001 .

[8]
authorNassar, H. et al.
titleNonreciprocity in acoustic and elastic materials.
journalNat. Rev. Mater.
volume5 (9), pages667 – 685
10.1038/s41578-020-0206-0 .

[9]
authorRasmussen, C., authorQuan, L. &
authorAlù, A.
titleAcoustic nonreciprocity.
journalJ. Appl. Phys.
volume129 (21) (year2021).
10.1063/5.0050775 .

[10]
authorWang, Z. & authorFan, S.
titleMagneto-optical defects in two-dimensional photonic
journalAppl. Phys. B
volume81, pages369–375
10.1007/s00340-005-1846-x .

[11]
authorFleury, R., authorSounas, D.,
authorSieck, C., authorHaberman, M. &
authorAlù, A.
titleSound isolation and giant linear nonreciprocity in a
compact acoustic circulator.
volume343 (6170), pages516–519
10.1126/science.1246957 .

[12]
authorCohen-Tannoudji, C., authorDiu, B. &
authorLaloe, F.
titleQuantum mechanics. Vol. 1
(publisherHermann, Wiley, NY, 1991, year2009).

[13]
authorCallen, H. B. & authorWelton, T. A.
titleIrreversibility and generalized noise.
journalPhys. Rev.
volume83 (1), pages34 – 40
10.1103/PhysRev.83.34 .

[14]
authorKubo, R.
titleThe fluctuation-dissipation theorem.
journalRep. Prog. Phys.
volume29 (1), pages255 – 284
10.1088/0034-4885/29/1/306 .

[15]
authorDinc, T. et al.
titleSynchronized conductivity modulation to realize
broadband lossless magnetic-free non-reciprocity.
journalNat. Commun.
volume8 (1) (year2017).
10.1038/s41467-017-00798-9 .

[16]
authorZhang, H., authorLi, R., authorBao,
Y., authorLiu, X. & authorZhang, Y.
titleTotal acoustic transmission in a honeycomb network
empowered by compact acoustic isolator.
journalSci. Rep. volume13 (1)
10.1038/s41598-023-28097-y .

[17]
authorFan, S., authorSuh, W. &
authorJoannopoulos, J.
titleTemporal coupled-mode theory for the fano resonance
in optical resonators.
journalJ. Opt. Soc. Am. A
volume20 (3), pages569–572
10.1364/JOSAA.20.000569 .

[18]
authorSuh, W., authorWang, Z. & authorFan,
titleTemporal coupled-mode theory and the presence of
non-orthogonal modes in lossless multimode cavities.
journalIEEE J. Quantum Electron.
volume40 (10), pages1511–1518
10.1109/JQE.2004.834773 .

[19]
authorHowe, M.
titleThe dissipation of sound at an edge.
journalJ. Sound Vib.
volume70 (3), pages407–411
10.1016/0022-460X(80)90308-9 .

[20]
authorElder, S., authorFarabee, T. &
authorDeMetz, F.
titleMechanisms of flow-excited cavity tones at low mach
journalJ. Acoust. Soc. Am.
volume72 (2), pages532–549
10.1121/1.388034 .

[21]
authorBourquard, C., authorFaure-Beaulieu, A. &
authorNoiray, N.
titleWhistling of deep cavities subject to turbulent
grazing flow: Intermittently unstable aeroacoustic feedback.
journalJ. Fluid Mech.
volume909 (year2020).
10.1017/jfm.2020.984 .

[22]
authorFaure-Beaulieu, A., authorXiong, Y.,
authorPedergnana, T. & authorNoiray, N.
titleSelf-sustained azimuthal aeroacoustic modes. Part
1. Symmetry breaking of the mean flow by spinning waves.
journalJ. Fluid Mech.
volume971, pagesA21 (year2023).
10.1017/jfm.2023.352 .

[23]
authorFaure-Beaulieu, A., authorPedergnana, T. &
authorNoiray, N.
titleSelf-sustained azimuthal aeroacoustic modes. Part
2. Effect of a swirling mean flow on the modal dynamics.
journalJ. Fluid Mech.
volume971, pagesA22 (year2023).
10.1017/jfm.2023.348 .

[24]
titleSupplementary notes to “Synchronization-based lossless
non-reciprocal scattering”  (year2023) .

[25]
authorPedergnana, T. & authorNoiray, N.
titleSuperradiant scattering by a limit cycle.
journalPhys. Rev. Appl.
volume20, pages034068
10.1103/PhysRevApplied.20.034068 .

[26]
authorJang, S.-H. & authorIh, J.-G.
titleOn the multiple microphone method for measuring
in-duct acoustic properties in the presence of mean flow.
journalJ. Acoust. Soc. Am.
volume103 (3), pages1520–1526
https://doi.org/10.1121/1.421289 .

[27]
authorBalanov, A., authorJanson, N.,
authorPostnov, D. & authorSosnovtseva, O.
titleSynchronization: From Simple to Complex
(publisherSpringer, Berlin, Heidelberg, year2009).

[28]
authorCoombes, S. & authorBressloff, P.
titleMode locking and arnold tongues in integrate-and-fire
neural oscillators.
journalPhys. Rev. E
volume60 (2), pages2086 – 2096
10.1103/PhysRevE.60.2086 .

§ APPENDIX

§.§ Scattering experiments

To test our theoretical predictions, we designed an aeroacoustic circulator, shown in Fig. <ref>a, which uses a low-Mach air flow to generate acoustic limit cycles. For this, a swirling flow with bulk speed $\overline{U}=U+V$, whose axial ($U$) and azimuthal ($V$) components can be separately tuned (see Figure <ref>), is imposed in the pipe to which the cavity is connected. The axial velocity $U$ is deduced from the mass ﬂow and the temperature in the channel, which are measured with a mass flow meter and a thermocouple. The azimuthal flow is induced using a swirler based on tangential injection of air in the upststream part of the pipe, and also monitored using a mass flow meter. Note that the velocity $V$ corresponds to the contribution of the mass flow through the swirler $\dot{m}_\mathrm{s}$ to the bulk flow in the pipe, and not to the local velocity in the swirler: $V=4\dot{m}_\mathrm{s}/\rho_0 \pi D^2$, where $\rho=1.15$ kg/m3 is the ambient density and $D=4$ cm is the pipe diameter. The ambient speed of sound in the experiments is equal to $c=343$ m/s. Three semi-anechoic waveguides (arms) with square cross-sections are connected to the cavity through compact circular holes. The ends of the arms are filled with acoustic foam and have been measured to reflect less than $13\%$ of the acoustic energy for all considered operating conditions. Three G.R.A.S. 46BD $1/4$” CCP microphones are flush mounted in each of the arms $x=-4.6$ cm, $x=-12.1$ cm and $x=-24.6$ cm, where $x$ is measured from the upstream inner cavity wall. The multi-microphone method was used to reconstruct the traveling waves in the arms at the incident wave frequency from the acoustic measurements [26]. Two microphones of the same type as in the waveguides were flush mounted on the downstream wall of the cavity at the same radial positions $r=9$ cm, separated by 120$^\circ$ degrees, to measure the acoustic limit cycle oscillations in the absence of forcing from any of the waveguides. The incident waves were imposed by a Beyma CP-850Nd compression driver mounted on one of the arms. The voltage of the signal to the drivers can be calibrated at each frequency to achieve a specified acoustic amplitude $s_\mathrm{in}$.

§.§ Acoustic field

The acoustic mode shown in Fig. <ref>b corresponds to the analytical solution of the Helmholtz equation for a circular cavity, $J_1(r Z_1/R)\cos\theta$, where $J_1$ is the Bessel function of the first kind of order $1$, $r$ is the radial coordinate and $Z_1$ is the first zero of the first derivative of $J_1$ with respect to $\zeta$ [23].

§.§ Fluid-dynamical data

Shown in Fig. <ref>c are isosurfaces of the azimuthal component of the vorticity fluctuations $\tilde{\Omega}_\theta$ at the acoustic frequency, which were extracted using stereoscopic particle image velocimetry (PIV) in an earlier study of a limit cycle at $U=25.9$ m/s and $V=0$, and which are included here for illustration purposes. The isosurfaces were produced from postprocessed data originally presented in <cit.>. The type of aeroacoustic instability used here is based on a low Mach flow [19, 20, 21], and as in the case of human whistling, the subtle interplay between the acoustic field and the rotational fluid motion results in finite-amplitude acoustic pressure oscillations in the cavity. To engineer our circulator, we designed a special whistle in the sense that its eigenmodes are spinning inside the cavity [22], in contrast to the classic police whistle or wind instruments like flutes, in which standing acoustic oscillations are found. Furthermore, by imposing a swirl on the cavity flow, we can provide self-sustained aeroacoustic waves that spin against the swirl, be it in the clockwise or in the counterclockwise direction, as explained in [23].

§.§ Bifurcation diagram

The whistling phenomenon observed in the experiments can be described in low-order models, such as the “$-$” branch  in Eq. (<ref>), which, for $\vert s\rangle=0$ and $\beta_->\gamma$, describes a limit cycle arising from a supercritical Hopf bifurcation. Idealized theory predicts that the square of the modal amplitude $a_-$ varies linearly with the bifurcation parameter $\beta_-$ in the vicinity of the bifurcation point. The bifurcation is governed by the system nonlinearity shown in Fig. <ref>. In the experiments, we have two free parameters ($U$ and $V$), which results in a two-parameter Hopf bifurcation, as shown in  Fig. <ref>. To generate this inset, we recorded the acoustic pressure oscillations in the absence of waves incident to the cavity ports with two microphones on the downstream wall of the cavity. These 10 s timetraces were recorded over an equally spaced grid of values of the bulk velocity $\overline{U}$ and the swirl number $\Sigma=V/U$. For each measurement, we then computed the RMS amplitude of each of the two microphones and took the arithmetic mean of these two values, which we denote by $\overline{p}$. We used the average of the RMS amplitudes from two non-diametrically opposed microphones instead of the RMS of a single microphone to avoid the situation that the nodal line of the aeroacoustic mode (compare Fig. <ref>b) lies on the microphone during the measurement. The operating point ($\overline{p}^2=0.61$ hPa2) for the scattering experiments, marked by the red dot, clearly lies in the limit cycle regime. The smaller cyan and green insets show the variation of the whistling intensity $\overline{p}^2$ (gray curves) along the cyan ($\overline{U}=23.1$m/s) and green ($\Sigma=0.49$) lines in Fig. <ref>c, respectively, and the black lines within the smaller insets are linear regressions to extrapolate the bifurcation points. In the cyan inset, the regression itself is trivial because the black line connects exactly two measuring points, but the estimation of the bifurcation point is nevertheless informative.

§.§ Synchronization dynamics

As described in the last section, we now consider a scenario where a limit cycle is reached at the conditions described above. Then, in addition to these self-sustained oscillations, using the compression driver mounted on the top arm, harmonic sound waves are sent through the channel towards the cavity. As we measure the acoustic pressure inside the cavity, at the same azimuthal position as the forcing waveguide, we can observe synchronization phenomena similar to those described in the book of Balanov et al. [27]. To analyze these dynamics, we performed repeated forced experiments at different incident wave amplitudes for fixed amplitude of the limit cycle. The results are presented in Fig. <ref>. Figure <ref>a illustrates the Arnold tongue [28] representing the region in which synchronization of the incident waves takes place. As we can see, for the smallest incident wave amplitude, the synchronized region is almost a double its size at the largest forcing amplitude. This is consistent with the classic paradigm that synchronization requires large amplitudes and small detuning. For a stronger forcing, hence, larger detuning is acceptable to satisfy this double constraint, leading to a wider region in which synchronization can take place. We now take closer look at two specific points, marked by the purple and green dots in Fig. <ref>a at a forcing amplitude of ${s_{\text{in}}}= 100$ Pa, we can visually compare the synchronized and asynchronous modes. For the latter, as we can see from Fig. <ref>b, there are two dominant peaks in the PSD while for the synchronized we observe just one. The non-synchronized mode has a peak at the self-sustained limit cycle around 800 Hz, another at the forcing frequency and a significantly less dominant harmonic peak at around 820 Hz. For the sychronized state, a peak is only observed at the forcing frequency. Figure <ref>c shows the pressure dynamics for a synchronized mode, which can be compared to those of an asynchronous mode as shown in Fig. <ref>d. While the synchronized mode shows roughly constant amplitude, significant beating oscillations accompany the asynchronous state, stemming from the interplay between the coexisting, detuned peaks of the forced and the self-sustained oscillations.

§.§ Zeeman bias modeling

As in [11], the broken symmetry of the theoretical model, expressed here by $\omega_+\neq \omega_-$ and $\beta_+\neq\beta_-$, corresponds in the experiments to the bias introduced by an azimuthal flow. Under the simplifying assumption that $f_\pm=\omega_\pm/2\pi$ and $\beta_\pm$ depend linearly on $\Sigma\neq 0$ and that for $\Sigma=0$, the system is perfectly symmetric, we obtain the variations of these parameters over the swirl number $\Sigma$ (see Figure <ref>). In this case, $a_+$ corresponds to the mode spinning with the swirling flow and $a_-$ to the mode spinning against it. Note that $f_\pm$ and $\beta_\pm$ split in opposite direction in our model, which is consistent with experiments and numerical simulations performed on the same cavity and presented in <cit.>. The bifurcation points extrapolated from the experiments and from the theoretical model both lie at $\Sigma=0.19$ with less than 1$\%$ relative deviation from each other, showing that the model and the experiments are consistent in their representation of the limit cycle's stability boundary.

§.§ Transmission without and with flow

If no flow is imposed on the cavity, then it behaves like a reciprocal, lossy and linear 3-port scatterer, as shown in  Fig. <ref>. This shows that, besides the whistling tone induced by the air flow, there are no non-linear effects present in the system.
If a flow is imposed on the cavity, a non-reciprocal lossless transmission can be achieved, as shown in Fig. <ref>.

cphv$R$=12.8 cm
dphv$D$=4 cm
ephv4.6 cm
fphv120$ ^{\circ}$
Bphv120$ ^{\circ}$
hphvair blower
iphvshop air
iphvshop air
jphvmass flow meter
Cphvmass flow meter
lphvcompression driver (source)
mphvacoustic foam
nphv70 cm
Aphv3 cm
Ephv3.2 cm

Lphvsemi-anechoic waveguides
Jphvcylindrical cavity
Mphvanechoic end  anechoic end
Kphvswirling air flow
Nphvfront view  back view
<phvside view

Sketch of the experimental setup used in this work (not for scale). A blower and shop air from the laboratory supply the wind channel of the experimental setup. The respective air flows are monitored using mass flow meters. At the swirler, air is injected tangentially to the wall of the circular cross-section wind channel, and it mixes with the main air from the blower. This configuration enables full control of the swirl intensity of the jet flowing through the cavity. This deep axisymmetric cavity is the central element of the 3-port circulator, which is equipped with three waveguides. Their angular positions are $120^\circ$ apart at the same radius. One of the waveguides is equipped with a compression driver to acoustically force one of the three circulator's ports. Microphones are placed along the waveguides and on the cavity for acoustic measurements. The wind channel has anechoic, catenoidal terminations which are filled with acoustic foam. Similarly, the ends of the waveguides attached to the cavity are also filled with anoechoic foam to achieve low acoustic reflection.

aphv2  3
vphv$s_\mathrm{in} $

mphv$\vert a_0\vert$
rphv$-\vert a_{\tiny 0}\vert$$0$$\vert a_0 \vert$
nphv$\vert a_\pm \vert$
sphv$\vert a_\pm \vert$
qphv $\gamma-\beta_-$
pphv $\gamma-\beta_+$
Uphv $a_\pm$
<phv$\vert s_\mathrm{in}\rangle$
>phv$\vert s_\mathrm{out}\rangle$

;phv$T_{j\shortrightarrow i}$

Some features of the theoretical model presented in the main text. Left inset: Typical damping non-linearity of the cavity modes considered in this work. In the absence of an incident wave, if the linear gain of the “$-$” mode, $\beta_-$, exceeds the linear damping $\gamma$, and if the non-linear component of the cavity gain decreases for increasing amplitude $\vert a_-\vert$, a stable limit cycle with amplitude $\vert a_0 \vert$ occurs, which corresponds to an equilibrium between gain and loss mechanisms. In this example, the mode $a_+$ is linearly stable because $\beta_+<\gamma$. Right inset: Sketch of the effective resistance associated with the non-linear cavity modes considered in this work. At large amplitudes $\vert a_\pm \vert$, the system exhibits a finite resistive term tending to $\gamma$, i.e., its response is essentially equivalent to that of a linearly stable resonance.

dphvac. pressure
ephv0  0.5  1
1.5  2
pphvbulk flow speed (m/s)
Iphvswirl number (-)
3phvwhistling intensity  (hPa2)
Fphv805 32
Ephv803.8 28.8
Dphv802.5 25.5
Kphv801.3 22.3
Hphv800 19
Gphv0  1/3  2/3
ophv$\Sigma$ (-)
Mphvfrequency $f_\pm$ (Hz)
Lphv$\beta_\pm$ (rad/s)
Sphv7.02    7.33    7.64   7.96     8.27
Tphv$V$ (m/s)
Uphv$U$ (m/s)
6phvisolation (dB)
>phvstable unstable
-phv$\tilde{\Omega}_\theta$ [s-1]

Bifurcation diagram of the circulator cavity. Shown is the whistling intensity, defined as the squared acoustic RMS amplitude $\overline{p}^2$, measured inside the cavity for different values of the bulk flow speed $\overline{U}=U+V$ and the swirl number $\Sigma=V/U$ (see Figure <ref> for a definition of $U$ and $V$). The operating point at $\overline{U}=23.1$ m/s and $\Sigma=0.49$ is marked by a red dot. The smaller insets show the variation of $\overline{p}^2$ (gray curves) along the cyan and green curves, respectively. The black lines in the smaller insets are linear regression lines to extrapolate the bifurcation points.

i[c] phv20
j[c] phv31.6
k[c] phv70
l[c] phv120
A[c]$f$ phv(kHz)
B[c]phvamp. $s_\mathrm{in}$phv(Pa)
C[c]$f$ phv(kHz)
phvPSD (dB)
E[c]$\Delta t$ phv(ms)
F[c]$p$ phv(Pa)
1[c]$\Delta t$

Synchronization dynamics of the nonlinear cavity mode. a Experimental Arnold tongue. Incident waves at fixed amplitudes  $s_{\text{in}}$ and frequencies $f$ were sent to the obstacle. By analyzing the relative peak height of the power spectrum density (PSD) of the forced mode, synchronized states can be identified. b A synchronized mode has an isolated peak at the forcing frequency. A non-synchronized mode has one peak at the self-sustained limit cycle frequency and another at the forcing frequency, as well as less dominant harmonics. c Pressure dynamics of a synchronized mode: quasi-sinusoidal oscillation at quasi-constant amplitude. d Non-synchronized states show beating oscillations due to interference of the self-sustained mode and the external forcing. The small insets in c and d show a couple of cycles of the same time traces.

dphvac. pressure
ephv0  0.5  1
1.5  2
pphvbulk flow speed $\overline{U}$ (m/s)
Iphvswirl number $\Sigma$ (-)
3phvwhistling intensity $\overline{p}^2$ (hPa2)
Fphv805 32
Ephv803.8 28.8
Dphv802.5 25.5
Kphv801.3 22.3
Hphv800 19
Gphv0  1/3  2/3
ophv$\Sigma$ (-)
Mphv $f_\pm$ (Hz)
Lphv$\beta_\pm$ (rad/s)
Sphv7.02    7.33    7.64   7.96     8.27
Tphv$V$ (m/s)
Uphv$U$ (m/s)
6phvisolation (dB)
>phvstable limit cycle
-phv$\tilde{\Omega}_\theta$ [s-1]

Modeled variation of the eigenfrequencies $f_\pm$, the linear gains $\beta_\pm$ and the decay rate $\gamma$ with respect to the swirl number $\Sigma$. The bias introduced by the swirl governs the Zeeman splitting of the eigenfrequencies that are degenerate at $\Sigma=0$. The values at the operating point are marked by red dots. The “$+$” and “$-$” branches correspond to the modes spinning with and against the swirling flow, respectively.

fphvfrequency (Hz)
lphv$ \Delta s $ (-)

Measured transmission coefficients in the no-flow case for different values of the incident wave amplitude $\Delta s$ (see Sec. <ref>) The black dotted curves correspond to $T_{1\shortrightarrow 3}$, the solid blue curves correspond to $T_{1\shortrightarrow 2}$ and the red lines correspond to 100$\%$ transmission.

iphvfrequency (Hz)
lphv$V$ (m/s)
kphv$U$ (m/s)

Sensitivity of the measured amplitude transmission coefficients at $\Delta s=9$ over different values of the axial ($U$) and azimuthal ($V$) flow components. The black dotted curves correspond to $T_{1\shortrightarrow 3}$, the solid blue curves correspond to $T_{1\shortrightarrow 2}$ and the red lines correspond to 100$\%$ transmission. The operating point is marked with a thick red box around the respective inset.

Supplementary notes to
“Synchronization-based lossless
non-reciprocal scattering”

Tiemo Pedergnana,$^{1}$ Abel Faure-Beaulieu$^{1}$, Romain Fleury$^{2*}$
and Nicolas Noiray$^{1*}$

${}^1$Department of Mechanical and Process Engineering, ETH
Zürich, Zürich, Switzerland

${}^2$Institute of Electrical and Micro Engineering, EPFL, Vaud, Switzerland

${}^*$Corresponding authors. E-mails<EMAIL_ADDRESS>
§ THEORY

§.§ Scattering matrix derivation

To model the non-linear scattering process, we propose a coupled-mode theory for the cavity modes $a=(a_+,a_-)^T$, which are forced by the incoming wave $\vert s_\mathrm{in}\rangle$, leading to an outgoing wave $\vert s_\mathrm{out}\rangle$. The following equations form the theory presented in this work:[For generality, we derive the scattering matrix for the general case where $\gamma_+$ and $\kappa_+$ are not necessarily equal to $\gamma_-$ and $\kappa_-$, respectively. In the main text, only the symmetric case with $\gamma_\pm=\gamma$, $\kappa_\pm=\kappa$ is considered. ]

\begin{eqnarray}
\dot{a}_\pm&=&i\omega_{_\pm}a_\pm  -\gamma_\pm a_\pm  +\frac{\beta_\pm a_\pm }{1+\kappa_\pm \vert a_\pm\vert^2}  +D^\dag_\pm \vert s_{\mathrm{in}}\rangle,\label{modal dynamics 2}\\
\vert s_{\mathrm{out}}\rangle&=&C \vert s_{\mathrm{in}}\rangle + D a, \label{input-output 2}
\end{eqnarray}

where a dot denotes the time derivative, $\omega_{\pm}$ are the eigenfrequencies, $\gamma_\pm$ are the decay rates corresponding to dissipative processes in the cavity and radiation losses, $\beta_\pm\in\mathbb{R}^+$ are the linear gain coefficients, $\kappa_\pm\in\mathbb{R}^+$ are the non-linearity constants, $C\in \mathbb{R}^{3\times 3}$ is the background scattering matrix and the matrix $D=[D_+,D_-]\in\mathbb{C}^{3\times 2}$ with column vectors $D_\pm$ describes the wave-mode coupling. When the system is harmonically excited from the $j^\mathrm{th}$ port, we have $\vert s_\mathrm{in}\rangle=s e^{i\omega t} \vert j\rangle$, where $s=\sqrt{\langle s_\mathrm{in} \vert  s_\mathrm{in} \rangle}\in\mathbb{R}^+$ is the incident wave amplitude and $\vert j \rangle$ is the unit vector in $j$-direction. The vector $a=(a_+,a_-)^T\in\mathbb{C}^2$ in Eq. (<ref>), where $(\cdot)^T$ is the transpose, comprises both cavity modes.

For simplicity, to compute the forced response, we write $\vert a_\pm \vert=A_\pm$ and we assume that the modes are perfectly synchronized with the excitation: $a_\pm=A_\pm e^{i(\omega t+\varphi_\pm)}$, where $\vert\cdot\vert\in \mathbb{R}^+$ is the absolute value and $A_\pm$ and $\varphi_\pm$ are the amplitudes and phases of the forced response, respectively. With this assumption, the modal dynamics (<ref>) become

\begin{eqnarray}
\big[\dot{A}_\pm+A_\pm i(\omega+\dot{\varphi}_\pm)\big]&=&\nonumber \\
&&\hspace{-2cm}\Big(i\omega_{0,_\pm} -\gamma_\pm +\frac{\beta_\pm}{1+\kappa_\pm A_\pm^2} \Big)A_\pm +\vert D_{\pm,j} \vert s e^{-i(\arg D_{\pm,j}+\varphi_\pm)}.
\end{eqnarray}

Separating real and imaginary parts yields

\begin{eqnarray}
\dot{A}_\pm&=&-\gamma_\pm A_\pm+\frac{\beta_\pm A_\pm}{1+\kappa_\pm A_\pm^2}+\vert D_{\pm,j} \vert s \cos{(\arg D_{\pm,j}+\varphi_\pm)}, \\
\dot{\varphi}_\pm&=&\Delta_\pm-\frac{\vert D_{\pm,j} \vert s\sin{(\arg D_{\pm,j}+\varphi_\pm)}}{ A_\pm},
\end{eqnarray}

where the terms $\Delta_\pm=\omega_{\pm}-\omega$ represent the detuning between the forcing and the respective oscillators. The unforced ($s=0$) limit cycles $A_\pm=A_{0,\pm}\in\mathbb{R}^+$, $\varphi_\pm=\varphi_{0,\pm}\in[0, 2\pi)$ are given by

\begin{eqnarray}
A_{0,\pm}&=&\sqrt{\frac{1}{\kappa_\pm}\Big(\frac{\beta_\pm}{\gamma_\pm}-1\Big)}, \label{Limit cycle formula}
\end{eqnarray}

and can only exist when the respective linear gain is greater than the dissipation: $\beta_\pm>\gamma_\pm$. Assuming this is the case, the linearized unforced modal dynamics around the limit cycles $A_{0,\pm}$,

\begin{eqnarray}
\dot{A}'_\pm&=&-\Bigg(\gamma_\pm+\frac{\beta_\pm (\kappa_\pm A_{0,\pm}^2-1)}{(1+\kappa_\pm A_{0,\pm}^2)^2}\Bigg)A'_\pm\nonumber\\
\end{eqnarray}

are asymptotically stable. If $\beta_\pm<\gamma_\pm$, the limit cycles $A_{0,\pm}$ given by Eq. (<ref>) do not exist and the only stationary solution is the trivial one: $A_\pm\equiv 0$. In that case, the linearized dynamics around the trivial solution,

\begin{eqnarray}
\dot{A}'_\pm&=&-(\gamma_\pm-\beta_\pm)A'_\pm, \label{Trivial solution stability}
\end{eqnarray}

are also asymptotically stable. If there does exist a limit cycle, then we see from Eq. (<ref>) that the trivial solution is linearly unstable. Since the unforced modal equations (<ref>) are decoupled, depending on the sign of $\beta_\pm-\gamma_\pm$, each of the modes $a_\pm$ has exactly one linearly stable stationary solution, which can be either a limit cycle or the trivial solution.

Having discussed the modes' unforced dynamics, we now compute their forced response by setting $\dot{A}_\pm=\dot{\varphi}_\pm=0$, leading to the following expression:

\begin{eqnarray}
\bigg[\Big(\vert D_{\pm,j} \vert s\Big)^2-(\Delta_\pm^2+\gamma_\pm^2)A_\pm^2\bigg](1+\kappa_\pm A_\pm^2)^2+2\gamma_\pm \beta_\pm A_\pm^2(1+\kappa_\pm A_\pm^2)\nonumber\\
-\beta_\pm^2 A_\pm^2=0,
\end{eqnarray}

which can be written as two cubic equations for $A_\pm^2$:

\begin{eqnarray}
0&=&\alpha_{0,\pm}+\alpha_{1,\pm} A_\pm^2 +\alpha_{2,\pm} A_\pm^4 +\alpha_{3,\pm} A_\pm^6, \label{forced response A}
\end{eqnarray}

\begin{eqnarray}
\alpha_{0,\pm}&=&\vert  D_{\pm,j} \vert^2 s^2\nonumber\\
\alpha_{1,\pm}&=&2 \vert D_{\pm,j} \vert^2 s^2 \kappa_\pm-(\Delta_\pm^2+\gamma_\pm^2)+\beta_\pm(2\gamma_\pm-\beta_\pm)\nonumber\\
\alpha_{2,\pm}&=&\kappa_\pm^2\vert D_{\pm,j} \vert^2 s^2-2\kappa_\pm(\Delta_\pm^2+{\gamma_\pm^2})+2\gamma_\pm\beta_\pm\kappa_\pm\nonumber\\
\alpha_{3,\pm}&=&-\kappa_\pm^2(\Delta_\pm^2+{\gamma_\pm^2})\nonumber
\end{eqnarray}

The linearly stable branches of $\varphi_\pm$ of the forced response are given by

\begin{eqnarray}
\varphi_\pm=-\arg D_{\pm,j} +\arcsin\Big(\frac{\Delta_\pm A_\pm}{\vert D_{\pm,j} \vert s}\Big). \label{forced response phi}
\end{eqnarray}

At infinitely large amplitudes $A_\pm$, the term $-\gamma_\pm A_\pm$ dominates the $\beta_\pm$-modulated term in the equations for $\dot{A}_\pm$. In this limit, therefore, the dynamics (<ref>) converge to

\begin{eqnarray}
\dot{a}_\pm&=&\Big(i\omega_\pm -\gamma_\pm  \Big)a_\pm +D^\dag \vert s_{\mathrm{in}}\rangle,\label{modal dynamics 2, large A}\\
\vert s_{\mathrm{out}}\rangle&=&C \vert s_{\mathrm{in}}\rangle + D a, \label{input-output 2, large A}
\end{eqnarray}

which are exactly equivalent to the classic temporal coupled-mode theory[Fan, S., Suh, W. & Joannopoulos, J. Temporal coupled-mode theory for the fano resonance in optical resonators. J. Opt. Soc. Am. A 20 (3), 569–572 (2003). <https://doi.org/10.1364/JOSAA.20.000569>;
Wang, Z. & Fan, S. Magneto-optical defects in two-dimensional photonic crystals. Appl. Phys. B 81, 369–375 (2005). <https://doi.org/10.1007/ s00340-005-1846-x>]. To model this limiting case in a manner which is consistent with non-reciprocal scattering, following the method from Pedergnana and Noiray[Pedergnana, T. & Noiray, N. Superradiant scattering by a limit cycle. Phys. Rev. Appl. 20, 034068 (2023). <https://doi.org/10.1103/PhysRevApplied. 20.034068>.], we superimpose a non-reciprocal optimal scattering matrix $S_*$ on a purely reflectional background $C$:

\begin{eqnarray}
0 & 1 & 0\\
0 & 0 & 1\\
1 & 0 & 0
\end{pmatrix},\quad
1 & 0 & 0\\
0 & 1 & 0\\
0 & 0 & 1
\end{pmatrix}.
\end{eqnarray}

The matrix $\sigma S_*-C$, where $\sigma\in[0,1]$ is the unitarity factor, has the following eigenvalues $\mu_i$ and corresponding eigenvectors $v_i$, $i=1,2,3$:

\begin{align}
\mu_1&=\sigma-1,\quad &v_1&=\frac{1}{\sqrt{3}}\begin{pmatrix}
1 \\
\end{pmatrix},\\
\mu_2&=\sigma e^{\frac{i\pi }{3}}+1,\quad &v_2&=\frac{1}{\sqrt{3}}\begin{pmatrix}
e^{-\frac{2\pi i}{3}}\\
e^{\frac{2\pi i}{3}}\\
\end{pmatrix},\\
\mu_3&=\sigma e^{\frac{-i\pi }{3}}+1=\mu_2^*,\quad &v_3&=\frac{1}{\sqrt{3}}\begin{pmatrix}
e^{\frac{2\pi i}{3}}\\
e^{\frac{-2\pi i}{3}}\\
\end{pmatrix}=v_2^*,
\end{align}

where $(\cdot)^*$ is the complex conjugate. Formally approximating the linear part of the scattering matrix by $\sigma S_*$ yields a low-rank approximation problem for $\sigma S_*-C$ by the coupling matrix $D$:

\begin{eqnarray}
\sigma S_*-C\approx D F^{-1} D^\dag, \label{Low-rank approximation}
\end{eqnarray}

\begin{eqnarray}
i(\omega-\omega_+)+\gamma_+& 0\\
0 & i(\omega-\omega_-)+\gamma_-
\end{pmatrix}. \label{Definition F}
\end{eqnarray}

In the recent work from Pedergnana and Noiray[Pedergnana, T. & Noiray, N. Superradiant scattering by a limit cycle. Phys. Rev. Appl. 20, 034068 (2023). <https://doi.org/10.1103/PhysRevApplied. 20.034068>.], consistency conditions are given to determine $D$ and perfect matching condition for a scatterer with one mode and two ports. The latter condition describes an optimal parameter configuration at which the dominant eigenspace of $\sigma S_*-C$ is best approximated by $D F^{-1} D^\dag$. The goal of the consistency conditions is to ensure that the linear scattering matrix derived in this section successfully approximates the dominant sub-spectrum of $\sigma S_*$. Because $D$ is the only unknown in our formulation of the wave-mode coupling, we consider the problem (<ref>) instead of directly approximating $\sigma S_*$ by $C+D F^{-1} D^\dag$. We formulate the following requirements for a scatterer with two modes and three ports:

(I) The singular (rank $2$) matrix $DF^{-1}D^\dag$ spans the subspace corresponding to the eigenvalues of $\sigma S_*-C$ with the largest absolute values $\vert \mu_{i} \vert$, $i=1,2,3$. In general, $DF^{-1}D^\dag$ does not share any eigenvalues with $\sigma S_*-C$.

(II) At perfect matching for a non-unitary process $\sigma<1$, the two out of three eigenvalues with largest absolute values and the corresponding eigenvectors of $DF^{-1}D^\dag$ coincide with the two of $\sigma S_*-C$ with largest absolute value.

We show that (R1) and (R2) imply that at perfect matching, $\sigma S_*$ and $C+DF^{-1}D^\dag$ also share the sub-spectrum corresponding to the two eigenvalues with largest absolute value. To understand how $D$ is determined by (R1) and (R2), we write out the right-hand-side of the approximation (<ref>) explicitly:

\begin{eqnarray}
D_+ & D_-
\end{pmatrix}\begin{pmatrix}
\dfrac{1}{g_+ \gamma_+} & 0\\
0 & \dfrac{1}{g_- \gamma_-}
\end{pmatrix}\begin{pmatrix}
D_+^\dag \\
\end{pmatrix}\\
&=&\dfrac{D_+ D_+^\dag}{g_+\gamma_+}+\dfrac{D_- D_-^\dag}{g_-\gamma_-}, \label{RHS of approximation}
\end{eqnarray}

where we have defined $g_\pm \gamma_\pm=i(\omega-\omega_\pm)+\gamma_\pm$. In comparison, the left-hand-side of the approximation (<ref>) can be expressed as

\begin{eqnarray}
\sigma S_*-C=v_1 v_1^\dag \mu_1+v_2 v_2^\dag \mu_2+v_3 v_3^\dag \mu_3. \label{relative s matrix eigendecomposition}
\end{eqnarray}

To approximate the dominant part of the spectrum of $\sigma S_*-C$, noting that $\vert\mu_{2}\vert=\vert\mu_{3}\vert>\vert\mu_1\vert$, we make the Ansatz

\begin{eqnarray}
D_+=\sqrt{\gamma_+\alpha(\sigma)}\, v_{2},\quad  D_-=\sqrt{\gamma_-\alpha(\sigma)}\, v_{3}, \label{D ansatz}
\end{eqnarray}

where $\alpha\in\mathbb{R}^+$. In the unitary case ($\sigma=1$), $\alpha=2$ and the perfect matching condition is $\omega_\pm=\omega\pm\gamma_\pm/\sqrt{3}$ [Wang, Z. & Fan, S. Magneto-optical defects in two-dimensional photonic crystals. Appl. Phys. B 81, 369–375 (2005). <https://doi.org/10.1007/ s00340-005-1846-x>]. By analogy, we assume the perfect matching condition for $\sigma<1$ has the form

\begin{eqnarray}
\omega_\pm=\omega\pm\chi(\sigma)\gamma_\pm, \label{perfect matching ansatz}
\end{eqnarray}

where $\chi\in\mathbb{R}^+$. Substituting the assumed relations (<ref>) and (<ref>) into Eq. (<ref>) yields, at perfect matching,

\begin{eqnarray}
DF^{-1}D^\dag&=&\alpha(\sigma)\bigg[\dfrac{v_2 v_2^\dag}{1-i\chi(\sigma)}+\dfrac{v_3 v_3^\dag}{1+i\chi(\sigma)}\bigg]. \label{RHS of approximation after substitution}
\end{eqnarray}

Comparing Eq. (<ref>) to Eq. (<ref>), we require that

\begin{eqnarray}
\dfrac{\alpha(\sigma)}{1-i\chi(\sigma)}=\mu_2\quad\leftrightarrow\quad \dfrac{\alpha(\sigma)}{1+i\chi(\sigma)}=\mu_3. \label{int. res. 1}
\end{eqnarray}

Since $\alpha\in\mathbb{R}$, $\mu_{2,3}=1+\sigma(1\pm i\sqrt{3})/2$ and $\vert \mu_{2,3}\vert^2 =\sigma^2+\sigma+1$, Eq. (<ref>) has the following unique solutions for $\alpha$ and $\chi$:

\begin{eqnarray}
\alpha&=&\dfrac{2(\sigma^2+\sigma+1)}{\sigma+2},\\
\chi&=&\frac{\sqrt{3}\sigma}{\sigma+2},
\end{eqnarray}

which, for $\sigma= 1$, reduce to the well-known values $\alpha=2$ and $\chi=1/\sqrt{3}$ found in literature on unitary scattering. Using the above results, the coupling matrix $D$ becomes

\begin{eqnarray}
\sqrt{\gamma_+} e^{-\frac{2\pi i}{3}} &\sqrt{\gamma_-}e^{\frac{2\pi i}{3}} \\
\sqrt{\gamma_+} e^{\frac{2\pi i}{3}}&\sqrt{\gamma_-}e^{-\frac{2\pi i}{3}}  \\
\sqrt{\gamma_+}&\sqrt{\gamma_-}
\end{pmatrix}. \label{coupling matrix D}
\end{eqnarray}

By computing the linear scattering matrix $C+DF^{-1}D^\dag$, it is straightforward to verify that for a perfectly matched sub-unitary scatterer ($\sigma<1$), its two eigenvalues with largest absolute value coincide with those of $\sigma S_*$. Indeed, at perfect matching,

\begin{eqnarray}
\sigma-1 &-2\sigma-1&\sigma-1 \\
\sigma-1&\sigma-1& -2\sigma-1\\
\end{pmatrix},
\end{eqnarray}

which has the eigenvalues $\lambda_1=-1$, $\lambda_{2,3}=\mu_{2,3}-1$ and corresponding eigenvectors $v_{i}$, $i=1,2,3$. The optimal scattering matrix $\sigma S_*$ also shares the same eigenvectors, with eigenvalues $\zeta_1=-\sigma$ and $\zeta_{2,3}=\lambda_{2,3}$. Therefore, for a perfectly matched unitary scatterer with $\sigma=1$, the optimal state is reached exactly: $C+DF^{-1}D^\dag=S_*$. This shows that the linear scattering matrix derived in this section successfully approximates the dominant sub-spectrum of $\sigma S_*$.

Using the above intermediate result and the Moore-Penrose Pseudoinverse[Pedergnana, T. & Noiray, N. Superradiant scattering by a limit cycle. Phys. Rev. Appl. 20, 034068 (2023). <https://doi.org/10.1103/PhysRevApplied. 20.034068>.], the non-linear scattering matrix is found to be

\begin{eqnarray}
S&=& \dfrac{\vert {s}_\mathrm{out}\rangle \langle s_\mathrm{in}\vert}{\langle s_\mathrm{in}\vert s_\mathrm{in}\rangle} \label{general scattering matrix}\\
A_{j,+} e^{i \varphi_{j,+}} \\
A_{j,-} e^{i \varphi_{j,-}}
\end{pmatrix} \langle j\vert, \label{Explicit scattering matrix}
\end{eqnarray}

where $A_{j,+}$ and $\varphi_{j,+}$, $j=1,2,3$, are obtained by solving the forced response equations (<ref>) and (<ref>) for $j=1,2,3$.

At large amplitudes $A_\pm >> A_{0,\pm}$, the non-reciprocal amplitude transmission coefficients $T_{j\shortrightarrow i}=\vert (C+DF^{-1}D^\dag)_{ij}\vert$ converge to

\begin{eqnarray}
\hspace{-0.31cm}&&T_{1\shortrightarrow 2}\approx\bigg\vert\dfrac{2(\sigma^2 + \sigma + 1)}{3(\sigma + 2)}\bigg( \dfrac{e^{-i 2\pi/3}}{1+i(\omega-\omega_+)/\gamma_+}+\dfrac{e^{i 2\pi/3}}{1+i(\omega-\omega_-)/\gamma_-} \bigg)\bigg\vert, \label{Trans12} \\
\hspace{-0.31cm}&&T_{1\shortrightarrow 3}\approx\bigg\vert\dfrac{2(\sigma^2 + \sigma + 1)}{3(\sigma + 2)}\bigg( \dfrac{e^{i 2\pi/3}}{1+i(\omega-\omega_+)/\gamma_+}+\dfrac{e^{-i 2\pi/3}}{1+i(\omega-\omega_-)/\gamma_-} \bigg)\bigg\vert. \label{Trans13}
\end{eqnarray}

Perfect matching is achieved for $\omega_\pm=\omega\pm\gamma_\pm\sqrt{3}\sigma/(\sigma+2)$, leading to $T_{1\shortrightarrow 3}=(2\sigma+1)/3$ and $T_{1\shortrightarrow 2}=\vert 1-\sigma\vert/3$.

§.§ Large amplitude limit

Equations (<ref>) and (<ref>) describe the asymptotically linear scattering process which occurs in the limit of infinitely large incident wave amplitudes. To understand how the non-linear terms in the modal dynamics (<ref>) affect the transmission coefficients in the large (but not infinite) amplitude limit, we use a perturbation series approach. For this, we introduce the small, auxiliary parameter  $\varepsilon > 0$. In the following, we assume forcing from port  “1”: $\vert s_\mathrm{in}\rangle=s e^{i \omega t}\vert 1 \rangle$, where the incident wave amplitude $s$ is assumed to be of order $\varepsilon^{-1}$, i.e., large. For $\vert a \vert \rightarrow \infty$, the coupled-mode equations (<ref>) and (<ref>) reduce to Eqs. (<ref>) and (<ref>). These limiting cases predict that both $a$ and $\vert s_\mathrm{out}\rangle$ depend linearly on $\vert s_\mathrm{in}\rangle$ in the large amplitude limit. Therefore, we assume that the leading order terms in $a$ and $\vert s_\mathrm{out}\rangle$ scale also with $s$, or $1/\varepsilon$, in the large amplitude limit. We later set $\varepsilon$ equal to $1$ and verify analytically that the assumed scaling ($a\propto s$) is indeed reproduced by our perturbation approach, applied to the non-linear coupled-mode theory described the previous section. To compute the perturbation series, we define the rescaled quantities

\begin{eqnarray}
\tilde{a}_\pm&=& \varepsilon a_\pm, \label{rescaled amplitude}\\
\vert \tilde{s}_\mathrm{in}\rangle&=& \varepsilon \vert {s}_\mathrm{in}\rangle, \label{rescaled incident wave}\\
\vert \tilde{s}_\mathrm{out}\rangle&=& \varepsilon \vert {s}_\mathrm{out}\rangle, \label{rescaled outgoing wave}
\end{eqnarray}

which, by the above assumptions, are all of order $O(1)$. Substituting Eqs. (<ref>)-(<ref>) into the non-linear coupled-mode theory (<ref>) and (<ref>) yields, after cancelling out $\varepsilon$,

\begin{eqnarray}
{\dot{\tilde{a}}_{\pm}}&=&{i\omega_{\pm} \tilde{a}_{\pm}} -{\gamma_\pm \tilde{a}_{\pm}}+\dfrac{\varepsilon^2\beta_\pm \tilde{a}_\pm }{\varepsilon^2+{\kappa_\pm \vert \tilde{a}_\pm\vert^2}}  +{D^\dag_\pm \vert \tilde{s}_{\mathrm{in}}\rangle},\label{expanded modal dynamics 2}\\
\vert \tilde{s}_\mathrm{out}\rangle&=&C \vert \tilde{s}_{\mathrm{in}}\rangle + D \tilde{a}. \label{expanded input-output 2}
\end{eqnarray}

Since $\varepsilon$ was assumed to be small, the rational function in Eq. (<ref>) can be expanded as a geometric series:

\begin{eqnarray}
\dfrac{\varepsilon^2\beta_\pm \tilde{a}_\pm }{\varepsilon^2+{\kappa_\pm \vert \tilde{a}_\pm\vert^2}} &=&\dfrac{\varepsilon^2\beta_\pm \tilde{a}_\pm}{\kappa_\pm \vert \tilde{a}_\pm\vert^2}\Big( 1-\frac{\varepsilon^2}{\kappa_\pm \vert \tilde{a}_\pm\vert^2}+\dots \Big). \label{geo series}
\end{eqnarray}

Next, we expand $\tilde{a}_\pm$ and the outgoing wave in terms of the perturbation parameter:

\begin{eqnarray}
\tilde{a}_\pm&=& \tilde{a}_{\pm,0}+\varepsilon^2 \tilde{a}_{\pm,1}+O(\varepsilon^4), \label{expansion amplitude}\\
\vert {s}_\mathrm{out}\rangle&=&  \vert \tilde{s}_\mathrm{out,0}\rangle+\varepsilon^2  \vert \tilde{s}_\mathrm{out,1}\rangle+O(\varepsilon^4). \label{expansion outgoing wave}
\end{eqnarray}

Using the general expression (<ref>), the expansion (<ref>) also implies that the scattering matrix can be expanded analogously:

\begin{eqnarray}
S &=& \dfrac{\big[\vert \tilde{s}_\mathrm{out,0}\rangle+\varepsilon^2  \vert \tilde{s}_\mathrm{out,1}\rangle+O(\varepsilon^4)\big]\langle s_\mathrm{in}\vert}{\langle s_\mathrm{in}\vert s_\mathrm{in}\rangle} \nonumber\\
&=&S_0+\varepsilon^2  S_1+O(\varepsilon^4)\,. \label{expansion S matrix}
\end{eqnarray}

Because the modal dynamics (<ref>) contain a quadratic non-linearity, only even powers of $\varepsilon$ appear in the geometric series (<ref>). Accordingly, only even orders of $\varepsilon$ contribute to the expansions (<ref>) and (<ref>).

By substituting Eqs. (<ref>)-(<ref>) into the rescaled system (<ref>) and (<ref>), we obtain the contributions of different orders in $\varepsilon$ to the scattering matrix. To zeroth order $\varepsilon=0$, as expected, this procedure yields the modal equations in the limit of infinitely large amplitudes:

\begin{eqnarray}
{\dot{\tilde{a}}_{\pm,0}}&=&{i\omega_{\pm} \tilde{a}_{\pm,0}} -{\gamma_\pm \tilde{a}_{\pm,0}}+{D^\dag_\pm \vert \tilde{s}_{\mathrm{in}}\rangle},\label{expanded modal dynamics 3}\\
\vert \tilde{s}_\mathrm{out,0}\rangle&=&C \vert \tilde{s}_{\mathrm{in}}\rangle + D \tilde{a}_{0}. \label{expanded input-output 3}
\end{eqnarray}

As standard in the TCMT literature, we assume synchronized conditions in the steady state, $\dot{\tilde{a}}_\pm=i\omega \tilde{a}_\pm$, which leads to

\begin{eqnarray}
\tilde{a}_{\pm,0}=\dfrac{D^\dag_\pm \vert \tilde{s}_\mathrm{in}\rangle}{i(\omega-\omega_\pm)+\gamma_\pm}. \label{Solution a0}
\end{eqnarray}

Therefore, the leading-order contribution to the scattering matrix at large amplitudes is the TCMT scattering matrix

\begin{eqnarray}
\end{eqnarray}

where the matrix $F$ is defined in Eq. (<ref>). Following the above steps and using the expansion (<ref>) we obtain, to order $O(\varepsilon^2)$,

\begin{eqnarray}
{\dot{\tilde{a}}_{\pm,1}}&=&{i\omega_{\pm} \tilde{a}_{\pm,1}} -{\gamma_\pm \tilde{a}_{\pm,1}}+\dfrac{\beta_\pm \tilde{a}_{\pm,0} }{{\kappa_\pm \vert \tilde{a}_{\pm,0}\vert^2}} ,\label{expanded modal dynamics 4}\\
\vert \tilde{s}_\mathrm{out,1}\rangle&=&D \tilde{a}_{1}. \label{expanded input-output 4}
\end{eqnarray}

Substituting Eq. (<ref>) into Eq. (<ref>) and solving for $\tilde{a}_{\pm,1}$ yields

\begin{eqnarray}
\tilde{a}_{\pm,1}=\dfrac{(\omega-\omega_\pm)^2+\gamma_\pm^2}{\kappa_\pm \langle \tilde{s}_\mathrm{in}\vert D_\pm D^\dag_\pm \vert \tilde{s}_\mathrm{in}\rangle}\dfrac{\beta_\pm D^\dag_\pm \vert \tilde{s}_\mathrm{in}\rangle}{[i(\omega-\omega_\pm)+\gamma_\pm]^2}. \label{Solution a1}
\end{eqnarray}

Using Eq. (<ref>) and setting $\varepsilon=1$, then, the next-to-leading order term in the scattering matrix expansion (<ref>) is found to be

\begin{eqnarray}
S_1=DF_1^{-1}(s) D^\dag, \label{order 1 S contribution}
\end{eqnarray}

where $s=\sqrt{\langle{s}_\mathrm{in}\vert {s}_\mathrm{in}\rangle}$ is the incident wave amplitude and

\begin{eqnarray}
\theta_+(s) & 0\\
0 & \theta_-(s)
\end{pmatrix}, \label{Definition F1}\\
\theta_\pm(s)&=&\dfrac{(\omega-\omega_\pm)^2+\gamma_\pm^2}{\kappa_\pm \langle {s}_\mathrm{in}\vert D_\pm D^\dag_\pm \vert {s}_\mathrm{in}\rangle}\dfrac{\beta_\pm}{[i(\omega-\omega_\pm)+\gamma_\pm]^2}.
\end{eqnarray}

As assumed above in the expansion (<ref>), $S_1$ is of order $O(\varepsilon^2)=O(s^{-2})$. For infinitely large incident wave amplitudes $s\rightarrow \infty$, $S_1$ and all higher expansion terms vanish. One could continue the perturbation series to show that the next-to-next-to-leading order term is of order $O(s^{-4})$ etc. Therefore, for large $s$, the term (<ref>) will dominate the higher-order terms in the expansion (<ref>), justifying our perturbative approach. Accordingly, all higher-order contributions to the scattering matrix in the large amplitude limit are neglected in the following.

Collecting the above results, we obtain the first-order correction to the scattering matrix due to the non-linearity in the modal dynamics (<ref>):

\begin{eqnarray}
S= C+D \big[F^{-1}+F_1^{-1}(s) \big]D^\dag + O(s^{-4}).
\end{eqnarray}

To obtain a simple analytical expression for the scattering matrix at optimally tuned conditions, we assume perfectly matched frequencies: $\omega_\pm=\omega\pm\gamma_\pm\sqrt{3}\sigma/(\sigma+2)$. As in the numerical examples shown in the main text, we set $\kappa_\pm=\kappa$. These simplifications lead to

\begin{eqnarray}
S_{21}&=& \dfrac{\sigma-1}{3} -\dfrac{(3-6\sigma-6\sigma^2)(\beta_++\beta_-)- 3\sqrt{3}i(1+2\sigma)(\beta_- - \beta_+)}{6 \kappa s^2 (\sigma^2 + \sigma + 1)} \nonumber \\ &&+O(s^{-4}), \label{T12lim} \\
S_{31}&=&-\dfrac{ (2\sigma+1)  }{3}-\dfrac{(3+12\sigma+3\sigma^2)(\beta_++\beta_-) + 3\sqrt{3}i(1-\sigma^2)(\beta_- - \beta_+) }{6 \kappa s^2 (\sigma^2 + \sigma + 1)}\nonumber\\ &&+O(s^{-4}). \label{T13lim}
\end{eqnarray}

We can now rewrite the approximate amplitude transmission coefficients as follows:

\begin{eqnarray}
T_{1\shortrightarrow 2}&=& \sqrt{\bigg[\dfrac{1-\sigma}{3} -\dfrac{(6\sigma^2+6\sigma-3)(\beta_++\beta_-)}{6 \kappa s^2 (\sigma^2 + \sigma + 1)}\bigg]^2+\bigg[\dfrac{3\sqrt{3}(1+2\sigma)(\beta_- - \beta_+)}{6 \kappa s^2 (\sigma^2 + \sigma + 1)}\bigg]^2} \nonumber\\ &&+O(s^{-4}), \label{T12lim3} \\
T_{1\shortrightarrow 3}&=&\sqrt{\bigg[\dfrac{ (2\sigma+1)  }{3}+\dfrac{(3+12\sigma+3\sigma^2)(\beta_++\beta_-)}{6 \kappa s^2 (\sigma^2 + \sigma + 1)}\bigg]^2 + \bigg[\dfrac{3\sqrt{3} (1-\sigma^2)(\beta_- - \beta_+) }{6 \kappa s^2 (\sigma^2 + \sigma + 1)}\bigg]^2}\nonumber\\ &&+O(s^{-4}). \label{T13lim3}
\end{eqnarray}

In the following, we assume that the second brackets under the square roots in Eqs. (<ref>) and (<ref>), respectively, is small, which is generally true for small asymmetry $\beta_+ \approx \beta_-$. This leads to

\begin{eqnarray}
T_{1\shortrightarrow 2}&\approx& \bigg \vert \dfrac{1-\sigma}{3} -\dfrac{(2\sigma^2+2\sigma-1)(\beta_++\beta_-)}{2 \kappa s^2 (\sigma^2 + \sigma + 1)}\bigg \vert+O(s^{-4}),\nonumber\\ \label{T12lim4} \\
T_{1\shortrightarrow 3}&\approx&\bigg \vert \dfrac{ (2\sigma+1)  }{3}+\dfrac{(1+4\sigma+\sigma^2)(\beta_++\beta_-)}{2 \kappa s^2 (\sigma^2 + \sigma + 1)}\bigg \vert+O(s^{-4}).\nonumber\\ \label{T13lim4}
\end{eqnarray}

We also assume that

* $\sigma>(\sqrt{3}-1)/2\approx 0.366$, so that the term $2\sigma^2+2\sigma-1$ in Eq. (<ref>) is positive.

* That the $s^{-2}$-modulated term in Eq. (<ref>) is less than or equal to $(1-\sigma)/3$, such that the  respective expression inside the vertical lines is non-negative.

The first of these assumptions is a minor restriction, since values of $\sigma \in [0.4, 1]$ represent a broad range of dissipative conditions (see also the next section for a quantitative discussion of this point). The second of the above assumptions is generally true for large enough $s$. We note that, in the parameter range considered, $T_{1\shortrightarrow 2}$ is decreased and $  T_{1\shortrightarrow 3}$ is increased by the non-linear effects described the modal dynamics (<ref>) for $\beta_\pm >0$. For dissipative conditions ($\sigma<1$), the non-linear interference with the limit cycle compensates–at least partially–for the dissipation incurred at port “3”, while simultaneously decreasing transmission to port “2”. To summarize, the above results demonstrate analytically that the non-linear wave-mode coupling envisioned in this work enhances the non-reciprocity of the 3-port scatterer at large incident wave amplitudes, while simultaneously decreasing transmission losses.

§.§ Unitarity and reversiblity

We model dissipation, or irreversible losses, in the infinite amplitude limit (<ref>) and (<ref>) using the intrinsic decay rates $\gamma_{\mathrm{i},\pm}\in\mathbb{R}^+$. Specifically, we assume that the decay rates $\gamma_\pm$ can be decomposed into a part $\gamma_{\mathrm{r},\pm}$ related to (reversible) decay through the ports (radiation losses) and another part $\gamma_{\mathrm{i},\pm}$ representing irreversible losses:

\begin{eqnarray}
\gamma_\pm=\gamma_{\mathrm{r},\pm}+\gamma_{\mathrm{i},\pm}.
\end{eqnarray}

The decay through the ports is reversible because the wave energy transmitted through the ports is available to the environment, while the irreversible losses are irrecoverable. In our theory, the presence of intrinsic losses is marked by a broken time-symmetry: while in forward time, the perturbations $a'_\pm$ decay like $\exp{(-\gamma_\pm t)}$, their time-reversed counterparts $(a'_\pm)^*$ grow like $\exp{(\gamma_{\mathrm{r},\pm} t)}$ in backward time, i.e., only the reversible part of the decay rate $\gamma_{\mathrm{r},\pm}$ contributes to the growth of time-reversed perturbations. Following the work of Zhao, Guo and Fan[Zhao, Z., Guo, C. & Fan, S. Connection of temporal coupled-mode-theory formalisms for a resonant optical system and its time-reversal conjugate. Phys. Rev. A 99 (3) (2019). <https://doi.org/10.1103/PhysRevA.99.033839>], the balance of reversible energy is expressed as

\begin{eqnarray}
D^\dag D&=&2\begin{pmatrix}
\gamma_{\mathrm{r},+}&0 \\
\end{pmatrix} \nonumber\\
\gamma_+-\gamma_{\mathrm{i},+}&0 \\
\end{pmatrix}.\label{Definition internal decay rates}
\end{eqnarray}

To confirm the broken time symmetry, we apply a time-reversal transformation to the exponential decay of the perturbations $a'$ in an initially energized cavity, which is governed by Eq. (<ref>), in the absence of forcing $\vert s_\mathrm{in}\rangle=0$. This transformation interchanges $\vert s_\mathrm{in}\rangle$ and $\vert s_\mathrm{out}\rangle$ in Eq. (<ref>). In the time-reversed scenario, time-reversed perturbations $(a')^*$ grow exponentially in backward time with a priori undetermined growth rates $\tilde{\gamma}_\pm\in\mathbb{R}^+$. The time-reversed dynamics are therefore given by

\begin{eqnarray}
\dot{a}'^*(-t)=(i\omega_0 +\tilde{\gamma}_\pm)a'^*(-t) -D_\pm^T \vert s_\mathrm{out}\rangle^*(-t), \label{time-reversed linearized modal dynamics}
\end{eqnarray}

where $(\cdot)^T$ is the transpose. At the complex frequencies $\omega=\omega_0-i\tilde{\gamma}_\pm$, we have

\begin{eqnarray}
\tilde{\gamma}_+&0 \\
\end{pmatrix} a'^*&=& D^T \vert s_\mathrm{out}\rangle^*\\
&=& D^T(D^*a'^*)\\
\Rightarrow D^\dag D &=&2\begin{pmatrix}
\tilde{\gamma}_+&0 \\
\end{pmatrix}. \label{auxiliary eq. 1}
\end{eqnarray}

Comparing Eqs. (<ref>) and (<ref>) yields

\begin{eqnarray}
\tilde{\gamma}_\pm=\gamma_\pm-\gamma_{\mathrm{i},\pm}=\gamma_{\mathrm{r},\pm},
\end{eqnarray}

demonstrating that time-symmetry is broken ($\gamma_\pm\neq \tilde{\gamma}_\pm$) for $\gamma_{\mathrm{i},\pm}\neq 0$. Using Eqs. (<ref>) and (<ref>), a relation between the intrinsic decay rates $\gamma_{\mathrm{i},\pm}$ and the unitarity factor $\sigma$ is established:

\begin{eqnarray}
\dfrac{\gamma_{\mathrm{i},\pm}}{\gamma_\pm}=\dfrac{1-\sigma^2}{\sigma+2}\geq0. \label{relation irrev. and sigma}
\end{eqnarray}

The non-unitarity and the irreversible losses of the scattering process are intrinsically connected. In particular, unitary scattering ($\sigma=1$) corresponds to perfect reversibility ($\gamma_{\pm}=\gamma_{\mathrm{r},\pm}$). It was mentioned in the previous section that $1>\sigma>0.4$ describes a broad range of dissipative conditions. We see from Eq. (<ref>) that, at $\sigma=0.4$, the irreversible decay rates $\gamma_\mathrm{i,\pm}=0.35 \gamma_\pm$ make up more than a third of the total decay rates $\gamma_\pm$ in the limit of infinitely large amplitudes.