Quantum-limited directional amplifiers with optomechanics

Directional amplifiers are an important resource in quantum information processing, as they protect sensitive quantum systems from excess noise. Here, we propose an implementation of phase-preserving and phase-sensitive directional amplifiers for microwave signals in an electromechanical setup comprising two microwave cavities and two mechanical resonators. We show that both can reach their respective quantum limits on added noise. In the reverse direction, they emit thermal noise stemming from the mechanical resonators and we discuss how this noise can be suppressed, a crucial aspect for technological applications. The isolation bandwidth in both is of the order of the mechanical linewidth divided by the amplitude gain. We derive the bandwidth and gain-bandwidth product for both and find that the phase-sensitive amplifier has an unlimited gain-bandwidth product. Our study represents an important step toward flexible, on-chip integrated nonreciprocal amplifiers of microwave signals.

In recent theoretical work, Ranzani and Aumentado [16,17] analyzed general conditions for nonreciprocity in parametrically coupled systems, and showed that nonreciprocity arises due to dissipation in ancillary modes and multi-path interference. Metelmann and Clerk [18] have shown that any coherent interaction can be made directional by balancing it with a dissipative process. Indeed, this insight led to a demonstration of nonreciprocity using optomechanics in the optical domain [19], and theoretical investigations into minimal implementations of directional amplifiers [20].
While implementing the balance between a direct coherent coupling between the cavities and a dissipative interaction is challenging experimentally, Refs. [25][26][27] have recently studied and demonstrated nonreciprocal transmission between two cavity modes where two mechanical resonators each mediate both coherent and dissipative coupling. Here, building on this concept, we propose directional amplifiers using exclusively optomechanical interactions. Microwave tones on the red and blue sidebands enable so-called beam-splitter and two-mode squeezing interactions (cf. Fig. 1), leading to a total of eight controllable terms in the Hamiltonian. We identify and analyze a simple directional phase-preserving amplifier that uses four tones and a directional phase-sensitive amplifier using six tones. While the gain-bandwidth product of the phase-preserving amplifier is limited to the cavity decay rate, the phase-sensitive amplifier has an unlimited gain-bandwidth product. We show that both amplifiers can reach their quantum limits of a half and zero added quanta, respectively, and emit thermal noise from the mechanical resonators in the reverse direction, a necessary consequence of impedance matching and directionality. We show how the reverse noise can be reduced through additional sideband cooling without interfering with directionality or amplification. Our concrete proposal bridges the gap between previous theoretical studies and experimental realization and therefore represents an important step toward on-chip integrated nonreciprocal amplifiers of microwave signals.
Model.-We consider an optomechanical plaquette, comprising two microwave cavities coupled via two mechanical arXiv:1705.00436v2 [cond-mat.mes-hall] 6 Dec 2017 resonators, as shown in Fig. 1(a). The cavities are driven close to the motional sidebands. After a standard treatment, which includes linearizing the Hamiltonian, neglecting counterrotating terms, and going into a rotating frame [37,38], we arrive at the time-independent Hamiltonian (ℏ = 1) (1) where ( ) is the annihilation operator for the th cavity mode (mechanical resonator), = − 0, and = + 0, are field-enhanced optomechanical coupling strengths, ± is the amplitude of the coherent state produced in cavity due to a pump at frequency c, ± (Ω + ), and 0, are the vacuum optomechanical couplings. Since the couplings , depend on the pumps, their amplitude and phase can be controlled. The interactions are represented in Figs. 1(a), 2(a), 3(a) as red ( ) and blue ( ) lines. Further details can be found in the Supplementary Information (SI), including a discussion about the limits of validity of the rotating-wave approximation.
We have chosen the couplings [Eq. (5)] due to the following reasons. First, an even number of blue and red tones ensures equivalent arms of the circuit. Second, amplification requires blue tones. Third, a directional amplifier with four blue tones cannot be impedance matched to the signal source (cf. SI). Last, swapping hopping and amplifier interactions in one arm of the circuit cannot lead to directional amplification [42].
The condition  2 <  1 ensures that the system does not exceed the parametric instability threshold. In the limit of large gain, we obtain our first main result, the scattering matrix with vanishing reverse gain | 2 → 1 (0)| 2 , but forward gain which can in principle be arbitrarily large, as long as the RWA is valid (cf. SI). At the same time, thermal noise from the mechanical resonators is suppressed by increasing  1 , as is demonstrated in Fig. 2(e), where we plot the noise added to the signal [30,37,43], where we sum over all noise sources, with associated thermal occupation , and scattering amplitude to the second cavity → 2 . Using Eq. (7), and denoting thermal cavity (mechanical) occupations by c, ( m, ), the noise on resonance  DPPA = 1 4 1 m,1 + m,2 + 1 + ( 1 +  2 ) 2 4 1  2 c,2 + 1 2 .
(9) As a result, for large  1 ≳  2 , and vanishing thermal occupation of the cavity input, we reach the quantum limit of half a quantum of added noise,  DPPA → 1∕2 [37,43].
Another important figure of merit is noise emerging from cavity 1, characterized by the output noise spectral density, out 1,out ( ′ )⟩, which we plot in Fig. 2(f). Ultimately, the reason for building directional amplifiers is to reduce this figure. On resonance, DPPA 1,out (0) = ( m,1 + m,2 + 1)∕2. Strategies to reduce this figure are discussed below.
The off-resonance behavior of the DPPA is remarkably rich and depends on the dimensionless quantities ∕Γ m, ,  1 ,  2 . We plot forward gain, reverse gain, added noise, and the noise spectrum at cavity 1 as functions of frequency at cooperativities  1 = {1, 3, 10, 30}, Fig. 2(c-f).  2 is chosen such that when increasing  1 both gain and bandwidth are enhanced.
We show in the SI that for Γ m, = Γ m and = and in Fig. 2(c)]. As the gain gets large and  1 ,  2 dominate all other dimensionless parameters, the bandwidth approaches Γ = ( 1 −  2 )∕ 1 , leading to the gainbandwidth product limit ≡ Γ √  → 2 , independent of ∕Γ m . Close to resonance, the reverse scattering amplitude ∕Γ m (cf. SI), such that the product of isolation bandwidth and gain is Γ m . Since the gain bandwidth is larger than the isolation bandwidth, there is large reverse gain off resonance [cf. Fig. 2(d)], and noise from cavity 2 dominates the noise spectral density at cavity 1 [cf. Fig. 2 With increasing effective mechanical linewidth Γ m (through additional sideband cooling), the isolation bandwidth grows, suppressing reverse gain off resonance (cf. red, dashed curve in Fig. 2 and Ref. [18]). In the SI we calculate how off-resonant terms renormalize the parameters of the DPPA.
Directional phase-sensitive amplifier (DPSA).-We now turn to an implementation of a DPSA, which necessitates six tones. Essentially, we replace the amplifier interaction in the DPPA by a phase-sensitive quantum non-demolition (QND) interaction that couples one quadrature of cavity 2 to only one quadrature of the mechanical resonator [44,45], choosing and the same as for the DPPA [Eq. (5)], illustrated in Fig. 3(a,b). Since the QND interaction requires | 2 | = | 2 |, and we require symmetric amplifier arms, two cooperativities suffice to characterize the six tones. In the DPSA, we need to ensure that the measured mechanical quadratures agree, = arg( 11 m,1 (0) 21 ∕ 21 ) = arg(± 12 m,2 (0) 22 ∕ 22 ), and that the information from both emerges in the same cavity quadrature, = arg( 21 21 ) = arg(± 22 22 ). and determine the quadratures involved in amplification. The two remaining phases are an arbitrary mechanical phase, and the plaquette phase.
While there is no parametric instability of the kind that limits back-action evading measurements [46,47], we show in the SI using a Floquet technique [48,49] that counterrotating terms induce an instability threshold for finite sideband parameter (similar to Ref. [50]), and the RWA is only valid for sideband parameters that are bigger than the cooperativities. This is not out of reach [51], but needs to be taken into account in experimental design.
The isolation, detuning, and impedance-matching conditions coincide with those of the DPPA, and we obtain another central result, the scattering matrix (on resonance) where we defined noise scattering intensity  ≡ The amplifier is phase-sensitive and directional, as only the phase quadrature of the second cavity, 2 , inherits the amplified signal from the phase quadrature of the first cavity, 1 . We calculate the noise added to the signal as before The crucial difference to the DPPA is that the noise stemming from reflection of fluctuations at cavity 2 can also be suppressed, such that in the limit  2 ≫  1 ≫ 1 added noise vanishes.
Backward propagating noise and sideband cooling.-The noise emitted in the reverse direction is of central importance for directional amplifiers. For both DPPA and DPSA, the output noise spectral density of cavity 1 on resonance is out 1,DPSA (0) = out 1,DPPA (0) = ( m,1 + m,2 + 1)∕2. Due to impedance matching and directionality fluctuations incident on the cavities do not appear in 1,out [cf. Eqs. (7) and (11)]. The commutation relations of 1,out then imply i.e, mechanical fluctuations have to appear in the output instead. The lowest possible value for 1,out is 1∕2, attainable for zero thermal noise quanta in the mechanical resonators.
However, even in state-of-the-art dilution refrigerators, the required temperatures are out of reach. One way to mitigate backward noise emission is to add another microwave mode to the setup that can replace the fluctuations in the output of cavity 1, essentially realizing a circulator. Without modifying the theory above, one can either increase the resonance frequency of the mechanical modes, which is mainly a technological challenge, or one could resort to external sideband cooling with an auxiliary mode. The latter can achieve m → 0 [52][53][54], and has the added benefit of enhancing mechanical linewidths [cf. red (dashed) curve in Figs. 2 and 3 and SI]. Whilst this could be done with an additional cavity mode for each resonator, implementing a circuit with four cavity modes coupled to two mechanical resonators is a formidable technical challenge. A problem arises when cooling with only one additional mode, since it can lead to a coupling of the mechanical resonators via the extra cooling mode, thereby changing the topology of the system thus spoiling directionality. This can be mitigated by detuning each pump by several mechanical linewidths (cf. SI), making cooling with only one additional mode feasible.
Conclusion.-We have presented quantum-limited, nonreciprocal amplifiers using an optomechanical plaquette comprising two cavities and intermediate mechanical resonators [25,26]. Such devices carry great promise, since they can be integrated into superconducting circuits and amplify near or at the quantum limit, whilst protecting the signal source. For an introduction to optomechanics and input-output theory, we refer to Refs. [S1, S2]. The optomechanical Hamiltonian for the plaquette is given by (ℏ = 1) 0, is the bare coupling between the th cavity and the th mechanical resonator. Each cavity is driven by up to four tones, at frequencies c, ± (Ω + ), i.e., on the red and blue motional sidebands due to the mechanical resonators, as illustrated in Fig. 1 in the main text. The drives generate coherent states in the cavities, such that the annihilation operators for the cavity modes can be written as a sum of coherent parts and fluctuations ̂ ̂ = − c, where are c-numbers that describe the coherent state, and ̂ is a bosonic operator with zero mean. In the following, we will rename ̂ →̂ for a cleaner notation, and drop the hats. Going into a rotating frame with respect to In the main text, we neglect all time-dependent terms for simplicity. When calibrating pumps and phases, their effect will have to be included. Below, we answer three questions about counterrotating terms. First, do they change the topology of the circuit and therefore make directionality impossible? We find that to a good approximation this is not the case, as we analyze in section "Off-resonant terms". Second, do off-resonant terms change the effective parameters of the circuit? As in previous studies the answer is yes, detailed in Sections "Off-resonant terms" and "Stability". Third, can off-resonant terms lead to an instability? Again the answer is yes, for sufficiently strong driving or sufficiently low sideband parameter, as described in Section "Stability". We also find that the RWA theory is valid for large but finite sideband parameter, after taking into account the modification of effective parameters. We describe the system with input-output theory and quantum Langevin equations [S1-S3] The coupling matrices , defined in Eq. 3, exhibit the symmetry since we are using both annihilation operators and their conjugates. The minus sign in the frequency is due to our choice of Fourier transform, [ ( )] † = † (− ). This is the only symmetry, since we know that there are 8 free complex parameters (the 8 driving amplitudes), and has 4 2 = 16 complex entries.
Taking the modulus, we find that isolation requires (S10) Note that sending Φ → −Φ interchanges the two equations in Eq. (S9), which means that isolation now occurs in the opposite direction. In order to avoid "double isolation", where the forward and backward transmission vanish at the same plaquette phase (as per the previous sentence, this occurs for Φ = 0 or Φ = ), we need 1 ∕Γ m,1 ≠ 2 ∕Γ m,2 . The phase at which isolation is obtained is given in Eq. (6) in the main text. For these choices, the scattering amplitude from 1,in to 1,out (reflection) is given by Impedance matching is attained at zero reflection, namely when the same as in Ref. [S4].

DIRECTIONAL PHASE-PRESERVING AMPLIFIER WITH ONLY BLUE TONES
Here we analyze the optomechanical plaquette with only pumps on the upper motional sidebands. While directional phasepreserving amplification is still possible, the signal cannot be impedance matched. We choose coupling amplitudes as follows The analysis proceeds similarly to the DPPA in the main text. Φ is the only physically relevant phase. We choose the same detuning parameterization as before 1 = Γ m,1 and 2 = − Γ m,2 , in which case the plaquette phase takes the same form as before [Eq. (6)]. However, if we look at the optical scattering matrix, we find (S14) The cooperativities are obey 4 2 + 1 >  > 0, where the first condition is required for stability and the second arises by definition, such that impedance matching is not possible, and input will be reflected and amplified. This property is already highly undesirable in a directional amplifier, since the main point of a directional amplifier is protection of the system that emits the signal.
Amplification is obtained if either or both of the cooperativities approach 2 2 + 1∕2.  1 → 2 2 + 1∕2 leads to a lot of noise being emitted from cavity 1, both due to the reflection and due to amplified mechanical noise (not shown above). Hence, we consider  2 = 2 2 + 1∕2 − , for some small > 0, and  1 = 1∕2. For 2 ≫ , we can then write (S15) For large , the signal source is only subject to the reflected signal and noise at cavity 1, leading to the approximate scattering matrix (S16) The suppression of mechanical noise works well for cavity 1, but cavity 2, where the signal emerges, is subject to amplified mechanical noise. As a consequence, the amplification is not quantum-limited, in addition to the lack of impedance matching.
Interchanging  1 ↔  2 leads to 1,out i.e., the amplification is quantum-limited, but the signal source is subject to amplified mechanical and optical noise.
We will now analyze the limits ≫ Γ m and ≪ Γ m , but note that they are only valid as long as  1 is smaller than the ratio of ∕Γ m and Γ m ∕ . In the limit ∕Γ m ≫ {1,  1 ,  2 }, we obtain which yields an amplitude bandwidth and intensity bandwidth Γ∕2. The gain-bandwidth product, defined as ≡ Γ √  evaluates to This implies that the gain-bandwidth product increases with gain. When  1 Γ m ∼ the approximation above breaks down.
In the opposite limit, ≪ Γ m , we instead have which implies that the bandwidth is close to ( 1 −  2 )∕ 1 and thus ≈ √ 4 2 ∕ 1 . For large gain, ≈ 2 , as before. The isolation bandwidth must be calculated separately. It is the range of frequencies over which sufficient isolation is attained. What sufficient means has to decided depending on the purpose. Close to = 0, to lowest order in , the reverse gain departs linearly from zero, namely † Thus, the isolation bandwidth is of order Γ m ∕ √ , independent of .

BANDWIDTH AND GAIN-BANDWIDTH PRODUCT OF DPSA
For the DPSA as discussed in the main text, we obtain Here, since both 2→1 ∝ √  2 and 1→2 ∝ √  2 , we can immediately conclude that the bandwidth is independent of the gain, and the gain-bandwidth product therefore unlimited. For and therefore such that the gain bandwidth Γ 1→2 = 2 1 Γ m .
To gain information about the departure from isolation, we expand the reverse gain around = 0. To first order (note that we do not assume { 1 ,  1 Γ m } ≪ here) Thus, the isolation bandwidth is again of order Γ m ∕ √  (but note that  takes different forms for DPSA and DPPA), independent of .

OFF-RESONANT TERMS
There are two kinds of off-resonant terms contained in Eq. (S4). One are "counterrotating terms", with a time-dependence (2Ω ), which are usually negligible. The other are off-resonant terms rotating at frequency Ω = Ω 2 + 2 − Ω 1 − 1 , which can have an appreciable effect [S4, S5]. For a related study see Ref. [S6]. The Hamiltonian describing the off-resonant terms is wherẽ 1 = 2 0, 1 ∕ 0, 2 (and the same for ↔ and 1 ↔ 2). Including off-resonant , the Hamiltonian is no longer timeindependent, but rather periodic, with period 2 ∕Ω. The resulting explicitly time-dependent Langevin equations can be mapped to stationary ones by use of a Floquet formalism [S7, S8], where we write system operators as Fourier series, for instance ( ) = ∑ exp( Ω ) ( ) ( ). We obtain Langevin equations without explicit time-dependence (diagonal in frequency space) where in Eq. (S33)(a), = 1 corresponds to the +-sign and = 2 to the −-sign. Note that in our conventions [ ( ) ( )] † = (− ) (− ). In principle, Eqs. (S33) constitute an infinitely large set of coupled linear equations (equivalently, an infinitely large matrix to invert). However, the coupling between different Fourier modes is suppressed by the mechanical and optical susceptibilities. In particular, the mechanical susceptibilities are strongly peaked (Γ m, ≪ Ω), such that it is a good approximation to let ( ≠0) = 0. Since in this approximation = (0) , we will omit the superscript (0) for in the following. Another consequence of the approximation is ( ) = 0 for | | > 1.
We can think of (±1) as four extra, "virtual" modes, as is done in Ref. [S4]. They are given by (1) whence we write down the new equation of motion for in Fourier space, e.g., The two types of terms that appear due to the off-resonant terms is one proportional to 1 that describes off-resonant cooling or heating, which can be incorporated into the susceptibility of the mechanical resonator, but also one that couples the first mechanical resonator to the second. The latter process only occurs when there is a drive on the red sideband of one resonator and one on the blue sideband of the other resonator. This is most easily understood when looking for example at the process underlying the term m,1 ( − Ω) † 2̃ * 11̃ 12 .̃ 12 is an interaction that creates a phonon in resonator 2 and a photon in cavity 1, but the process is off-resonant, meaning that the frequency of the photon created is approximately c,1 − Ω.̃ 11 shifts the frequencies the other way, such that this term mediates a resonant interaction between the mechanical resonators, with an off-resonant intermediate state. In contrast,̃ 11 would create a phonon in resonator 1 and a photon in cavity at frequency c,1 + Ω. Thus the term̃ 11̃ 12 would produce a phonon with frequency Ω 1 − 2Ω, a process that is strongly suppressed.
The spurious coupling between the resonators trivially vanishes in four-tone schemes, where one of or is zero for all , .
Eqs. (S38) do not have a contribution from the drives on the second cavity, because the blue and red drives are balanced, such that the dynamical backaction cancels. In the end, in an approximation where we neglect the frequency dependence of the self energies [i.e., c, ≈ ( ∕2 + Ω) −1 ], which is valid for small frequencies around resonance ≪ Ω, , the effect of the complex self-energies can be subsumed as a change of damping and detuning parameters.
We stress again that the other important conclusion from this analysis is that the off-resonant terms do not change the topology of the circuit, neither in the DPPA nor in the DPSA, such that the theory presented in the main text applies to a very good approximation. The embedding into a Floquet ansatz explains how the extra modes and their properties arise in the description of Ref. [S4]. The above calculation can be repeated for the proposed off-resonant cooling with a third cavity mode. This involves two drives detuned from the sidebands of the mechanical resonator (or one drive roughly in the middle of the two sidebands). For instance, consider pumps at frequencies c,3 − Ω + Δ , where c,3 is the frequency of the third cavity mode, Δ 1 = Δ, Δ 2 = −Δ, and Δ ≫ Γ m, . This leads to similar off-resonant terms as above, namely,

OFF-RESONANT COOLING WITH AUXILIARY MODE
This adds another contribution to the self-energy, similar to above The detuning choice ensures that phonons cannot hop from one resonator to the other via the auxiliary mode. The final result of this treatment is a new, enhanced damping rate Γ eff . However, the noise strength in the Langevin equations is unchanged. We can writė where the new effective noise has correlators ⟨̃ in ( )̃ in ( ′ )⟩ = thermal ( − ′ ) × (Γ m ∕Γ eff ) = thermal ( − ′ )∕Λ, where we have introduced a cooling parameter Λ = Γ eff ∕Γ m . In effect, we have modified parameters Γ eff = ΛΓ m and eff = thermal ∕Λ, parametrize by Λ, with the remainder of the Langevin equations unchanged. The effect of cooling with an additional drive is illustrated for the DPSA in Fig. S1. In this figure, we show the plaquette with the extra cooling mode, and plot the gain and noise in both directions for fixed coupling constants, but for varying levels of cooling Λ. The reason for keeping the coupling rates rather than the cooperativities unchanged is that high cooperativities become unattainable for strongly broadened mechanical resonators. We observe that the auxiliary cooling negatively affects gain, and forward added noise, which is due to the effective decrease in cooperativity. On the other hand, the isolation bandwidth is increased, since it depends on the mechanical linewidth. Most importantly though, the backward noise can be strongly reduced.

STABILITY Theory
It is important to know in which regimes the rotating-wave approximation (RWA) is valid. In order to numerically analyze the system beyond the RWA, we use the Hamiltonian Eq. (S4), but now keep all of the terms. In the resulting Hamiltonian, there are two frequencies present, Ω 2 + 1 and Ω 2 + 2 . Our approach is taken from Ref. [S7], but now extended to incorporate two frequencies. Note that other approaches are possible [S9]. In fact, we will choose the two frequencies to beΩ 1 = Ω 1 + 1 + Ω 2 + 2 andΩ 2 = Ω 1 + 1 − Ω 2 − 2 , since it will lead to a more compact Floquet matrix. If we restrict the theory to only the second Fourier frequencyΩ 2 we recover the theory from the two sections above. While this does capture changes in effective parameters, it does not show instabilities that are present in the full description. Collecting all system variables into a vector ⃗ = ( 1 , † 1 , 2 , † 2 , 1 , † 1 , 2 , † 2 ) , we can write the equations of motion aṡ where ( ) is sometimes called Langevin matrix and = diag( . Since the Hamiltonian has two Fourier frequencies, we can write the Langevin matrix in terms of its Fourier components as well For the Hamiltonian we consider, the non-zero Fourier components ( , ) are those with | |, | | ≤ 1. We can do the same to ⃗ , which defines ⃗ ( , ) ( ). However the Fourier components of ⃗ will still contain fluctuations, so they are not time-independent. We additionally perform a Fourier transform of all Fourier components of ⃗ , so that we can finally write the original equation of motion (S42) as (∀ , ) While in principle there are infinitely many coupled equations (one for each combination of , ∈ ℤ), we have to truncate at a certain number to make the problem tractable. From the truncated matrix we calculate the scattering matrix between the Fourier modes. Of particular interest here is the scattering between the zeroth Fourier modes, for we are interested in a signal on resonance. Forward and reverse gain are defined in the same way as in the main text, except that they now refer to the equivalent elements in relation to the scattering matrix between the zeroth Fourier modes. In general, we can distinguish four different regimes. For sideband parameters Ω ∕ , (Ω 1 − Ω 2 )∕ that are very large in comparison to the cooperativities  = ∕( Ω ), the RWA theory is fully valid. As the sideband parameters decrease, there is a regime where the effective parameters of the systems start to change. Since isolation relies on fine-tuning parameters, it is very sensitive to such a change of parameters. The RWA theory can be restored when working with renormalized parameters or when numerically optimizing the plaquette phase, as we demonstrate below. As the nonresonant terms become even stronger, i.e., through increasing cooperativities or a lower sideband parameter, the system has a qualitatively different response. Finally, there is a regime where the system invariably becomes unstable (cf. also Ref. [S9]).
In the whole section we take 1 = 2 = and Γ m,1 = Γ m,2 = Γ m = 10 −2 and choose units such that = 1 for simplicity. Our proposed additional sideband cooling can increase Γ m beyond that value. We note that a mechanical damping rate Γ m = quantitatively changes our conclusions below, since it enhances the detrimental effect of the counterrotating terms, but that the results are qualitatively the same.
Optimizing coupling rates at finite sideband parameters Here we show by example (cf. Fig. S2) that for moderate sideband parameters and cooperativities, optimizing the plaquette phase leads to near ideal behavior. In the middle panel we have numerically optimized the plaquette phase. As comparison, the right panel shows the ideal theory. Parameters are  1 = 4,  2 = 16, Ω 1 ∕ = 5, Ω 2 ∕ = 20. Although the differences may appear subtle, we can quantify them (all on resonance): For unoptimized parameters, reverse gain is 0.23, whereas forward gain is 61. After optimizing, we have reverse gain of the order of 2 × 10 −17 , which is essentially 0 within numerical errors, and forward gain 55. The ideal theory predicts 0 and 56.

Qualitative and quantitative deviations from ideal theory
For larger coupling rates we enter the regime where qualitative differences appear and the quantitative difference increase further. This is illustrated in Fig. S3. Optimizing the plaquette phase still recovers isolation, as before, but forward gain takes a qualitatively different form and is considerably reduced in comparison with ideal theory.

Instability threshold
In this subsection we demonstrate a method to find the instability threshold. We find that for a sideband parameter Ω∕ ≫ 1 the eigenvectors are centered around a specific Fourier frequency and only mix with adjacent frequencies, due to the fact that the coupling is highly off-resonant. In the converse case, the eigenvectors are spread out over many Fourier frequencies. Thus, in the resolved sideband regime, we have to include only few Fourier frequencies. Furthermore, due to the truncation at a certain Fourier index, there are eigenvectors that are concentrated at the edges of the Fourier domain that are not eigenvectors of the infinitely big Increasing the coupling strength further, we enter the regime where qualitative differences appear. The left panel shows forward and backward gain of the DPSA when using the coupling parameters of the ideal RWA theory. In the middle panel we have numerically optimized the plaquette phase. As comparison, the right panel shows the ideal theory. Parameters are  1 = 6,  2 = 36, Ω 1 ∕ = 5, Ω 2 ∕ = 20. While in Fig. S2 the differences were only quantitative, here the forward gain in the middle panel shows a strong deviation from ideal theory (right panel). Again we list forward and reverse gain on resonance: unoptimized parameters, reverse/forward gain of 12.7 and 105, after optimizing, reverse/forward gain are ∼ 8 × 10 −17 and 44.2. As a comparison, the ideal theory predicts 0 and 88. matrix. After solving the eigenvalue problem, we therefore select only those eigenvectors that have an appreciable support at the zeroth Fourier component. We then look for the instability in this restricted set of modes. We observe that there is an instability that occurs for certain cooperativities. For an analysis in a related system, see Ref. [S9].
We plot the instability threshold (in units of ) as a function of resonator frequency in Fig. S4. Note however that the plots in this section have been obtained for driving parameters given by the ideal theory. We can make some progress through optimizing the driving parameters in the presence of the off-resonant terms, as discussed in the previous subsections.  S4. We plot the cooperativity at which the instability occurs as a function of the two sideband parameters, Ω 1 ∕ and Ω 2 ∕ . The onset of stability is defined as the cooperativity  1 (note  2 =  2 1 ) at the real part of one of the eigenvalues becomes positive. The set of eigenvalues is restricted to those belonging to eigenvectors with support at the zeroth Fourier component.
In order to appropriately map the Langevin equations onto a master equation, we have to take into account that there are two coherent interactions and two baths. Therefore, we have to repeat this procedure twice, making the two parts of the coherent interaction