Circulation by microwave-induced vortex transport for signal isolation

Magnetic fields break time-reversal symmetry, which is leveraged in many settings to enable the nonreciprocal behavior of light. This is the core physics of circulators and other elements used in a variety of microwave and optical settings. Commercial circulators in the microwave domain typically use ferromagnetic materials and wave interference, requiring large devices and large fields. However, quantum information devices for sensing and computation require small sizes, lower fields, and better on-chip integration. Equivalences to ferromagnetic order -- such as the XY model -- can be realized at much lower magnetic fields by using arrays of superconducting islands connected by Josephson junctions. Here we show that the quantum-coherent motion of a single vortex in such an array suffices to induce nonreciprocal behavior, enabling a small-scale, moderate-bandwidth, and low insertion loss circulator at very low magnetic fields and at microwave frequencies relevant for experiments with qubits.


I. INTRODUCTION
In order to scale up quantum computing devices, the supporting hardware for qubits must appropriately scale down. Essential components include elements for signal routing and isolation, such as circulators, as well as lownoise amplifiers, filters, and cooling systems. While development in improving miniaturization and effectiveness of these components progresses apace [1][2][3][4][5][6][7][8][9], traditional circulators remain bulky and require high magnetic fields for operation.
Creating a circulator requires nonreciprocal behavior induced by implicit or explicit breaking of time-reversal symmetry. This is then combined with microwave engineering to enable unidirectional transmission with little to no loss in a clockwise or counterclockwise fashion over a large operating bandwidth. The nonreciprocal behavior observed in commercial ferrite junction circulators, which are commonly utilized in smaller systems, is a result of the Faraday effect. The interaction between a magnetic field, e.g., from a permanent magnet, and the central ferrite discs of the device create circulation via wave interference and magnetically induced gyration [10]. As a result, commercial ferrite junction circulators are fundamentally limited in size by the signal wavelength. However, their passive operation necessarily eliminates any additional control hardware. Nonetheless, with their large size and reliance on strong fields, the prospect of employing ferrite circulators in large quantum computing schemes is a challenge for the field, though some progress in wavelength-size devices [11] has been made. This requires finding small-volume, on-chip solutions that operate at low magnetic field and enable the isolation of quantum devices and corresponding precise qubit control and measurement. * brr215@umd.edu One partial solution to this problem is generating nonreciprocity via the quantum Hall effect, where large magnetic fields break time-reversal symmetry and quantum Hall edge modes are utilized in passive devices to observe circulation [12,13]. Other proposed solutions are actively controlled and take a nonmagnetic or small-field approach, using frequency conversion in driven systems and irreversible dynamics (loss) or synthetic gauge fields in a Floquet basis to generate nonreciprocity [3,4,[14][15][16][17][18]. However, due to their active nature, they necessarily require additional control hardware that is undesirable in larger-scale systems. To date, the only passive, low-field, subwavelength approach uses particle tunneling in small Josephson junction arrays (JJAs) or quantum phase slip junctions, utilizing their behavior in magnetic fields-the Peierls phase in tunneling-to break time-reversal symmetry [19,20]. Both of these proposed devices suffer from challenges in implementation due to susceptibility to charge noise, difficulties in implementing robust quantum phase slip devices, and poor performance in operational bandwidth.
Here we address these challenges by moving from a Cooper pair tunneling device to a persistent current (vortex) tunneling device via insertion of additional Josephson junctions. Our approach, shown schematically in Figs. 1a and 1b, builds upon the three-island Josephson junction loop proposed by Koch et al. [19] and further explored by Müller et al. [20]. In contrast to those works, we suggest operating in a non-charge-conserved, intermediate regime of Josephson energy and charging energy being similar. By keeping the Josephson energy similar to the charging energy, we can maintain quantum-coherent vortex dynamics rather than having the microwave response behavior be dominated by spin waves, as occurs at large Josephson energy. We show a reduction in sensitivity to charge noise while enabling large bandwidth operation.
To understand this performance, we draw an anal-arXiv:2010.04118v2 [quant-ph] 14 Jun 2021 ogy to bulk superconductors, where the Meissner effect leads to persistent currents around holes in a superconductor. In the limit of large tunneling rates, these currents are fixed to each hole. This is related to the dual picture of a lattice of vortices with large effective mass trapped around each hole. However, as the kinetic term increases by decreasing capacitance, and phase-slips are made available by the introduction of tunnel junctions (Josephson junctions), the vortices are able to move. Some devices, such as the flux qubit [21,22], have demonstrated the ability to prepare a superposition of different vortex states for a small circuit [23,24]. We suggest that the successful operation of our proposed device arises by exploiting the behavior of such coherent vortices in a magnetic field, analogous to the behavior of the ferrite system where the XY model of the junctions acts as an effective ferromagnet. Here, we propose circulation is achieved via the rotation of these persistent current vortices: the voltage induced by an incoming signal generates a vortex in the circuit, this vortex then rotates about the circuit, leading re-radiation into each port.
This paper is organized as follows. Section II introduces the central circuit model to be explored, motivating this choice by presenting an investigation of vortices in the model and their potential role in device functionality. Section III couples the central circuit to external circuits to investigate transmission properties of the system. A theoretical analysis in the single excitation limit is presented, resulting in the single photon scattering matrix (S-matrix), followed by conditions for ideal circulation. Numerical results related to bandwidth, circulation, and noise performance are shown in Section IV, and we conclude and provide an overview of results in Section V. In the appendices that follow, we explore an alternative circuit model to the one presented in Section II, followed by a numerical examination of two flux points in addition to the one investigated in Section IV. We then discuss various parameter effects and provide additional details to supplement the noise analysis presented in Section IV.

A. Model circuit
Here we investigate a minimal model of a circuit for our proposed system which consists of three superconducting islands connected to each other and a central grounded island via Josephson junctions, as shown in Fig. 1b, all threaded by an external flux Φ. We note that the central ground can, in principle, be replaced by a capacitive island with C 0 C J , the nominal capacitance of each island. This alternative circuit is shown in Fig. 1c and explored in detail in Appendix A. For the circuit in Fig. 1b, the overall Hamiltonian is: The operatorsN andφ are the vectors of number and phase operators for each island of the circuit, with individual elements satisfying the commutation relation [φ j ,N k ] = iδ jk . E J is the Josephson energy of each junction and E Σ is the circuit's charging energy, defined as E Σ = e 2 2CM where C M is the largest eigenvalue of the circuit's capacitance matrix C, defined in Eq. (2). K −1 is the dimensionless inverse capacitance matrix appropriate for the circuit, found by factoring C M from C −1 . As it is written, H CC only stipulates the central circuit be a superconducting, three-island system containing identical Josephson junctions. By specifying the elements of K −1 and the tunneling potential V (φ, Φ) we establish the particular central circuit for the system explored in this work. This allows us to compare to prior work without substantial computational difficulty.
We define K −1 by inverting the capacitance matrix C, where C S = C C + C G + 3C J with C J the capacitance of each Josephson junction, C C the coupling capacitance, and C G the residual capacitance to ground. For this matrix, C M = C C + C G + 4C J is the largest eigenvalue. Factoring C M from C −1 we can then find: where δ = CJ CC+CG+CJ . The inductive terms, corresponding to tunneling through the Josephson junctions, are described by where we take i+1 = 4 and i = 1 to be the same. We note that this potential is equivalent to the XY model where the phase ϕ represents the angle of the U (1) degrees of freedom.
The acquired Peierls phase factor, 2πA, is equal to the line integral of the vector potential in a counterclockwise fashion around each third of the device: This A corresponds to the frustration of the circuit, i.e., the flux penetration per loop of the device in terms of the magnetic flux quantum.
In what follows, we denote the eigenstates of H CC as |ε m with eigenenergy ε m , and label the ground state |G . These are parametric functions of A, and periodic in A.

B. Entering the vortex regime
The behavior and properties of the central circuit, shown in Fig. 1b, depend on the relative size of the two characteristic energy scales in the system: the energy associated with Cooper pair tunneling between islands, E J , and the energy associated with putting a single electron on an island, E C ∼ E Σ . In the two extreme energy limits, E J E C or E J E C , the relevant excitations are charges or vortices, respectively [25][26][27]. Note that charges have the same properties in the chargeddominated regime that vortices have in the tunnelingdominated regime [26,28]. In the charge-dominated regime, where E J E C , charges are well-defined, massive particles that feel a Lorentz force in the presence of an electromagnetic field and acquire an Aharanov-Bohm phase (represented in the charge basis as a Peierls phase) when moving around an area containing a magnetic field [29,30]. However, in this regime, external voltages and other sources of charge disorder, such as offset charge, have a dramatic effect on circuit behavior due to the very small capacitances necessary to achieve large E C .
In the classical tunneling limit, where E J E C , vortices are the topological excitations of the XY model that behave as massive (and essentially classical) particles. They respond to the presence of external current by experiencing a force perpendicular to the current flow (a Lorentz force), and the presence of external charge on the islands affects a moving vortex the same way a magnetic field affects a moving charge (an Aharanov-Casher effect for vortices) [26][27][28]31]. In this way, charges and vortices are dual to each other, with these mirrored properties being manifestations of the charge-vortex duality that exists in JJAs. However, the vortex domain has two characteristic excitations: the topological excitations, as just described, as well as plasma oscillations (spin waves in the XY model language) that correspond to the response of a classical circuit of inductors and capacitors with values set by the expectation value of the Josephson inductance for the nominal ground state (see, e.g., Ref. [32]). The plasma oscillations do not break timereversal symmetry; in essence, the Meissner effect traps flux in the vortices, and thus, all time-reversal symmetrybreaking terms arise only from vortex motion. When E J E C , this motion becomes exponentially slow due to the challenge of tunneling between allowed current configurations, and thus, one cannot make a circulator in this regime.
We wish to identify an energy regime using this duality, where time-reversal symmetry is broken in the microwave regime and charge noise is not relevant. We thus focus on E J ∼ E C as a regime where vortices can move and charge noise may be reduced, enabling the potential for noise-resistant nonreciprocity in the form of unidirectional signal transmission that is robust in the presence of random variation in offset charge. In this way, the model may act as an ideal circulator.
This intermediate regime has been partially explored by Koch et al. [19] and Müller et al. [20]. However, in those works, the system was defined with a chargeconserving three-island Josephson junction loop, and total charge on the small capacitance region was fixed. In what follows, we leverage the connection of each island to ground via an additional Josephson junction to enable large charge number fluctuations and suppress charge noise. This leads to circulator performance driven by vortex behavior rather than Cooper pair behavior.

