Graph-based analysis of nonreciprocity in coupled-mode systems

In this work we derive the general conditions for obtaining nonreciprocity in multi-mode parametrically-coupled systems. The results can be applied to a broad variety of optical, microwave, and hybrid systems including recent electro- and opto-mechanical devices. In deriving these results, we use a graph-based methodology to derive the scattering matrix. This approach naturally expresses the terms in the scattering coefficients as separate graphs corresponding to distinct coupling paths between modes such that it is evident that nonreciprocity arises as a consequence of multi-path interference and dissipation in key ancillary modes. These concepts facilitate the construction of new devices in which several other characteristics might also be simultaneously optimized. As an example, we synthesize a novel three-mode unilateral amplifier design by use of graphs. Finally, we analyze the isolation generated in a common parametric multi-mode system, the DC-SQUID.


INTRODUCTION
Reciprocity is the symmetry of a physical system with respect to the exchange of source and detector. In devices that demonstrate nonreciprocal behaviour, such as electromagnetic isolators and circulators, the transmission loss and/or phase delay depends on the direction of propagation [1,2]. In general, reciprocal symmetry can be violated in coupled-mode systems. While the most well-known examples are ferrite-based [1,3], other recent examples can be found in magneto-optic materials [4], Josephson junction Jaynes-Cummings lattices [5,6], parametric [7][8][9][10][11][12][13] and nonlinear devices [14][15][16], and SQUID amplifiers [17]. Most recently a phononic circulator involving two counter-propagating sound waves inside a cavity coupled to three external ports has been described [18].
In the case of two-mode systems, nonreciprocity can be obtained only by cascading multiple elementary transformations in a sequence [9]. However, when the number of modes is greater than two, this is unnecessary and nonreciprocity can be obtained in a single lumped-element device without sequencing, by synthesizing complex coupling coefficients between modes [5,17]. This can be accomplished by use of time-dependent (parametric) couplings [17] or by breaking reciprocal symmetry with magnetic fields [5]. In Reference [17] in particular, Kamal et al. showed that gain isolation in Superconducting Quantum Interference Device (SQUID) amplifiers is a consequence of the coupling between the signal mode and several mixing products driven by multiple parametric pumps derived from the Josephson junction nonlinearities embedded in the SQUID itself. In [5] nonreciprocity was observed in a completely different three-mode system consisting of three transmission lines coupled through a loop of Josephson junctions, where the condition for non-reciprocity was again determined by the phase relationship between the coupling coefficients. This system is reciprocal if and only if the sum of the phases of the coupling coefficients along the junction loop is equal to 0 or π. We might infer from these previous works that the presence of more than two coupled modes, as well as a specific phase relationship between the mode-coupling coefficients, are general properties of various nonreciprocal systems. This suggests that it should be possible to define a set of general requirements for breaking reciprocal symmetry in coupled-mode systems, independent of specific physical implementation.
In this work we associate an abstract graph to a coupled-mode system and use it to derive general conditions for nonreciprocity. We show that the scattering parameters between any two modes can be written as a sum of different terms, each corresponding to a different directed path connecting the modes within the graph. We find that nonreciprocity is a consequence of asymmetric interference between oppositely directed sets of such paths. For instance, reciprocity is guaranteed if the sum of the phases of the mode coupling coefficients along any closed loop in the graph is 0 or π, in analogy with similar conditions found in synthetic gauge field lattices [6]. Moreover amplitude nonreciprocity can be obtained only in the presence of finite dissipation into ancillary modes. Meanwhile, maximum isolation is possible only if specific constraints between the dissipation rates, mode frequency detunings and coupling coefficients are satisfied. In order to connect these considerations to physical device concepts, we introduce a new type of triple-pumped parametric amplifier in which a three-mode system with a closed-loop graph topology is used to obtain both forward gain and reverse isolation without the need for an external circulator, in contrast to traditional parametric amplifiers. Lastly, we apply the graph technique as a tool to analyze gain isolation in DC-SQUID devices.

