Circulation by microwave-induced vortex transport for signal isolation

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 $\Phi$. We note that the central ground can, in principle, be replaced by a capacitive island with $C_{0}\gg C_{\text{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 operators $\underline{\hat{N}}$ and $\underline{\hat{\varphi}}$ are the vectors of number and phase operators for each island of the circuit, with individual elements satisfying the commutation relation $[\hat{\varphi}_{j},\hat{N}_{k}]=i\delta_{jk}$. $E_{\text{J}}$ is the Josephson energy of each junction and $E_{\Sigma}$ is the circuit’s charging energy, defined as $E_{\Sigma}=\frac{e^{2}}{2C_{\text{M}}}$ where $C_{\text{M}}$ is the largest eigenvalue of the circuit’s capacitance matrix $\underline{\underline{C}}$, defined in Eq. (2). $\underline{\underline{K}}^{-1}$ is the dimensionless inverse capacitance matrix appropriate for the circuit, found by factoring $C_{\text{M}}$ from $\underline{\underline{C}}^{-1}$. As it is written, $H_{\text{CC}}$ only stipulates the central circuit be a superconducting, three-island system containing identical Josephson junctions. By specifying the elements of $\underline{\underline{K}}^{-1}$ and the tunneling potential $V(\underline{\hat{\varphi}},\Phi)$ 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 $\underline{\underline{K}}^{-1}$ by inverting the capacitance matrix $\underline{\underline{C}}$,

where $C_{\text{S}}=C_{\text{C}}+C_{\text{G}}+3C_{\text{J}}$ with $C_{\text{J}}$ the capacitance of each Josephson junction, $C_{\text{C}}$ the coupling capacitance, and $C_{\text{G}}$ the residual capacitance to ground. For this matrix, $C_{\text{M}}=C_{\text{C}}+C_{\text{G}}+4C_{\text{J}}$ is the largest eigenvalue. Factoring $C_{\text{M}}$ from $\underline{\underline{C}}^{-1}$ we can then find:

where $\delta=\frac{C_{\text{J}}}{C_{\text{C}}+C_{\text{G}}+C_{\text{J}}}$.

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 $\varphi$ represents the angle of the $U(1)$ degrees of freedom.

The acquired Peierls phase factor, $2\pi 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_{\text{CC}}$ as $\ket{\varepsilon_{m}}$ with eigenenergy $\varepsilon_{m}$, and label the ground state $\ket{G}$. These are parametric functions of $A$, and periodic in $A$.

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_{\text{J}}$, and the energy associated with putting a single electron on an island, $E_{\text{C}}\sim E_{\Sigma}$. In the two extreme energy limits, $E_{\text{J}}\ll E_{\text{C}}$ or $E_{\text{J}}\gg E_{\text{C}}$, the relevant excitations are charges or vortices, respectively [26, 27, 28]. Note that charges have the same properties in the charged-dominated regime that vortices have in the tunneling-dominated regime [29, 27]. In the charge-dominated regime, where $E_{\text{J}}\ll E_{\text{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 [30, 31]. 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_{\text{C}}$.

In the classical tunneling limit, where $E_{\text{J}}\gg E_{\text{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) [29, 27, 28, 32]. 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. [33]). The plasma oscillations do not break time-reversal symmetry; in essence, the Meissner effect traps flux in the vortices, and thus, all time-reversal symmetry-breaking terms arise only from vortex motion. When $E_{\text{J}}\gg E_{\text{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_{\text{J}}\sim E_{\text{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. [20] and M$\ddot{\text{u}}$ller et al. [21]. However, in those works, the system was defined with a charge-conserving 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.

To explore the vortex picture, we set $E_{\text{J}}=4E_{\text{C}}$ to ensure that there are sufficient phase fluctuations to warrant a quantum treatment of the vortex dynamics, and numerically diagonalize $H_{\text{CC}}$ [given by Eq. (1)] for a range of $A$ from 0 to 0.5, as the eigenvalues are symmetric in $A\rightarrow-A$ and periodic (same for $A+1$ as for $A$). We denote the corresponding energy eigenstates $\ket{\varepsilon_{n}}$ and eigenvalues $\varepsilon_{n}$, with implicit dependence on $A$ understood. We note for our chosen parameters ($E_{\text{J}}=4E_{\text{C}}\approx 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 $\ket{n_{1},n_{2},n_{3}}$ for Cooper pairs on each island, with $|n_{i}|\leq 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\approx 0.2$ and $A\approx 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_{\text{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 persistent-current 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 $\{\ket{\varepsilon_{0}}(=\ket{G}),\ket{\varepsilon_{1}},\ket{\varepsilon_{2}}\}$, denoted $I^{kl}_{i}=\bra{\varepsilon_{k}}\hat{I}_{i}\ket{\varepsilon_{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\neq l$.

We write the argument of the $I^{mn}_{i}$ matrix elements as $\phi_{i}^{mn}$. The differences $\phi_{i}^{mn}-\phi_{i+1}^{mn}$ only take the values $\{0,\pm 2\pi/3\}$, as expected by the three-fold symmetry of the system. We now consider their role in circulation. Examining a superposition $c_{m}\ket{\varepsilon_{m}}+c_{n}\ket{\varepsilon_{n}}$ with $m>n$, we see that the $i^{\text{th}}$ loop-current operator has an expectation value at time $t$ of

where $\phi_{i}^{mn}$ is the argument of the off-diagonal matrix element $I_{i}^{mn}$, $\theta_{mn}$ is the argument of $c_{m}^{*}c_{n}$, and $\omega_{mn}=\frac{(E_{m}-E_{n})}{\hbar}$. 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 $\phi_{i}^{mn}-\phi_{i+1}^{mn}=\pm 2\pi/3$ at a frequency $\omega_{mn}$. This is shown schematically in Fig. 3a.

Accordingly, we define the directional function for the difference in $\phi_{i}^{mn}$ variables modulo $2\pi$ 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. $\Theta^{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.

We now consider what happens if an oscillating voltage is applied to the $j^{\text{th}}$ island of the central circuit via a coupling capacitor. We denote the coupling as $\hbar\Omega_{j}(t)\hat{B}_{j}$, with $\underline{\hat{B}}=\underline{\underline{K}}^{-1}\underline{\hat{N}}$ 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 $\Omega_{j}(t)\propto\delta(t)$ as

where we use $P$ as in a dimensional setting this quantity has units of power, scaling as $\frac{2eI_{\text{C}}}{C_{\text{M}}}$, 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 $\omega$ 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 voltage-voltage correlation function,

whose Fourier transform will be identified later as the transmission matrix of the circuit.

Working this out explicitly,

and

where $\gamma\rightarrow 0^{+}$ is taken to ensure causal ordering and $P_{ij}(\omega)=\tilde{P}_{ij}(\omega)+\tilde{P}^{*}_{ij}(-\omega)$. We note that a narrowband signal will preferentially excite states near the signal frequency, as is apparent in $\tilde{P}_{ij}(\omega)$’s denominator.

Examining Eq. (11), we can interpret the linear response behavior as the evaluation of $\langle\hat{I}_{i}\rangle$ for a state of the form $\ket{G}-i\sum_{n}\outerproduct{\varepsilon_{n}}{\varepsilon_{n}}(\sum_{j}\Omega_{j}\hat{B}_{j})\ket{G}$ for $\Omega_{j}\rightarrow 0$. Thus we can conclude that, insofar as $\hat{B}_{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^{\text{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.

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.1, is threaded by an external flux $\Phi$ 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_{\text{TL}}$ [34] describes each transmission line as containing an identical continuous spectrum of modes over frequency $\nu$. $H_{\text{R}}$ [20] is the contribution from the system’s three identical resonators, consisting of a single low-lying mode with a characteristic frequency $\omega_{\text{R}}$. The operators $\hat{a}^{\dagger}_{k}$ and $\hat{a}_{k}$ create or destroy a photon in the $k^{\text{th}}$ resonator and satisfy the usual commutation relation: $[\hat{a}_{j},\hat{a}^{\dagger}_{k}]=\delta_{jk}$. Similarly, $\hat{b}^{\dagger}_{k}(\nu)$ and $\hat{b}_{k}(\nu)$ create or destroy a photon in the $k^{\text{th}}$ transmission line and satisfy the commutation relation: $[\hat{b}_{j}(\nu),\hat{b}^{\dagger}_{k}(\nu^{\prime})]=\delta_{jk}\delta(\nu-\nu^{\prime})$.

$H_{\text{CC-R}}$ and $H_{\text{R-TL}}$ [34, 20] 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 $\kappa$, to be identical across all three interactions, and will later optimize to ensure the best circulator performance. The operator $\hat{B}_{k}$, as defined in Section II.4, is a dimensionless operator corresponding to the dimensionless voltage on the $k^{\text{th}}$ island, defined as the $k^{\text{th}}$ row of $\underline{\hat{B}}\ =\ \underline{\underline{K}}^{-1}\underline{\hat{N}}$.

Lastly, $H_{\text{CC}}$ is as defined in Section II.1 and describes the central superconducting circuit. Here we also included applied voltages $V_{\text{g}k}$ (shown in Fig. 1a) that can lead to charge offset noise $N_{\text{g}k}$ on each island.

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}(\nu)$ are time-dependent complex probability amplitudes for each state, and we abbreviate, with a slight abuse of notation, the ground state $\ket{0,0,0,G,0,0,0}$ as $\ket{G}$, and the excited circuit state with no photons in any resonator or transmission line $\ket{0,0,0,\varepsilon_{n},0,0,0}$ as $\ket{\varepsilon_{n}}$. The summation over $n$ is a summation over only the excited states of the circuit— it excludes the ground state as this state is not in the manifold being considered.

By using the Schr$\ddot{\text{o}}$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 [34, 35] 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_{\text{in}}(\omega)$ and $f_{\text{out}}(\omega)$ are the in-going and out-going single photon mode amplitudes in frequency space that reference either the initial or final time $t_{\text{s}}$ defined by:

The single photon S-matrix is found by solution in the frequency domain for the coefficients in our ansatz, $\ket{\Psi_{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 $\underline{\underline{T}}(\omega)$ 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.4. Having established the constituent pieces of this single photon S-matrix, we can now move forward to determine the necessary conditions for the desired 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 circuit-specific term in the S-matrix, $\underline{\underline{T}}(\omega)$, 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 $\underline{\underline{T}}(\omega)$ 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 $\alpha=iT_{jj+1}$ and $\beta=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 $\omega_{\text{T}}$, these conditions are

and

Rewriting $\underline{\underline{S}}(\omega_{\text{T}})$ for these optimal values, we have the following matrix equation:

This shows that near resonance ($\omega\approx\omega_{\text{T}}$), the S-matrix will be dominated by the behavior of $\underline{\underline{T}}(\omega)$. The first terms that contribute to nonreciprocal behavior can be understood by expanding $(\mathbb{1}+i\eta\underline{\underline{T}}(\omega))^{-1}$ in $\eta<|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 $\alpha=|\alpha|e^{i\Delta\theta/6}$, where the phase is written in terms of the difference accumulated between a clockwise and counterclockwise path $\Delta\theta$, we have $|T_{12}T_{23}T_{31}\ -\ T_{21}T_{32}T_{13}|=2|\alpha|^{3}|\sin(\Delta\theta/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. [36]. Additionally, these results specify the phase differences corresponding to clockwise (CW) or counterclockwise (CCW) transmission, yielding:

We now have all we need to explore our system.

We now determine performance via the S-matrix property that the normalized output power at the $i^{\text{th}}$ port provided input at the $j^{\text{th}}$ port is given by $|S_{ij}|^{2}$. Setting the external circuit parameters $g$, $\kappa$, and $\omega_{\text{R}}$ for the center frequencies $\omega_{\text{T}i}$ 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 $\kappa$ and $\omega_{\text{R}}$. Table 1 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 1, 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 $|\alpha|$) 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 $\Delta\theta=-\pi$ where $\omega_{\text{T4}}\approx 5.46\text{ 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 $\omega\approx 5.46\text{ GHz}$ in the clockwise direction. Fig. 5b demonstrates the system is lossless and nonreciprocal, and at $\omega\approx 5.46\text{ GHz}$, matched— all characteristics of an ideal circulator’s S-matrix [10]. Note that due to the three-fold rotational symmetry of the system, these results are reproduced under all cyclic permutations.

Thus far, our analysis has assumed high symmetry: identical, static offset charges $N_{\text{g}i}$ on each island due to the presence of applied voltages $V_{\text{g}i}$ 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_{\text{g}1}=N_{\text{g}2}=N_{\text{g}3}\approx 0.4$ Cooper pairs, $A\approx 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.1. 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 behavior with equal, fixed offset charges $N_{\text{g}i}$. We then introduce an additional offset charge $\Delta_{i}$ that is unique for each island and randomly sampled from a normal distribution with a mean of $N_{\mu}=0$ and a standard deviation $N_{\sigma}$ in units of Cooper pair number. As detailed and defined in Appendix D, we assess performance in the presence of charge noise by examining $\text{ER}^{\Upsilon}_{\text{0}}$ and $\Upsilon$ as a function of $N_{\sigma}$, 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_{\sigma}$, we observe an initial slight decrease in performance indicated by the increase in $\widetilde{\text{ER}}^{\Upsilon}_{\text{0}}$ (the median of the distribution of $\text{ER}^{\Upsilon}_{\text{0}}$), as shown in Fig. 5c. However, at $N_{\sigma}\approx 0.25$ Cooper pairs we observe saturation to $\widetilde{\text{ER}}^{\Upsilon}_{\text{0}}\approx$ -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_{\sigma}$ exceeds $\approx 0.25$, the distribution is wide enough to encompass this rescaling. The device yield $\Upsilon$ displays similar behavior with increasing $N_{\sigma}$, as shown in Fig. 5e, where the likelihood of large transmission also saturates at $\approx 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 behavior where all three loops of the central circuit are equal in area and penetrated by a flux $A\approx 0.2478$. We then introduce variation in the flux penetration of each loop, $\Delta_{i}^{\prime}$, that is unique for each loop but maintains the total area and flux penetration of the ideal case, i.e., $\Delta_{1}^{\prime}+\Delta_{2}^{\prime}+\Delta_{3}^{\prime}=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 $\Delta_{i}^{\prime}$ to be a fraction of $A$, the original flux penetration per loop, randomly sampling from a normal distribution with a mean of $R_{\mu}=0$ and a standard deviation $R_{\sigma}$. As detailed and defined in Appendix D, we assess performance in the presence of flux noise by examining $\text{ER}^{\Upsilon}_{\text{0}}$ and $\Upsilon$ as a function of $R_{\sigma}$, shown in Figs. 5d and 5f, respectively. As stated above, $\text{ER}^{\Upsilon}_{\text{0}}$ and $\Upsilon$ 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_{\sigma}$, we observe a decrease in performance indicated by the increase in $\widetilde{\text{ER}}^{\Upsilon}_{\text{0}}$ (the median of the distribution of $\text{ER}^{\Upsilon}_{\text{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_{\sigma}\approx 0.025$ to $\widetilde{\text{ER}}^{\Upsilon}_{\text{0}}\approx$ -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.