C. Quantum behavior of vortices in the intermediate tunneling regime
To explore the vortex picture, we set E J = 4E C to ensure that there are sufficient phase fluctuations to warrant a quantum treatment of the vortex dynamics, and numerically diagonalize H CC [given by Eq. (1)] for a range of A from 0 to 0.5, as the eigenvalues are symmetric in A → −A and periodic (same for A + 1 as for A). We denote the corresponding energy eigenstates |ε n and eigenvalues ε n , with implicit dependence on A understood. We note for our chosen parameters (E J = 4E C ≈ 45 GHz) and eventual optimum operating point (near A = 0.2478), that the lowest energy spin-wave excitation occurs at 20 GHz, far higher than the low-energy excitations (in the few GHz regime) we consider here. As a result, we expect the predominant dynamics to be that of vortices, not spin waves.
We approximate the full Hilbert space using a charge basis |n 1 , n 2 , n 3 for Cooper pairs on each island, with |n i | ≤ 4. We checked truncation and found no difference by incrementing the maximum allowed charge by one. This creates a Hilbert space dimension of 9 3 and is straightforward to solve on a personal computer. The energy eigenvalues of the circuit as a function of A are shown in Fig. 2a. We note the presence of level-crossings at A ≈ 0.2 and A ≈ 0.26 within the first four energy levels due to the high symmetry of the system.
We now consider the low-energy space of this system, and examine the presence and behavior of localized vortices. For low-frequency excitations (in the few GHz regime), and starting from the circuit ground state, we can expect to only explore superpositions of the first few eigenstates. To see the vortices, we identify the operators associated with the persistent currents around each of the small loops in our central circuit, shown in Fig. 1b: where we take i + 1 = 4 and i = 1 to be the same. To find the true current in the system, one multiplies Eq. (6) by I C , the single-junction critical current. Fig. 1d shows a schematic representation of each of these loop currents for our model of interest. A local vortex is identified by a nonzero loop current in one or two of the three loops. Note that if all three loop currents are equal, there is only a persistent current circulating about the exterior. For the energy eigenstates, Fig. 2b shows that the expectation values of each of the three loop currents are the same, and thus the eigenstates do not correspond to localized vortex states, though each eigenstate has a distinct persistent current around the exterior of the circuit. We also see that the level crossings manifest themselves in the subsequent persistentcurrent analysis, proving them to be regions of flux that create dramatic changes in behavior due to changes in the nature of the ground state, as seen in Figs. 2b, 2c, 3d, and 3e. Our chosen operating point for circulation, shown with a vertical grey line, is away from these crossings.
While there are no localized vortices for energy eigenstates, an incoming microwave photon can lead to the excitation of superpositions of circuit states. To examine this, we consider the matrix elements of the three loop currents over the first three eigenstates {|ε 0 (= |G ), |ε 1 , |ε 2 }, denoted I kl i = ε k |Î i |ε l , whose magnitudes are shown in Figs. 2b and 2c. The symmetry of the circuit requires that the absolute value of each of the loop-current operator matrix elements be independent of i, that is, the loop to which they correspond. However, as we now show, superpositions of eigenstates exhibit localized vortices due to nontrivial phases of the off-diagonal matrix elements I kl i for k = l. We write the argument of the I mn i matrix elements as φ mn i . The differences φ mn i − φ mn i+1 only take the values {0, ±2π/3}, as expected by the three-fold symmetry of the system. We now consider their role in circulation. Examining a superposition c m |ε m + c n |ε n with m > n, we see that the i th loop-current operator has an expectation value at time t of where φ mn i is the argument of the off-diagonal matrix element I mn i , θ mn is the argument of c * m c n , and ω mn = (Em−En) . Thus, we find a superposition of two states has, in general, a local vortex, and that these local vortices circulate (countercirculate) about the system when φ mn i − φ mn i+1 = ±2π/3 at a frequency ω mn . This is shown schematically in Fig. 3a.
Accordingly, we define the directional function for the difference in φ mn i variables modulo 2π as to allow us to plot the behavior of potential local vortices for different values of m, n, and A, including the existence or absence of vortex circulation and what direction. Θ mn takes the value +1 for clockwise circulation, -1 for counterclockwise circulation, and 0 for no circulation (absence of a local vortex). The directionality of vortex circulation associated with the two-state superpositions 0 − 1 and 0 − 2 as a function of A are shown in Figs. 3d and 3e, respectively.