RECIPROCITY IN COUPLED-MODE SYSTEMS
Resonant and parametrically coupled modes are typically described by the Langevin equations of motion [19]. In this work, we utilize a normalized matrix form of these equations to build a graph representation of multimode coupling dynamics. We begin by considering a set of N r coupled resonators, described by a set of N m internal mode amplitudes. In general, we allow the number of modes to exceed the number of resonators N m ≥ N r such that several modes may reside within a single resonator. This system can be a set of parametrically or resonantly coupled oscillators, but we are not going to assume a particular physical implementation for now, in order to derive results that are as general as possible. We will call ω j and γ j the natural oscillation frequency and total dissipation rate for the resonator in which mode j resides. The complex coupling rate between modes j and k will be denoted by g jk . In order to describe both mode frequency conversion and amplification, we divide the modes into two subsets of p internal mode amplitudes b 1...p and q = N m − p internal mode amplitudes b † p+1...Nm in the frequency basis, with each set of modes corresponding to the driven response at frequencies ω s 1...p and −ω s p+1:Nm . We assume that g jk = g * kj in the case of frequency conversion between modes j and k (for 1 ≤ j, k ≤ p or p + 1 ≤ j, k ≤ N m ) and g jk = −g * kj in the case of parametric amplification (for 1 ≤ j ≤ p and p + 1 ≤ k ≤ N m , or 1 ≤ k ≤ p and p + 1 ≤ j ≤ N m ). For example, the coupling matrix for a 2-mode frequency converter or resonantly coupled oscillator system corresponds to the case p = 2, q = 0, while a conventional parametric amplifier couples two modes, one at a positive frequency (p = 1) and the other at a negative frequency (q = 1) [20]. With this general prescription for defining our mode basis, we can perform an input/output analysis of any system of coupled resonators and calculate the corresponding scattering matrix connecting the vector of input fields b in to the vector of output fields b out (see Appendix A). The scattering matrix can in general be expressed as where M is a N m × N m matrix given by and where we have defined an overall normalization prefactor and environmental coupling matrix, In general γ j ≥ γ ext j , with the equal sign if the internal dissipation rate of mode j is zero. The diagonal elements of M are complex normalized detunings between the driven response frequencies and the natural resonator frequencies and include the mode dissipation rates, while the off-diagonal elements are normalized coupling coefficients. In principal, the normalizations above are unnecessary, but they allow one to cast the Langevin equations of motion (see Appendix A) in a simple matrix form (as M , the "Langevin matrix"), emphasizing the underlying structure of the connections between modes. We show below that this structure, as revealed in the graph representation of M , allows one to draw immediate conclusions about the (non)reciprocity in multimode coupled systems.
We start by introducing a formal definition of reciprocity. The scattering matrix (Equation 1) describes a fully reciprocal system (transmission between any two modes is reciprocal) if it obeys the constraint [21,22], Here the matrix U φ corresponds to a set of phase shifts in the mode basis where we can assume without loss of generality that N ℓ=1 φ ℓ = 0. The matrix U φ corresponds therefore to a redefinition of the mode phases. A system that violates condition (7) is nonreciprocal and is characterized by an asymmetric phase and/or amplitude transmission coefficient between at least one pair of modes. Since the environmental coupling matrix K in Equation 6 is diagonal, it follows from Equation 1 that reciprocity is guaranteed if and only if the Langevin matrix obeys the same similarity transformation as in Equation 7, for some set of phase shifts {φ ℓ }. After substituting β jk = ±β * kj in Equation 9 we obtain where ∠β jk is the phase of the coupling coefficient β jk between modes j and k and n jk is an integer. The system of Equations 10 has at most a set of N m (N m − 1)/2 independent relations (if all off-diagonal elements of M are nonzero) in N m − 1 unknown phases {φ ℓ }. The system is therefore overdetermined and the reciprocity condition (Eqn. 7) is, for parametrically coupled systems, satisfied only for specific choices of coupling phases. In fact, this underscores the peculiar fact that in any parametrically coupled multi-mode system (N m ≥ 3) it is actually difficult to satisfy the reciprocity condition (10) without tuning the pump phases specifically.