D. Connecting vortex behavior with the microwave response
We now consider what happens if an oscillating voltage is applied to the j th island of the central circuit via a coupling capacitor. We denote the coupling as Ω j (t)B j , withB = K −1N the vector of dimensionless voltages on each island. Using linear response theory, we can evaluate the expected loop currents at a later time for a perturbation Ω j (t) ∝ δ(t) as where we use P as in a dimensional setting this quantity has units of power, scaling as 2eIC CM , and can be interpreted as the power in the circuit induced by the drive. Using the Green's function relations, we then also expect that a monochromatic excitation at frequency ω gives the Fourier transform of this value for the expected long-time response of the different loop-current operators. We will also examine, for the microwave response, the voltagevoltage correlation function, whose Fourier transform will be identified later as the transmission matrix of the circuit.
Working this out explicitly, where γ → 0 + is taken to ensure causal ordering and . We note that a narrowband signal will preferentially excite states near the signal frequency, as is apparent inP ij (ω)'s denominator.
Examining Eq. (11), we can interpret the linear response behavior as the evaluation of Î i for a state of the form |G − i n |ε n ε n | ( j Ω jBj ) |G for Ω j → 0. Thus we can conclude that, insofar asB j produces superpositions of energy eigenstates, we can expect the creation of moving vortices within the array upon an incoming signal. As we expect sinusoidal signals near the 0 − 1 and 0 − 2 transition frequencies, we plot the terms of P i1 (t) corresponding to the first and second excited states in Figs. 3b and 3c [denoting the n th term of P i1 (t) as p n i1 (t)] to illustrate the basic expected dynamics of the system.
Noting that −P 1j (−t) corresponds to the expected voltage at island j due to a current at loop 1, we can see the emergence of circulator dynamics: a voltage creates vortices, which rotate, and re-radiate into the leads. Crucially, we note that the observed rotation of vortices corresponds to highly nonreciprocal behavior, and thus we can expect generically that this central circuit will provide the opportunity for building a circulator. However, substantial practical details must now be worked out. We do so below, focusing on the linear response (low power) regime.

III. CENTRAL CIRCUIT COUPLED TO EXTERNAL CIRCUITS
A. The full circuit model We now consider the central, three-island superconducting circuit capacitively coupled at each island to a resonator and voltage source, as shown in Fig. 1a. These resonators are included to enable effective impedance matching between incoming microwave signals and the central circuit, and will be tuned in frequency and coupling to enable optimal performance. Each resonator port is opened to an external source of microwaves, i.e., a transmission line, allowing signals to flow into and out of the system. The central circuit, specified above in Section II A, is threaded by an external flux Φ via the application of an external magnetic field.
The Hamiltonian for this model is composed of terms corresponding to the transmission lines, resonators, and central circuit, as well as each circuit-resonator interaction and resonator-transmission line interaction. Generally, this may be written as: Taking the weak-coupling limit, the rotating-wave approximation, and the Markov approximation, each term of the Hamiltonian in Eq. (13) can be expressed as: H TL [33] describes each transmission line as containing an identical continuous spectrum of modes over frequency ν. H R [19] is the contribution from the system's three identical resonators, consisting of a single low-lying mode with a characteristic frequency ω R . The operatorsâ † k and a k create or destroy a photon in the k th resonator and satisfy the usual commutation relation: [â j ,â † k ] = δ jk . Similarly,b † k (ν) andb k (ν) create or destroy a photon in the k th transmission line and satisfy the commutation relation: [19,33] are the interactions between each circuit island and resonator and between each resonator and transmission line, respectively. By treating all resonators and all transmission lines as identical, we take the coupling constants, g and κ, to be identical across all three interactions, and will later optimize to ensure the best circulator performance. The operatorB k , as defined in Section II D, is a dimensionless operator corresponding to the dimensionless voltage on the k th island, defined as the k th row ofB = K −1N .
Lastly, H CC is as defined in Section II A and describes the central superconducting circuit. Here we also included applied voltages V gk (shown in Fig. 1a) that can lead to charge offset noise N gk on each island.

B. Single excitation S-matrix analysis
For this analysis, we work in the low-power, single excitation limit, wherein we consider only one photon as input to or output from the system and require that the excitation energy in the system only exist in one component of the system at any one time-in one of the resonators, one of the transmission lines, or as an excited state of the central circuit. Upon enforcing this limit, we establish the S-matrix of the system using input-output theory. First, we write a generic state in this single excitation manifold as: where d k , c n , and f k (ν) are time-dependent complex probability amplitudes for each state, and we abbreviate, with a slight abuse of notation, the ground state |0, 0, 0, G, 0, 0, 0 as |G , and the excited circuit state with no photons in any resonator or transmission line |0, 0, 0, ε n , 0, 0, 0 as |ε n . The summation over n is a summation over only the excited states of the circuitit excludes the ground state as this state is not in the manifold being considered. By using the Schrödinger equation, we establish equations of motion for the time-varying complex probability amplitudes appearing in Eq. (15). With these equations of motion, we can utilize input-output theory [33,34] to develop a relation connecting in-going and out-going mode amplitudes via the S-matrix. Proceeding with this analysis results in the usual input-output relation, where f in (ω) and f out (ω) are the in-going and out-going single photon mode amplitudes in frequency space that reference either the initial or final time t s defined by: The single photon S-matrix is found by solution in the frequency domain for the coefficients in our ansatz, |Ψ 0 . In frequency space we have: where the inverse susceptibility of the resonators combined with the central circuit is given by: The elements of the term T (ω) are defined relative to starting in the ground state of the central circuit: where we drop the tilde to indicate the complex-valued Fourier transform of the voltage-voltage linear response function given by Eq. (10) in Section II D. Having established the constituent pieces of this single photon Smatrix, we can now move forward to determine the necessary conditions for the desired ideal circulation.

C. Conditions for ideal circulation
Having established the S-matrix using input-output theory, we would like to determine what conditions exist on parameters of the circuit such that this matrix matches that of an ideal circulator's S-matrix: Examination of Eqs. (18) and (19) indicates that any nonreciprocal behavior must originate from the circuitspecific term in the S-matrix, T (ω), as all other components are trivially diagonal. Focusing on this coupling matrix, we will now work to find conditions for nonreciprocity. Since the central circuit of our model is symmetric with respect to cyclic permutation, the diagonal elements of T (ω) are all equal, as are the off-diagonal elements that are cyclic permutations of each other, i.e., T 12 = T 23 = T 31 and T 21 = T 32 = T 13 . With this in mind, we rewrite Eq. (20) as: with α = iT jj+1 and β = iT jj . We note that this problem can be solved exactly-doing so reveals that optimal performance will occur when frequency shifts of the resonators due to the central circuit are compensated, and when the coupling to the central circuit is matched to the resonator decay. Explicitly, for a signal at target frequency ω T , these conditions are and Rewriting S(ω T ) for these optimal values, we have the following matrix equation: This shows that near resonance (ω ≈ ω T ), the S-matrix will be dominated by the behavior of T (ω). The first terms that contribute to nonreciprocal behavior can be understood by expanding (1 + iηT (ω)) −1 in η < |T |.
We see terms in the series that correspond to clockwise or counterclockwise circulation beginning and ending at the same port, namely T 21 T 32 T 13 and T 12 T 23 T 31 , respectively. If the central circuit is reciprocal, the quantity |T 12 T 23 T 31 − T 21 T 32 T 13 |, which represents the destructive interference between a signal going clockwise and counterclockwise around the circuit, should be zero. However, if the system is ideally nonreciprocal, they instead constructively interfere. Writing α = |α|e i∆θ/6 , where the phase is written in terms of the difference accumulated between a clockwise and counterclockwise path ∆θ, we have |T 12 T 23 T 31 − T 21 T 32 T 13 | = 2|α| 3 | sin(∆θ/2)|. Thus for optimal behavior, We now have three conditions on our system that, when satisfied, yield ideal circulation in either the clockwise or counterclockwise direction: With these conditions, as long as the central circuit has a target frequency for which the phase difference condition in Eq. (27a) holds, we may choose one of the external circuit parameters and let the conditions given in Eqs. (27b) or (27c) determine the others. Note that these results are consistent with previous work, e.g., Ref. [35]. Additionally, these results specify the phase differences corresponding to clockwise (CW) or counterclockwise (CCW) transmission, yielding: ∆θ = −π, 3π, −5π... , CW Transmission π, −3π, 5π... , CCW Transmission .
We now have all we need to explore our system.

IV. NUMERICAL RESULTS
We proceed now with a numerical analysis, examining the input-output properties of the entire system using the analytical results established in Sections III B and III C. Broadly speaking, we wish to enforce the conditions for ideal circulation given in Eq. (27) and then examine the S-matrix given altogether by Eqs. (18), (19), and (20). As there are a variety of different regimes of the central circuit, we investigated several values of external flux A, described in Appendix B. The best performing region was identified near A ≈ 0.2478, and we focus on this in what follows.
To evaluate performance, we compute the elements of T (ω) [Eq. (20)] using the energy eigenstates and eigenvalues established numerically in Section II C. With the elements of T (ω), we find the phase difference ∆θ as a function of potential incoming signal frequency ω. Examining the phase difference allows us to identify target frequencies where ∆θ = ±π, ±3π..., and use Eqs. (27b) and (27c) to set external circuit parameters that fix κ, ω R , and g. We can then examine relevant S-matrix properties, as well as introduce variation in the offset charge on each island and flux penetration per loop to test the system's susceptibility to charge and flux disorder, i.e., variations of these parameters on timescales well below the circuit's bandwidth.
There are numerous points for which the phase difference condition is satisfied, as indicated in Fig. 4c. Of particular interest are instances like the third and fourth points, found at ∆θ = −π where ω T3 ≈ 4 GHz and at ∆θ = −π where ω T4 ≈ 5.46 GHz, respectively. At the third point, the first derivative of the phase difference is zero, indicating much more stability with regard to fre- quency than at other target frequencies. At other points, the first derivative is much larger and a small change in frequency results in significantly larger changes to ∆θ. The fourth point is situated between the first and second excited energies of the central circuit. Due to the close proximity of this target frequency to two resonances of the circuit, it may be easier to inject energy into the circuit here than at other nonreciprocal points. In both cases, intuition suggests that these target frequencies, ω T3 and ω T4 , may allow for increased bandwidth of the system. This possibility will be explored upon examination of elements of the S-matrix in Section IV A.
A. Bandwidth, nonreciprocity, and circulation We now determine performance via the S-matrix property that the normalized output power at the i th port provided input at the j th port is given by |S ij | 2 . Setting the external circuit parameters g, κ, and ω R for the center frequencies ω Ti with i = 1, 2, 3, 4, the target frequencies identified in Fig. 4c, we fix g to a reasonable but large value of g = 1.6 GHz and use Eq. (27) to determine κ and ω R . Table I shows the resulting parameters for the first four target frequencies in Fig. 4c. We remark that for large bandwidth operation, g needs to be large, but making stronger coupling to the resonator has limits in other relevant parameters, as discussed in Appendix C.
For each set of external circuit parameters specified by the target frequencies in Table I, we wish to examine the bandwidth of the resulting transmission between two ports. For each case, the direction of maximum transmission is identified by the phase difference in Fig 4c, and we plot the normalized output power (in dB), displaced by the target frequency of each for comparison. These results are shown in Fig. 5a.
We define the bandwidth as a 3 dB window: the width of the signal in frequency space over which transmission is greater than 50 %. Examination of Fig. 5a indicates that the parameters for the first and second points result in bandwidths of approximately 23 MHz and 17 MHz, respectively. At the third and fourth points, however, we see much larger bandwidths. At the third point we find a bandwidth of approximately 144 MHz, while the parameter set for the fourth point yields a very large bandwidth of approximately 1.34 GHz. These considerably larger bandwidths appear to confirm the intuition described previously: features in the phase difference like that at the third point in Fig. 4c can lead to increased bandwidth, but even larger bandwidth is achieved when nonreciprocal points are situated near resonances-e.g., the fourth point-providing greater ease in inserting energy into the system (that is, a large |α|) and a larger range of frequencies over which it is possible to do so. Since this large bandwidth is desirable in systems such as these, we proceed now using the external circuit parameters for the fourth point (at ∆θ = −π where ω T4 ≈ 5.46 GHz).
Upon fixing external circuit parameters to maximize bandwidth, we examine nonreciprocal behavior and ideal circulation in the system. To do so, we compare transmission from a single port to each of the three outputs, plotting the normalized output powers (in dB) at each port provided input at a single port to observe any circulator behavior. Fig. 5b shows the resulting transmission to each port provided input at port one.
We see in Fig. 5b that microwave signals only propagate in one direction; there is a clear preferred direction of signal transmission, with full transmission at ω ≈ 5.46 GHz in the clockwise direction. Fig. 5b demonstrates the system is lossless and nonreciprocal, and at ω ≈ 5.46 GHz, matched-all characteristics of an ideal circulator's S-matrix [10]. Note that due to the threefold rotational symmetry of the system, these results are reproduced under all cyclic permutations.