GRAPH REPRESENTATION OF THE LANGEVIN MATRIX
The conditions for Equations 10 to have a unique solution, and therefore for the system to be reciprocal, can be visualized by drawing a graphical representation of the matrix M (see Figure 1). The graph has N m vertices. We associate a weighted directed branch ("edge") between vertices j and k to every element β jk in the matrix M . The diagonal elements ∆ j and −∆ * k of M are the weights of self-loops in the graph. The conditions for the entire system to have a unique solution can be obtained by summing Equations 10 along every closed loop L, yielding a simple result, where k L is an integer. Equation 11, reminiscent of Kirchoff's first law in electrical circuits, is analogous to the reciprocity condition in Jaynes-Cummings lattices described in [6]. We stress, however that Equation 11 is applicable to any linear coupled-mode system bearing arbitrarily complex connections and is not confined to the nearest-neighbor, planar graph couplings shown here. Moreover, in the case of parametrically coupled devices in particular, the vertices of our graphs represent modes that may physically correspond to energy that is differentiated by both space and frequency. As such, the graphs are a very general representation of a coupled mode system and can describe coupling of energy at several different frequencies as well as physical ports.
In order for a system to be nonreciprocal, the condition 11 must be violated. This, however, only expresses a minimal requirement for phase nonreciprocity. We can further define a condition for which a system will also demonstrate amplitude nonreciprocity. Specifically, we may ask which conditions allow the amplitude isolation ratio |S jk /S kj | to be different than 1. From (1) it follows that the scattering coefficient S jk scales with the corresponding element of the inverse Langevin matrix, M −1 jk . The problem of calculating the scattering coefficients is therefore reduced to the problem of inverting the Langevin matrix which, as we show below, can be conveniently cast in terms of directed graphs.
The inverse Langevin matrix can be expressed in terms of the transpose of its cofactor matrix [23] The scattering coefficient is then given by where δ jk is the Kronecker delta function. The elements of the cofactor matrix C kj are equal to the determinant of the matrix obtained from M after removing the k-th row and the j-th column (i.e., the kj minor), multiplied by (−1) j+k [23].
The determinant of M and the cofactor C kj can be easily calculated from the graph in Figure 1(a) by the following procedure (see Refs [24,25]). Given a graph, we call a subgraph G a "permutation subdigraph" if every vertex in G has one and only one outgoing edge and one and only one incoming edge [25]. The total weight w(G) of a subgraph G is the product of the weights of all of its edges times (−1) c+Nm , where c is the number of cycles in G (by cycle we mean a sequence of distinct connected vertices, where the first and last vertex coincide). The determinant of the matrix M is the sum of weights of every possible permutation subdigraph of the graph associated to M [24,25], as shown in Figure 1 The cofactor matrix C kj is the coefficient of the M kj element of the matrix M in the Laplace expansion of the determinant |M |. Therefore, in order to calculate C kj , we need to compute only the permutation subdigraphs in the expansion of |M | that contain the kj edge, with the kj edge removed. In other words, for every path p FROM vertex k TO vertex j, we consider the collection of all the permutation subdigraphs G p r of the graph of M containing p, with the kj edge removed as in Figure 1(c). The cofactor C kj is then given by where p varies over all the paths connecting vertex j to vertex k and r varies over all the subdigraphs containing the path p. The minus sign is due to the fact that, by removing the kj edge, we reduced the number of cycles in the subdigraph by 1. An example of this procedure is outlined in Figure 1(b). The graphs {G p r } in the expression for the cofactor elements have a simple physical interpretation: they represent the possible scattering mechanisms that connect modes j and k.
Since, from (13), the ratio between the scattering coefficients S jk /S kj is equal to the ratio between the corresponding cofactors C kj /C jk , the problem of analyzing isolation in a given multi-mode system can be reduced to evaluating the differences between the specific graphs that connect j to k and vice versa.
As a more concrete example, we apply this procedure to a simple three-mode system in which all modes are coupled via frequency conversion (β jk = β * kj ) and are therefore described by the following matrix (see Figure where ∆ j and β jk are normalized detunings and couplings as before, and we have dropped the overall normalization prefactor for convenience. The isolation between ports 1 and 2 can be calculated by use of the graphs representing the cofactors (see Figure 2(b) and (c)) yielding the expression, where, in the last line, we have removed all normalizations.
From (18) we observe that if the total dissipation of mode "3," γ 3 = 0, then S 12 = S * 21 and there is no isolation between modes 1 and 2. In this case mode 3 serves as an internal lossless mode and the system behaves as a simple reciprocal two-mode system. Here, dissipation plays a crucial role in violating reciprocity, as it breaks the symmetry between scattering elements which are otherwise related by complex conjugation. Assuming finite dissipation in mode 3, we find the condition for maximum isolation I 12 = 0, In other words, to achieve maximum isolation, the sum of the phases of the couplings along the loop connecting the three modes, φ 13 + φ 21 + φ 32 , must be equal to the detuning angle, tan −1 [γ 3 /(2(ω s 3 − ω 3 ))]. At zero detuning in mode 3, ω s 3 = ω 3 , this loop phase must be π/2 for perfect isolation. Equivalently, one can interpret this isolation condition, Eqn. 19, as interference between the two possible coupling paths connecting modes 1 to 2 in the graph (Figure 2).

THREE-MODE DIRECTIONAL AMPLIFIERS
Traditional single-stage lumped-element parametric amplifiers, such as the one shown in Figure 3(a), operate in reflection-mode such that a circulator is needed in order to achieve unidirectionality. In [26] a single-stage amplifier with two parametric pumps at the same frequency was introduced, but in addition to forward gain it also had unity reverse gain. The graph methodology Three-mode Delta amplifier discussed in this work, which consists of one flux-biased DC-SQUID connected to 3 resonators at frequencies ω1,2,3 and its graph representation. The SQUID is equivalent to a nonlinear inductor modulated by the magnetic flux through the SQUID loop. One pump at frequency ωp1 = ω2 + ω1 is used to obtain parametric amplification between mode 1 at frequency ω1 and mode 2 at frequency ω2. The second pump at frequency ωp2 = ω3 − ω1 generates frequency conversion between mode 1 and mode 3 at frequency ω3. A third pump at frequency ωp3 = ωp1 + ωp2 closes the loop and generates parametric amplification between modes 2 and 3. In this configuration non-reciprocal trans-gain can be obtained between modes 1 and 2. (c) and (d) subgraphs used to calculate the forward and reverse gain between modes 1 and 2.
introduced here allows us to revisit the problem of gain and isolation from another perspective. In particular, our discussion above identified the importance of introducing a multi-edge loop to obtain nonreciprocity. Here we analyse a closed-loop amplifier which we term the "Delta" amplifier, shown schematically in Figure 3(b) both as a device concept and as its corresponding graph.
As drawn, the Delta amplifier requires three independent pumps to generate three pairs of coupling edges {β jk }. Two of these pairs create amplification (i.e., β 23 = −β * 32 , β 12 = −β * 21 ), while the third pair couples through frequency conversion (β 31 = β * 13 ). One concept for how this might be implemented is shown schematically in Figure 3(b) as three resonators that share a common current anti-node in which a flux-driven SQUID is operated as a parametrically modulated inductor. In principle, one might drive two of these couplings using the same pump frequency [27] [28], but for clarity we discuss only the case in which each of the three couplings is driven independently, which has the advantage of increased flexibility due to the increased number of degrees of freedom. Regardless of particular implementation, we can recognize the presence of a closed loop connecting all three modes, as well as the integration of amplification as a coupling process. From this bare description we recognize that with three independent pumps we are constrained to optimizing only three of the nine possible scattering elements. In laboratory use, we can choose to optimize for amplification in the "forward" direction (|S 21 | > 1), isolation in the "reverse" direction (S 12 ∼ 0), and low input reflection coefficient (S 11 ∼ 0) to prevent unwanted reflections of signals propagating away from a device-under-test.
We begin by writing down the expressions for the forward and reverse transmission from Equation 13, where |M | is the determinant of the coupling matrix represented in the graph shown in Figure 3(b). Assuming ideal isolation, S 12 = 0, we find an expression similar to Equation 19, which provides an overall constraint on the coupling amplitudes and phases as a function of the complex detuning of mode 3. Under this primary constraint, we then set the input reflection coefficient to zero, giving a constraint on the determinant, Taken together, these constraints yield an expression for the forward transmission gain, If the device dissipation is symmetric and external coupling dominates, {γ k } = {γ ext k }; in the limit of small detunings (ω s k − ω k ) ≪ γ k this reduces to Note that this form gives the gain at a fixed frequency where the input match is perfect (Equation 23). High gain can be obtained when |g 23 | ≈ √ γ 2 γ 3 . The remaining two couplings can then be chosen in order to satisfy Equations (22) and (25). We have calculated the resulting S-parameters for the Delta amplifier in Figure 4 under various constraints on the pump detunings for the case of uniform detunings γ k = γ. In practice, one might fix the pump frequencies, locking the three detunings together, and measure the S-parameters shown in Figure 4(d). For this example calculation, we selected values for the pump amplitudes to fix the nominal power gain to ∼20 dB while obtaining significant isolation. Finally, we find that the calculated output S 22 shows reflection gain at the output port, and is consistent with the limitations of having only three parameters (the independent pumps) to optimize three of the four scattering parameters between ports 1 and 2. However, this might not be a limitation in practical use where one is more concerned with the backaction of the amplifier on a quantum circuit on the input side (port 1).

DIRECTIONALITY OF DC-SQUID AMPLIFIERS
DC-SQUID amplifiers are a type of superconducting amplifier with noise temperatures approaching the quantum limit [29][30][31]. Since DC-SQUIDs also have inherently low power dissipation in addition to high gain and directionality [32], they have been the subject of several applied superconductivity efforts to produce reliable The SQUID can be modelled as a parametric coupled-mode system. The modes at frequencies ω ± nωJ are coupled through internally generated pump signals at multiples of the Josephson frequency ωJ . At each frequency there are two doubly degenerate modes, corresponding to the common and differential excitation of the junction phases.
amplifiers for microwave quantum information measurements. In [17], Kamal et al. showed that the directional-ity of the SQUID is a consequence of multiple parametric amplification and frequency conversion processes due to frequency mixing between the input microwave signal and internally-generated parametric pumps at multiples of the Josephson frequency. This approach casts the DC-SQUID as a multi-mode coupling problem and, as such, is also amenable to a graph description. Here we describe the parametric coupling in the SQUID following the description in Reference [17] by use of graphs to map the connections between modes. In particular, we will direct our attention to the interference and dissipation that must be present to generate nonreciprocity in the presence of gain. A DC-SQUID consists of a pair of Josephson junctions inside a superconducting loop as in Figure 5. Each junction is characterized by a phase difference φ 1,2 across its terminals, but it is more convenient to recast the dynamics of the DC-SQUID in terms of the common and differential phase difference φ C,D = φ 1 ± φ 2 [17]. φ C,D can be expressed as the sum of fast oscillating terms at multiples of the Josephson frequency ω J = 2πV dc /Φ 0 , where Φ 0 is the magnetic flux quantum, and slow oscillating terms at the signal frequency ω. The fast oscillating terms act as parametric pumps causing mode conversion between ω and the mixing product frequencies ω ± kω J . The lowest-dimensional model that can capture nonreciprocity involves two parametric pumps and a total of 10 modes at 5 frequencies ω ± 2ω J , ω ± ω J , ω. Each mode is doubly-degenerate, because for each frequency there are a common and a differential mode excitation. Some scattering subgraphs to describe the SQUID are shown in Figure 6, where, for example, β cc ω,ω−ωJ indicates the coupling rate between the common modes at frequency ω and ω − ω J . The input and output modes are shown at the very bottom and top respectively. Calculation of the mode coupling rates as a function of the internally generated pump amplitudes is detailed in Appendix B.
The number of subgraphs to describe the SQUID is too large for analytical calculations. However if the SQUID bias current is high enough, we can look for an approximate expression for the SQUID isolation for small values of the bias parameter ǫ = I c /I b , where I c is the critical current of the junctions and I b is the bias current. The pump amplitudes and the coupling coefficients can then be calculated for small ǫ and we can simplify our calculations by considering only the subgraphs of leading order in ǫ. These subgraphs are shown in Figure 6. The order of the graph is given by the product of its edges and therefore depends on the length of the path connecting the input to the output modes as well as the order of the edge weights (coupling rates). The order of the coupling rates increases with the difference between the mode drive frequencies (for example β cc ω,ω−kωJ ∼ ǫ k ). Moreover a small (∼ ǫ 3 ) direct coupling exists between modes at the same frequency, as explained in Appendix B.
The dissipation rates for the common and differential modes when the SQUID is terminated into ideal infinite transmission lines, are given by with β L = (2LI c )/Φ 0 , β c = 2πI c R 2 C/Φ 0 , ω c = 2πI c R/Φ 0 , L is the inductance of the SQUID loop, and R and C are the junction shunt resistance and capacitance. Since the SQUID is a nonresonant system, the mode frequencies ω i in Equation (3) are zero. We can approximately compute the SQUID isolation by use of the lowest-order subgraphs shown in Figure 6. The first conclusion we can draw, based on our previous discussion, is that non-zero dissipation rates in (28) are essential to obtain nonreciprocity in the SQUID. A change in the impedance terminations at the various mode frequencies ω ± kω J will change the amount of nonreciprocity, viz., modes that are terminated in a short or an open circuit (corresponding to zero dissipation rates) will not contribute to the SQUID isolation. Moreover, it is important that γ C = γ D . In fact, by taking into account the symmetries between the mode couplings, the weight of the 2 nd order subgraphs in Figure 6 is given by apart for a common factor equal to the product of the internal modes normalized detunings. If γ C = γ D = γ, then 1/∆ * c− + 1/∆ d+ = 2ω/(ω 2 − ω 2 J + γ 2 ) and w ǫ 2 becomes an even function of iγ C,D . The same is true for the 3 rd -order subgraphs in Figure 6. For the 4 th -order subgraphs, every term proportional to products of complex detunings 1/∆ * C ∆ D , 1/|∆ C,D | 2 is also an even function of iγ C,D if the dissipation rates are equal. As we discussed previously, if the scattering coefficients are an even function of the dissipation rates, there is no explicit dependence on the imaginary terms and therefore no isolation, because all the coupling rates are complex-conjugated upon inversion of the direction of propagation. Isolation can be restored by increasing the resonant coupling rate β cd ω,ω . In a real device we expect such a rate to be higher, due to the stray capacitive coupling between the input and output as well as other parasitics [32]. Finally, in Figure 7 we show the computed isolation and power gain of the SQUID as a function of the bias points for the exact numerical solution and for the lowest-order approximations. We find that graphs up to 4 th order in ǫ are enough to accurately describe the isolation properties of the SQUID. The exact power gain was also computed numerically by extracting the impedance matrix from the scattering matrix of the SQUID and then calculating the gain as discussed in [17]. In principle the power gain may also be computed from the subgraphs, but the number of necessary graphs is large, even at lowest order. In this case, the graphs help identify the elements of a complex multimode scattering problem that contribute to the isolation properties.

CONCLUSION
In this work we introduced a graph theoretical representation that can be used as a combinatorial accounting tool to analyse nonreciprocity in coupled-mode systems. An abstract graph associated to the coupled-mode system was used to compute the input/output scattering coefficients between modes and derive general constraints that need to be satisfied in order for nonreciprocity to occur. In order for a multi-mode system to be reciprocal the sum of the phases along any loop has to be 0 or π. If this condition is violated, phase and/or amplitude nonreciprocity is present. Since the elements of the scattering matrix can be computed from individual sub-graphs containing paths that connect modes, this condition for reciprocity can be interpreted as interference between different permutations of connecting paths. We also find that dissipation in the remaining disjoint vertices/modes is crucial for amplitude nonreciprocity/isolation. Specifically, it breaks the symmetry between forward and back-ward multi-mode scattering processes between two modes such that one may obtain constructive loop interference in one direction, and destructive interference in the reverse direction.
Although the results obtained using this graph methodology can be obtained using more conventional approaches, graphs offer a new perspective that can be useful for immediately identifying the roles of various elements in a complex multi-mode scattering problem. This can be leveraged as a design tool, allowing one to reduce a set of desired device characteristics to a mini- mal description that can be utilized to design new amplifier/converter concepts such as the Delta amplifier. We also showed, with the complex example of the DC-SQUID, that graphs aid in identifying the critical elements that create isolation in the presence of forward gain. Finally, graphs may provide a useful approach to the question of nonreciprocity in recent multi-mode parametrically coupled hybrid systems such as the electricalmechanical-optical (three-mode) transduction bridge [33] and electro-mechanical amplifiers [34]. In general, problems like these are particularly well-suited to a graphbased approach since many of the interesting properties one is interested in (gain, isolation, gain-bandwidth product, etc.) actually arise from the structure of the coupling network (the topology of the graph) and not the physical particulars of a given implementation.
(37) In Equations 36 and 37 ω j and γ j are the natural frequency and total dissipation rate for the resonator in which mode j resides, while the "*" indicates complex conjugation, and g jk is the complex coupling coefficient between modes j and k. We further assume that g jk = g * kj if 1 ≤ j, k ≤ p or p + 1 ≤ j, k ≤ N m and g jk = −g * kj for 1 ≤ j ≤ p and p + 1 ≤ k ≤ N m , or for 1 ≤ k ≤ p and p + 1 ≤ j ≤ N m . For example, the coupling matrix for a 2-mode frequency converter or resonantly coupled oscillator system corresponds to the case p = 2, q = 0, while a conventional parametric amplifier couples two modes, one at a positive frequency (p = 1) and the other at a negative frequency (q = 1) [20]. With this general prescription for defining our mode basis, we can perform an input/output analysis of the system (Equation 32) and calculate the vector of output (scattered) fields b out by means of the relation [35]: In general, the phase factor ξ j depends on the nature of the coupling to mode j. For simplicity, we set e iξj = −1, as one would obtain with a small series capacitance or weakly reflecting mirror. Other choices of coupling will have the effect of modifying the identity matrix on the r.h.s. of Equation 1. By substituting Equation 33 and 34 into Equation 32 and by use of the input/output boundary conditions (Equation 38), we can obtain the expression for the scattering matrix in Equation 1.

Appendix B: Coupled-modes analysis of the DC-SQUID
A DC-SQUID consists of a pair of Josephson junctions inside a superconducting loop as in Figure 5. Each junction is characterized by a phase difference φ 1,2 across its terminals. Defining the common and differential phase difference φ C,D = φ 1 ± φ 2 , the dynamical equations describing the system are [17]: where ω B = 2πI b R/Φ 0 , and ω c = 2πI c R/Φ 0 . I c is the junction critical current and I b is the bias current of the SQUID, Φ 0 is the magnetic flux quantum, L is the loop geometric inductance and R is the equivalent shunt resistance across the SQUID. The system (39) can be linearized around the following bias point: where δ c , δ d are small perturbations. In [17] it was shown that by applying a harmonic balance procedure, the phases δ c,d can be expressed as (43) where the components p c,d k are internally generated pumps at frequency kω J = 2πkV dc /Φ 0 , where V dc is the static DC voltage across the SQUID. The pumps cause parametric frequency mixing between the various signal modes s c,d k . In [17] the harmonic balance equations were solved for the case M = 3 to show how nonreciprocity arises from the interference between the different parametric conversion paths. As shown in the main text, we can employ the graph representation of a coupledmode system to write explicitly the interference conditions needed to maximize the amount of amplitude nonreciprocity between two modes. In the M = 3 case we can describe the SQUID as a 10-mode system. The 10 modes have frequencies ω ± 2ω J , ω ± ω J , ω and are doubly degenerate (for each frequency there are a common and a differential excitation). In order to calculate the scattering parameters, we need to find the coupling coefficients between the modes as a function of the internal pump amplitudes. This can be obtained by linearising Equations (39) for small signal mode-amplitudes and equating terms at the same frequency in order to obtain a set of coupled mode equations 32. After taking symmetries into account, the system can be described by the following five coupling coefficients g cc ω,ω−ωJ = iǫ(cos φ e − ipi c 2 cos φ e − p d 2 sin φ e ) (44) g cd ω,ω−ωJ = iǫ(i sin φ e − p c 2 sin φ e − ip d 2 cos φ) (45) g cc ω,ω−2ωJ = iǫ(ip c 1 cos φ e − p d 1 sin φ e ) (46) g cd ω,ω−2ωJ = iǫ(ip d 1 cos φ e − p c 1 sin φ e ) (47) g cd ω,ω = iǫ(Im(p d 1 ) cos φ e − Re(p c 1 ) sin φ e ), (48) where for example g cc ω,ω−ωJ corresponds to the coupling between the common modes at frequency ω and ω − ω J .