B. Resilience to charge and flux noise
Thus far, our analysis has assumed high symmetry: identical, static offset charges N gi on each island due to the presence of applied voltages V gi as well as identical loop areas in the central circuit (ensuring identical flux penetration per loop). In particular, given our expectation that charge offset is not a strong effect, we arbitrarily set N g1 = N g2 = N g3 ≈ 0.4 Cooper pairs, A ≈ 0.2478, and subsequently tuned the external circuit parameters for ideal circulation in the prior parts of this section. However, sources of noise that break this high symmetry are inevitable. Each island of the superconducting circuit may be subject to random variation in the offset charge that will necessarily affect performance given the intermediate energy regime in which we have chosen to work. Similarly, static flux variations across the device due to imperfections in manufacture of the central circuit or local impurities, unexpected magnetic field gradients, and other phenomena will likely degrade the ideal circulation we have observed in Section IV A. Therefore, we wish to address to what degree this is the case by introducing these sources of noise into our analysis. We note that an additional source of disorder is that of the Josephson junctions themselves; we examine this fabrication disorder in Appendix D 2.
To test the system's susceptibility to local charge noise, we begin with the system tuned for ideal circulator be-  Table I. For the first, second, and fourth points, transmission direction is such that |Si,i−1| 2 is plotted. For the second point, |Si,i+1| 2 is plotted. (b) A comparison of normalized output powers at each port assuming input at port one. At ω ≈ 5.46 GHz we observe full transmission in a clockwise fashion to the second port with zero transmission in the opposite direction, indicative of circulator behavior. Parameters and offset charge used are the same as in Fig. 4. External circuit parameters are set according to the fourth point, shown in Table I. The median of ER Υ 0 (defined in Appendix D 1 as the extinction ratio for devices with less than 1 dB of insertion loss), ER Υ 0 , is plotted as a function of the standard deviation of (c) the charge noise distribution Nσ and (d) Rσ, the distribution of fractional flux penetration per loop. Error bars indicate the spread from the first to third quartile of ER Υ 0 for each value of Nσ and Rσ. Device yield Υ (defined in Appendix D 1 as the fraction of devices with less than 1 dB of insertion loss) is plotted as a function of (e) Nσ and (f) Rσ. Parameters for ideal tuning in plots (c)-(f) are the same as in (b).
havior with equal, fixed offset charges N gi . We then introduce an additional offset charge ∆ i that is unique for each island and randomly sampled from a normal distribution with a mean of N µ = 0 and a standard deviation N σ in units of Cooper pair number. As detailed and defined in Appendix D, we assess performance in the presence of charge noise by examining ER Υ 0 and Υ as a function of N σ , shown in Figs. 5c and 5e, respectively. These quantities convey the distribution of the minimum of the extinction ratio for those noise samples resulting in less than 1 dB of insertion loss and the likelihood that this is the case using a sample size of 1200 noise samples.
Upon increasing N σ , we observe an initial slight decrease in performance indicated by the increase in ER Υ 0 (the median of the distribution of ER Υ 0 ), as shown in Fig. 5c. However, at N σ ≈ 0.25 Cooper pairs we observe saturation to ER Υ 0 ≈ -15 dB. This saturation is reasonable given that the behavior of the circuit is unchanged upon rescaling by any integer number of Cooper pairs. Once N σ exceeds ≈ 0.25, the distribution is wide enough to encompass this rescaling. The device yield Υ displays similar behavior with increasing N σ , as shown in Fig. 5e, where the likelihood of large transmission also saturates at ≈ 0.25 to about 50 %. It is very promising that even in the presence of essentially arbitrary, static charge noise, the system retains its nonreciprocity.
To test the system's susceptibility to flux noise, we begin with the system tuned for ideal circulator behav-ior where all three loops of the central circuit are equal in area and penetrated by a flux A ≈ 0.2478. We then introduce variation in the flux penetration of each loop, ∆ i , that is unique for each loop but maintains the total area and flux penetration of the ideal case, i.e., ∆ 1 + ∆ 2 + ∆ 3 = 0, to account for tuning the flux optimally for a given device (where the total field can be varied but not, perhaps, the individual loop fields). We take ∆ i to be a fraction of A, the original flux penetration per loop, randomly sampling from a normal distribution with a mean of R µ = 0 and a standard deviation R σ . As detailed and defined in Appendix D, we assess performance in the presence of flux noise by examining ER Υ 0 and Υ as a function of R σ , shown in Figs. 5d and 5f, respectively. As stated above, ER Υ 0 and Υ convey the distribution of the minimum of the extinction ratio for those noise samples resulting in less than 1 dB of insertion loss and the likelihood that this is the case using a sample size of 1200 noise samples.
Upon increasing R σ , we observe a decrease in performance indicated by the increase in ER Υ 0 (the median of the distribution of ER Υ 0 ), as shown in Fig. 5d. This degradation in performance is more substantial than in the charge noise case-decreasing more rapidly and saturating at R σ ≈ 0.025 to ER Υ 0 ≈ -1 dB. Note that this flux noise model encompasses both random variations in loop area, e.g., imperfections in manufacture, as well as static flux variations across the device that may occur. We do not consider variation in total flux penetration, as this quantity can be experimentally calibrated.

V. CONCLUSIONS
In this work, we have presented a theoretical analysis of a compact superconducting circuit operating as a circulator. Our investigation shows that operating in a non-charge-conserved regime with mobile vortices in the XY model, where the charging energy and Josephson en-ergy scales are similar, leads to performance with a reduced dependence on static charge noise. Furthermore, the presence of both clockwise and counterclockwise vortex circulation at nearby energies appears to enable wide bandwidth performance that exceeds one GHz.
However, a variety of open questions remain for further investigation. It is likely possible to generalize the understanding achieved in this three-port model to a circulator with greater than three ports, e.g., a doublejunction circulator, which could be explored in future work. In addition, in the context of quantum computing and hardware, circulators are typically needed to handle both lower power signal routing as well as higher power device isolation. In this work, our single photon analysis suggests the model is suitable for such low power signal routing. However, a crucial operational question is the performance of the system at higher powers, beyond the single photon limit, thereby expanding the model's potential use cases. We anticipate that further exploration of the vortex-tunneling picture may elucidate regimes of operation where photon-vortex-photon terms enable high power operation, but may also require more complex circuit designs. Another key question is whether ordered vortex states, such as those arising in larger arrays, can exhibit robust edge modes for vortex circulation that may further reduce the effects of noise and disorder in circulator performance. Additionally, implementation of this type of system in an experimental setting will help improve our understanding of quantum vortex dynamics and may lead to other new avenues of exploration (e.g., Ref. [36]) for these and other topological excitations. The implementation of the grounded node present in the proposed central circuit (as shown in Fig. 1b) is a difficult experimental task. A simple modification-floating the central node-provides an alternative circuit that is easier to realize. In this appendix, we explore this alternative circuit (as shown in Fig. 1c), highlighting the connection between this model and the original and showing that similar behavior can be achieved without the presence of a central grounded node.
By replacing the central ground with a floating island, we introduce two additional degrees of freedom-the charge number and phase on this new island,N 0 and ϕ 0 , respectively. To achieve similarity with the original model, the new central island should mimic a ground node, which we can accomplish by increasing the island's capacitance via its physical size. However, in doing so, we can no longer neglect the self-capacitance of the now large central island as well as the associated mutual capacitance it has with the outer three islands.
In all, the modified Hamiltonian is given by: where nowN = {N 1 ,N 2 ,N 3 ,N 0 } and N g = {N g1 , N g2 , N g3 , 0}. The capacitance matrix C is now a 4x4 matrix with the form: where C S = C C + C G + C µ + 3C J and C F = C J + C µ with C 0 the self-capacitance of the center island, C µ the mutual capacitance between it and the outer islands, and all other capacitances as defined in Section II A of the main text. We express the eigenvalues of the capacitance matrix in Eq.(A2) as and define the quantity λ = C 0 C J +(C J +C µ )(4C J +C µ ) so that we may write the inverse capacitance matrix C −1 in the compact form: where δ = λ λ3λ4 , η = λ1(CJ+Cµ) λ3λ4 , and ξ = λ1(λ1−3CJ) λ3λ4 . Notably, by replacing the central ground with a floating island, this alternative central circuit is now chargeconserving. We therefore define the total charge operator N T =N 0 + 3 i=1N i , which commutes with the Hamiltonian [H,N T ] = 0 and whose expectation value is constant, yielding the total charge in the system N T = n. With this conserved quantity, we can eliminate a dynamical degree of freedom.
In particular, by considering the unitary transforma-tionÛ we can eliminate the charge number and phase operators of the central island from the Hamiltonian. Under this transformation, the charge number operators of the outer islands are unaffected while the charge number operator of the central island transforms asN 0 →N 0 − 3 i=1N i . Similarly, the phase operators of the outer islands transform asφ i →φ i +φ 0 while that of the central island is unchanged. Note that the total charge operator transforms asN T →N 0 , i.e., in this new frameN 0 represents total charge and is the conserved quantity.
Performing the unitary transformation in Eq.(A5) on the Hamiltonian in Eq.(A1), we find that the Hamiltonian in this frame takes a form The form of the Hamiltonian in this frame is identical to that of our original model [H CC given by Eq.(1)], where now the charging energy of the circuit E Σ , its dimensionless inverse capacitance matrix K −1 , and the offset charges N gi all take on slightly modified forms. E Σ is defined as , where λ 1 = C C + C G + C µ + 4C J plays the role of C M from the main text but modified by the presence of C µ . We find that K −1 takes a familiar form where ρ = δ − 2η + ξ and we see upon comparison to Eq.(3), we just have δ → ρ. Finally, the new offset charges are uniformly shifted from their original value: Note that in the limit where C 0 → ∞ and C µ → 0, these modified components revert to the forms specific to H CC (the original model's Hamiltonian). However, provided a relatively large C 0 and small C µ (e.g., C 0 = 30C J and C µ = 0.25C J suffice), we find similar transmission performance in a zero-noise setting as compared with that of the original circuit with a central ground. Importantly, we see that aside from the more straightforward parameter-related modifications to E Σ and K −1 , the presence of this new floating central island is equivalent to our original model with the addition of an offset charge that acts symmetrically on each island.
In Section IV B of the main text, we examined the effects of random charge disorder on performance by introducing a shift in offset charge that was randomly sampled and unique to each island. Having uncovered the equivalence between our original model and this alternative circuit, we can examine the effect of replacing the central ground with a floating island by instead looking at uniform charge disorder, i.e., introducing a shift in offset charge that is randomly sampled but identical for each island, in the original circuit with a central ground. This is shown in Fig. 6 without implementing any transmission threshold (as described Appendix D 1).
Regarding performance, we see in Fig. 6 a saturation to ER 0 ≈ −10 dB at N σ ≈ 0.25. We can understand the increase in N σ displayed on the x-axis in Fig. 6 as an increase in total charge N 0 , since the uniform shift we introduce is linearly dependent on N 0 as shown by Eq. (A8). Generally, this total charge can be taken to be a constant scalar quantity, as it is typically fixed to some value in the process of cooling down the system in an experimental setting. However, these results suggest continued performance even with the addition of a central island and some variation in N 0 . While it may be possible in an experimental setting to compensate for the shift in offset charge due to the additional island, it is promising that without any intervention, we find good performance.  Figs. 7a and 8a). In contrast to the large bandwidth circulator performance found at A ≈ 0.2478, we find that the chosen points in the first and third regions do not exhibit similar qualities. We understand this behavior as follows.
In region 1, we note that the second and third excited states are nearly degenerate while the first excited state maintains a nearly constant separation from the ground state, mirroring the ground state, as evidenced in Figs. 2a and 2b of the main text. Furthermore, in this region there are no circulating vortices between the ground and first excited states, in contrast to that with the second excited state, shown in Figs. 3d and 3e. Since the expected mode of circulator operation is via the tunneling of persistent current vortices, points for ideal transmission should not occur around the first excited state. This is confirmed by examining the phase difference at A ≈ 0.17, shown in Figs. 7b and 7c, where ideal nonreciprocal points do not appear until the second and third excited states. This is understood by noting that between the ground and both the second and third excited states there exists circulating and counter-circulating vortices, analogous to the vortex dynamics of the first and second excited states in region 2. We also note that at A ≈ 0.17, the energy level structure of the first, second, and third excited states (highlighted in Fig. 7a) is similar to that of the ground, first excited, and second excited states at A ≈ 0.2478 in region 2 (as investigated in the main text). However, the separation of the second and third excited states at A ≈ 0.17 is about half the size of the separation between the first and second excited states at A ≈ 0.2478. Examining performance at A ≈ 0.17, we find transmission is relatively large around the second and third excited states, shown in Figs. 7d and 7e, as compared to transmission at higher energy nonreciprocal points (not shown here). However, due to the closer proximity of the second and third excited states (as compared to the first and second excited states at A ≈ 0.2478), we find smaller bandwidth when compared to that found at A ≈ 0.2478.
In region 3, the energy level structure is markedly different than either of the other two regions. As shown in Fig. 8a, the level separations between the ground state, first excited state, and second excited state are roughly equivalent, and higher energy states are overall more accessible. Vortices are present between the ground and both the first and second excited states (as well as others). In fact, the first and second excited states exhibit  ) and (c), where the procedure for determining g, κ, and ωR is identical to that described Section IV A in the main text.
circulating and counter-circulating vortices that are qualitatively similar to that of the second region. For this reason, numerous points for ideal circulation around these resonances should occur, which is confirmed by the phase difference found at A ≈ 0.36, shown in Fig. 8b. However, despite the abundance of vortex dynamics and overall lower energy as compared to either the first or second regions, the resulting performance-shown in Figs. 8c, 8d, and 8e-is not as good as that found in region 2. This is perhaps due to the larger separation of energy levels, in particular between the first and second excited states, in contrast to the first and second excited states in region 2. Parameters are as in (b), where the procedure for determining g, κ, and ωR is identical to that described Section IV A in the main text.

Appendix C: Parameter Considerations
In the model we have considered, shown in Figs. 1a and 1b, there are a variety of central circuit and external circuit parameters that can affect system behavior.
Unfortunately, not all of these parameters have clear analytic expressions that characterize their influence on device performance. When considering performance metrics like transmission bandwidth and noise sensitivity, we must strike a balance between these parameters to find optimum performance. In this appendix, we discuss some of these parameters and our understanding of their effect on the system's performance.
Regarding the central superconducting circuit, of tremendous influence is the ratio of Josephson energy E J to charging energies E C and E Σ . As touched on in Section II B, the relative size of these energies impacts the relevant excitations and their quantum behavior, however, it also impacts the circuit's noise sensitivity. As E J increases with fixed plasma frequency (E C , E Σ decreasing together), the circuit becomes more robust against charge noise, but more susceptible to flux noise. The opposite is true as E J decreases. Additionally, our investigation suggests that this ratio impacts the role of fabrication disorder, as discussed in Appendix D 2, whereby increasing E J leads to more susceptibility to noise in E J and less susceptibility to noise in various capacitances. Another important feature is the effect E J has on the separation of adjacent energy levels of the central circuit. As indicated by our investigations in Appendix B, we believe that the proximity of relevant energy levels impacts the bandwidth achieved in each region of flux. As E J increases with fixed plasma frequency, the separation of the circuit's energy levels decreases, limiting bandwidth performance. Considering all of these effects together, our goal with the results presented in the main text, where E J /E C = 4 and E J /E Σ ≈ 40, is to balance the benefits of increased bandwidth with the detriment of decreased resistance to charge noise, all while remaining in a regime where the quantum behavior of vortices may be exploited.
In addition to central circuit parameters, the parameters that characterize the connected external circuits (g, κ, ω R ) are also quite important, influencing bandwidth performance as well as placing physical constraints on the elements of the system. Generally speaking, increasing the size of either of the coupling strengths g or κ will provide some increase in bandwidth. However, doing so affects physically relevant quantities such as the capacitance and characteristic impedance of the resonators (C R ,Z R ) and transmission lines (C TL ,Z TL ). A microscopic derivation utilizing the weak coupling limit (C C C R , C TL ), as done in Refs. [19,37], yields expressions for both g and κ, the coupling strength between the circuit and resonators and the coupling strength between the resonators and transmission lines, respectively. For g we find where ω R is the resonator frequency, Z R the characteristic impedance of the resonators, and R K the von Klitzing constant. The capacitances C C and C M are as defined in Section II A of the main text. Increasing the size of g (and necessarily ω R ) increases Z R while pushing the ratio C R /C C down, impacting the validity of the weak coupling limit. When combined with the conditions for ideal circulation given by Eq. (27) and the resulting parameter specifications given in Table I, Eq. (C1) determines the characteristic impedance of the resonators required to achieve ideal circulation for each operating point. In the results discussed in Section IV A and shown in Fig. 5, the coupling strength was chosen to be g = 1.6 GHz for all operating points. For the fourth operating point, which has particularly large bandwidth, this value of g along with the other parameters given in the main text can be used in Eq. (C1) to give a value of Z R ≈ 700 Ωlarger than that of a typical coplanar waveguide, but not unreasonable [38][39][40].
In addition, we find for κ the expression where ω T is the target frequency, Z TL the characteristic impedance of the transmission lines, and C R the capacitance of the resonators. Increasing κ, given the required conditions for ideal circulation as provided in Eq. (27) of the main text, necessarily increases g, resulting in the same impact on Z R and C R /C C as described above. Therefore, while gains in bandwidth may be possible by increasing g or κ, it must be weighed against the resulting increases in these physically constraining parameters.
We can also use Eq. (C2) and the parameter values shown in Table I of the main text to obtain the required transmission line impedance for an increased g. Unsurprisingly, given the size of κ needed for a larger g, we also need the line impedance tuned to enable this strong coupling, resulting in Z TL ≈ 5 kΩ. While this is very large, this is a microwave engineering challenge that may have various solutions, such as inductively coupling, rather than capacitively coupling, the transmission lines to the resonators. Despite this, given the numerous parameters (and the necessarily vast parameter space), there may exist as-yet-undiscovered operating points that maintain or improve upon our reported bandwidth and noise performance but decrease these large impedance requirements. Our investigation indicates that the presence of nonreciprocal behavior relies on the tuning of parameters relevant to the central superconducting circuit. While the external circuit parameters place physical constraints on the system, we still find that nonreciprocity can be achieved with any reasonable set of g, κ, and ω R . We note that while in this work, we opted to first fix g and then examine requirements of the resonators, one could alternatively first choose the properties of the resonators and examine the resulting implications on g and other system parameters.
This metric provides a comparison between the desired direction of transmission and the undesired, reverse direction of transmission. For an ideal circulator, ER → −∞, thus in practice, the more negative this quantity, the more nonreciprocal the device. The minimum of the extinction ratio ER 0 (as a function of frequency) is what we will use to assess the system's response to noise. Similarly, we define the insertion loss IL between two ports of the circuit as The insertion loss quantifies the amount of transmission lost between two ports of the circuit. For zero transmission, IL → −∞. However, in an ideal circulator there is full transmission between two ports of the system, resulting in IL = 0. To help characterize the system's response to noise, we will examine the maximum of the insertion loss IL 0 , where the closer this value is to zero, the better the performance. The introduction of deviations from ideal parameters (noise) will necessarily lead to variations in device performance. These deviations parameterize a multidimensional landscape of device performance, whose dimension is characterized by the degrees of freedom associated with each type of disorder. For example, variations of the static gate voltages on each of the three islands (here translated to effective Cooper pair number) define a three-dimensional landscape characterizing performance for variations from the ideal offset charge. To showcase this landscape, we consider a two-dimensional slice that examines the variations in two of the three charge noise parameters, shown in Fig. 9. A landscape of this type exists for flux noise (where it is a three-dimensional space) as well as for fabrication disorder (where it is a six-dimensional space).
In what follows, rather than attempting to tune a given device to a perfect operating point, we instead consider a statistical analysis that samples from all dimensions of each landscape to find the fraction of devices with low insertion loss, and of those devices, the distribution of their minimum extinction ratios, specifically highlighting a median device outcome. This would correspond, in the case of static disorder, to finding high performing devices, and provides a yield of viable devices given a performance metric.
FIG. 9. The effect of offset charge number variations from the ideal on the minimum of the extinction ratio ER0 (in dB) between ports one and two. Shown is the two-dimensional slice for ∆3 = 0, that is, an unchanged Ng3, while introducing variations ∆1 and ∆2 from −1 to +1 in Ng1 and Ng2, respectively. Parameters and offset charge for ideal tuning are the same as in Fig. 4.

Statistical Analysis
Upon introducing either charge noise or flux noise (as described in the main text), sampled from a normal distribution with standard deviation N σ or R σ , respectively, we examine both ER 0 and IL 0 for a total of 1200 noise samples. As an example, Fig. 10a shows the pairwise plot of ER 0 and IL 0 for 1200 charge noise samples for which N σ = 0.25. Similarly, Fig. 10b shows the same pairwise plot but for 1200 flux noise samples for which R σ ≈ 0.012.
We then impose a performance threshold, examining only those noise samples for which there is less than 1 dB of insertion loss (IL ≥ −1), shown visually by the points within the shaded regions of Fig. 10. We define this subset of points as ER Υ 0 , characterizing this distribution by calculating its quartiles for increasing N σ or R σ . Figs. 5c and 5d of the main text display the median of this distribution ER Υ 0 with error bars from the first to third quartiles as a function of N σ and R σ , respectively.
We also define the device yield Υ as the probability for transmission to exceed the imposed threshold (IL ≥ −1). This is approximated by Υ = # of Samples with IL ≥ −1 Total # of Samples .
Figs. 5e and 5f display this probability for increasing N σ and R σ , respectively. Note that this procedure is also utilized for the fabrication disorder introduced and described below in Appendix D 2.  Fig. 4. Results shown are specifically for transmission between ports one and two, but are representative of the circuit overall due to circuit symmetry and the sampling from identical noise distributions at each island.

Fabrication Disorder
In addition to considering charge and flux disorder, we also examine the effects of fabrication disorder in the individual Josephson junctions of the model, shown in Fig. 1b. The high symmetry present in our analysis also extends to the junctions-we assume identical Josephson energies E J and junction capacitances C J (yielding identical charging energies E C = e 2 2CJ ) for all six junctions. By separately introducing disorder in these quantities, we can then assess the model's subsequent performance when subject to variations in junction parameters.
To test the system's fabrication tolerances, we begin with the system tuned for ideal circulator behavior with all six Josephson junctions having identical E J and C J . Analogous to the charge and flux noise cases, we then introduce a variation ∆ in either E J or C J that is unique for each junction. In both the E J and C J noise cases, we take ∆ to be a fraction of the original value, randomly sampling from a normal distribution with a mean of D µ = 0 and a standard deviation D σ . As detailed and defined above in Appendix D 1, we assess performance in the presence of either E J or C J noise by examining ER Υ 0 and Υ as a function of D σ , shown together in Figs. 11a and 11b, respectively. As stated above, ER Υ 0 and Υ convey the distribution of the minimum of the extinction ratio for those noise samples resulting in in less than 1 dB of insertion loss and the likelihood that this is the case using a sample size of 1200 noise samples.
Upon increasing D σ , we observe markedly different performance between the introduction of noise in E J versus that in C J , as shown in Fig. 11. In particular, we see that while the introduction of noise in the Josephson energies of the junctions degrades performance, varying the junction capacitances does not. This is best understood by considering that the original ideal parameters are such that E J = 4E C ≈ 40E Σ , thereby necessitating that any variations introduced in E J are bound to have a larger effect than those in C J (and subsequently E C and E Σ ). Nonetheless, Fig. 11 would indicate that variations of approximately 1.5 % in E J may be tolerated. Beyond 2 %, nonreciprocity is extinguished, where we see that ER Υ 0 ≈ 0 and Υ 20 %. In contrast, Fig. 11 shows that variations up to 5 % (and likely higher) in C J can be tolerated-for all distribution sizes ER   11. (a) The median of ER Υ 0 (defined in Appendix D 1 as the extinction ratio for devices with less than 1 dB of insertion loss), ER Υ 0 , is plotted as a function of the standard deviation of the EJ or CJ noise distribution Dσ as indicated by the legend. Error bars indicate the spread from the first to third quartile of ER Υ 0 for each value of Dσ. (b) Device yield Υ (defined in Appendix D 1 as the fraction of devices with less than 1 dB of insertion loss) is plotted as a function of Dσ. Parameters for ideal tuning are the same as in Fig. 5b.