Stabilizing two-qubit entanglement by mimicking a squeezed environment

It is well known that qubits immersed in a squeezed vacuum environment exhibit many exotic phenomena, including dissipative entanglement stabilization. Here, we show that these effects only require interference between excitation and decay processes, and can be faithfully mimicked without non-classical light using simple classical temporal modulation. We present schemes that harnesses this idea to stabilize entanglement between two remote qubits coupled via a transmission line or waveguide, where either the qubit-waveguide coupling is modulated, or the qubits are directly driven. We analyze the resilience of these approaches against various imperfections, and also characterize the trade-off between the speed and quality of entanglement stabilization. Our protocols are compatible with state of the art cavity QED systems.

It is well known that qubits immersed in a squeezed vacuum environment exhibit many exotic phenomena, including dissipative entanglement stabilization. Here, we show that these effects only require interference between excitation and decay processes, and can be faithfully mimicked without non-classical light using simple classical temporal modulation. We present schemes that harnesses this idea to stabilize entanglement between two remote qubits coupled via a transmission line or waveguide, where either the qubit-waveguide coupling is modulated, or the qubits are directly driven. We analyze the resilience of these approaches against various imperfections, and also characterize the trade-off between the speed and quality of entanglement stabilization. Our protocols are compatible with state of the art cavity QED systems.

I. INTRODUCTION
Squeezed states of electromagnetic radiation [1,2], have long been of interest, with applications ranging from quantum metrology and sensing [3][4][5][6][7][8][9], to qubit readout [10][11][12][13][14] and enhanced quantum gates [15,16]. Driving systems with squeezed vacuum noise can also lead to interesting dissipative physics. A single qubit driven by broadband squeezed vacuum noise is described by the master equation [17] ρ q = γD [cosh(r)σ − + sinh(r)σ + ]ρ q , whereσ ± are standard qubit raising/lowering operators, D[x]ρ = xρx † − x † x,ρ /2 is the usual Lindblad dissipator, γ is the qubit-bath coupling rate, and r is the standard squeezing parameter. While the steady state of Eq. (1) is mixed, it generates non-trivial dynamics, with consequences including the existence of two distinct transverse relaxation times [17,18] and modifications of the Mollow triplet fluorescence spectrum [19][20][21][22]. The dissipative dynamics becomes even richer when two qubits are coupled to a broadband two-mode squeezed vacuum (TMSV) environment, as the dynamics can now prepare and stabilize remote qubit entanglement [23]. The evolution is described bẏ where γ 1/2 are the qubit-bath coupling rates, and r is the squeezing parameter of the TMSV. The photon number pairing correlations of the TMSV environment are inherited by the qubits, leading to a pure entangled steady state having the form G U a u y s n Q b 4 r u e n A y G v g v B v 7 7 5 7 3 j 1 x t B 7 q H H 6 A n q I x + 9 R M f o L R q j C a L o B / r V 2 m n t e p + 8 L 9 5 X 7 9 u a 2 m 5 t a h 6 h K + Z 9 / w 1 2 q E 1 r < / l a t e x i t > Squeezed vacuum Figure 1. Synthetic squeezing. a) The qubit dynamics produced by biharmonic modulation of the qubit-cavity coupling are equivalent to that for injected squeezed vacuum resonant with the cavity. b) Circuit QED schematic of twoqubit remote entanglement generation via synthetic squeezing. SQUID-modulated coupling between the qubits and a transmission line implements the dynamics of Eq. (2).
The dissipative physics of both Eqs. (1) and (2) are often taken to be hallmarks of the interaction of matter with non-classical radiation. Unfortunately, realizing these "dissipative squeezing" ideas experimentally is severely limited by the difficulty in generating, transporting and injecting squeezed states with high efficiency. In this work, we point out an underappreciated fact: the non-trivial dissipative dynamics of both Eqs. (1) and (2) can be realized without the use of any non-classical radiation, but instead by directly engineering interference between relevant qubit excitation and relaxation pathways. As we show in detail, these can be achieved in a variety of systems where the the qubit frequencies and/or couplings to a vacuum reservoir can be temporally modulated (depicted in Fig. 1a)). We term these schemes "synthetic squeezing", and show that they provide a powerful route arXiv:2110.06201v2 [quant-ph] 22 Jul 2023 towards entanglement stabilization.
Our work discusses several methods for implementing synthetic squeezing. This is of more than just academic interest: it also provides a potentially powerful new method for preparing and stabilizing entangled states between remote qubits using only classical drives. We discuss two versions of this idea: a first method that uses time-modulated couplings between qubits and a photonic link (i.e. a waveguide or transmission line structure), and a second method that only uses local Rabi drives on each qubit. Schemes similar to our first approach were studied in Refs. [? ? ] but the equivalence to driving with squeezed light was not discussed. The second approach both generalizes and provides intuition for a protocol first introduced in Ref. [24]. Including the impact of intrinsic thermal dissipation, we show that two-qubit entanglement stabilization with concurrence above 90% is achievable in contemporary circuit QED devices.
While relatively unexplored in the few qubit case, we note that the general idea of mimicking squeezed reservoirs with parametric or Raman processes has been studied in the context of stabilizing bosonic states. This includes schemes for producing bosonic squeezed states in trapped ions [25], optomechanics [26] and circuit QED [27], and bosonic optomechanical entanglement [28,29]. Related master equations have also been studied in the setting of spin ensembles [30,31], though the specific equivalence to driving with squeezed light was not discussed.
The remainder of this paper is organized as follows. In Secs. II and III we present our synthetic squeezed dissipation schemes for both the single and two qubit cases. In Sec. V, we study the approach to steady-state of the two-qubit scheme, its connections to large-spin entanglement stabilization, and the robustness of the steady-state entanglement to thermal dissipation. Finally, in Sec. VI we present our concluding remarks.

II. SINGLE QUBIT SYNTHETIC SQUEEZED DISSIPATION VIA COUPLING MODULATION
We begin by describing our scheme for synthetic dissipative squeezing of a single qubit, i.e. Eq. (1). The starting point is a cavity-qubit system with a time dependent coupling, described by the Hamiltonian: This is the usual Rabi Hamiltonian describing a cavity of frequency ω c coupled to a qubit of frequency Ω, with a time-dependent coupling g(t). We consider phase-locked double-sideband cosine modulation of the coupling at frequencies ω c ± Ω, described by b) Modulated coupling a) Squeezed reservoir (1). a) Interference in the qubit excitation and decay channels when subject to a squeezed vacuum environment. b) The blue and red sideband processes driven by coupling modulation at ωc±Ω interfere such that both qubit decay (red sideband) and excitation (blue sideband) result in an emitted cavity photon at frequency ωc.
After going to the interaction frame for the qubit and cavity, we make a rotating-wave approximation (RWA) to drop all time-dependent terms that oscillate faster than ∆ = ω c − Ω, and obtain the Hamiltonian where we require that |α ± | 2∆/ḡ for the RWA. Assuming that the cavity decays at a rate κ |α ± |ḡ into a zero temperature environment, adiabatic elimination of the cavity yields a qubit-only master equatioṅ With the associations α + = α sinh(r) and α − = α cosh(r) we obtain Eq. (1) with decay rate γ = 4ḡ 2 α 2 /κ and squeezing parameter r set by |α + /α − | = tanh(r). We thus have our first (and simplest) example of how classical time-modulation can be use to simulate the dissipative effects of a squeezed environment. We stress that this protocol does not require any input squeezed light, nor does it generate any bosonic squeezing in the cavity mode. Note that modulation of tuneable couplings to control qubit cooling and heating processes has been studied previously [34,40], but two-tone modulation and the connection to squeezing was not discussed. It is worth emphasizing the heuristic reasons why interfering sidebands mimic the dissipative effects of squeezed vacuum. For a squeezed vacuum reservoir, the interference of qubit excitation and relaxation processes encoded in Eq.  . Two-qubit interference. a) and b) interference between the qubit decay and excitation processes mediated by the coupling modulation that results in the combined dis-sipatorsĴ a/b . c) For a cavity with a broad linewidth, the two dissipative qubit channels sample the cavity density of states at different frequencies, and effectively interact with distinct environments. d) In the transmission line setup, due to the modulated coupling both qubit 1 (2) decay and qubit 2 (1) excitation result in a TL excitation at Ω1 (Ω2).
2M + 1, this could have resulted from a qubit excitation process starting with 2M + 2 reservoir photons, or from a qubit relaxation process starting from 2M reservoir photons (see Fig. 2a)). In our synthetic scheme, an identical interference is achieved using only classical radiation, i.e. the coupling modulation. The interference is now between the blue sideband Raman process (driven by coupling modulation at frequency ω c + Ω) with the red sideband Raman process (driven by the modulation at frequency ω c −Ω). Both processes result in the generation of a cavity photon, as depicted graphically in Fig. 2b).

III. TWO-QUBIT SYNTHETIC SQUEEZED DISSIPATION VIA COUPLING MODULATION
The above ideas can easily be extended to two qubits (frequencies Ω 1 and Ω 2 ), allowing one to mimic the effects of a two-mode squeezed reservoir; we refer to this general idea as synthetic two-mode squeezing. We now use phase-locked coupling modulation between the qubits and two independent decay channels (mediated by bosonic modes). As in the single qubit case, the key ingredient is interference between red/blue sideband Raman processes that result in photon emission at the same frequency. However, in the two-qubit setup the red sideband for one qubit and the blue sideband for the other interfere, as depicted in Fig. 3a) and b). This results in the reduced qubit-only master equatioṅ describing the independent qubit dissipative channels mediated by the bosonic modes shown in Fig. 3a) and 3b). Hereḡ and α ± (β ± ) define the qubit-mode coupling modulation for mode a (b), in direct analogy to Eq. (5). By setting α − = β − = α cosh(r) and α + = β + = α sinh(r) we obtain Eq. (2) with γ 1 = γ 2 = 4ḡ 2 α 2 /κ. We now present three setups to engineer Eq. (8), leaving detail of the calculations to App. A, .
The first and simplest setup consists of two cavities, each coupled to both qubits, as described by the coupling Hamiltonian The coupling g 1 (t) is modulated such that the a cavity mediates the interference process of Fig. 3a) described by the dissipatorĴ a , and similarly g 2 (t) and the b cavity mediate Fig. 3b) andĴ b . The a and b modes need not be spatially localized: they could, e.g. , be the hybridized supermodes of two tunnel-coupled physical cavities, with the underlying localized cavity modes coupling to only one qubit each. If instead we interfere both red-sideband processes and both blue-sideband processes, as can be done with trivial modification of this setup, we would obtain the dissipatorŝ which leads to Eq. (2) with the state of qubit 2 flipped. The second setup uses a single cavity with a broad linewidth (κ) to supply the two distinct dissipative channels for the qubits, as depicted in Fig. 3c), with the cavity coupling to both qubits with the same modulated coupling, given by the Hamiltonian In this setup, as a result of the modulated coupling, one pair of sideband processes results in cavity photon emission at frequency ω a + δΩ, and the other pair results in emission at ω a − δΩ, with δΩ = (Ω 1 − Ω 2 )/2. These two mediated decay processes are independent and add incoherently, which leads to Eq. (8). Note that this setup cannot easily be modified to directly engineer Eq. (11). Finally, the third setup is two qubits connected to a long waveguide or transmission line (TL) with static and modulated coupling, as depicted in Figs. 1b) and 3d). Modulation of qubit 1 coupling at the qubit sumfrequency drives a Raman process where qubit 1 is excited and a photon of frequency Ω 2 is emitted into the TL. Due to the static coupling, qubit 2 also decays into the TL at frequency Ω 2 . These two processes interfere, as do the analogous processes at emission frequency Ω 1 when the qubits' roles are reversed, which results in qubit evolution described by Eq. (8). While this setup is sensitive to the distance between the qubits (as we do not want a waveguide-mediated Hamiltonian qubit-qubit interaction, see App. A 3), it can be used to generate and stabilize steady-state entanglement between two remote qubits. We note that several recent circuit QED experiments studied setups with tuneable qubit-transmission line couplings [35,36], a platform that could be ideal for the above implementation.

IV. TWO-QUBIT SYNTHETIC SQUEEZED DISSIPATION VIA LOCAL LINEAR DRIVING AND COLLECTIVE LOSS
The idea of generating remote entanglement by synthetic two-mode squeezing is appealing; however, not all architectures facilitate dynamically modulated qubitwaveguide couplings. In this section, we show that one can still use synthetic squeezing ideas to stabilize entanglement by simply combining the passive, collective loss generated by a photonic link with tailored local qubit Rabi drives. In the correct reference frame, the physics of this setup is almost directly analogous to driving two qubits with two-mode squeezed light (c.f. Eq. (2)). Our discussion here both generalizes and helps provide intuition for the dissipative entanglement protocol presented in Ref. [24].
Consider two qubits which are coupled to a common source of loss at a rate Γ. The dissipation is described by the master equationρ q = ΓD Ĵ ρ q with the collective loss operatorĴ where η allows for an asymmetry between the qubit-bath couplings. For clarity we take η = 1 in what follows; as we show later, analogous results hold for η = 1.Ĵ has a two-parameter family of dark states of the form where 0 ≤ α ≤ 1. These states are entangled for any α > 0, whereas for α = 0 we have the trivial dark state |00 . Any state (pure or impure) in the span of this dark state manifold is a possible steady state, making this dissipative evolution of little use if the goal is entanglement stabilization. As such, we would like to introduce additional simple dynamics to the qubits so that there is a unique entangled dark state (for some α > 0). We will solve this problem by making a direct mapping to entanglement via two-mode squeezed dissipation. This is possible because the dark states in Eq. (14) are unitarily equivalent to the "paired" family of entangled states associated with two-mode squeezed dissipation, i.e. states of the form in Eq. (3). Without loss of generality we set φ = 0; in this case we have the explicit mapping It corresponds to opposite local rotations of each qubit about the y-axisÛ where qubit 1 (2) takes the upper (lower) sign. The rotation angle θ is a function of α (or equivalently r), and is given by Picking out a non-trivial dark state is now mapped in the new frame to ensuring the steady state has a particular non-zero value of the squeezing parameter r. This latter task is something we know how to do, namely by using the synthetic squeezing master equation in Eq. (2).
In what follows, we choose a particular value of r = r 0 , corresponding to the entangled state we would like to stabilize; we will also viewÛ and |Ψ q to be functions of r rather than α. The next step in our mapping to synthetic two-mode squeezing is to examine the form of the collective loss dissipatorĴ in the new frame defined byÛ [r 0 ]. We havê Here,Ĵ 1 andĴ 2 are precisely the two squeezing dissipators needed for synthetic squeezing (c.f. Eq. 2)) Further,Ĵ Z is a correlated dephasing term Hence, in the new frame defined by our chosen value of r = r 0 , the original collective loss dissipator has a direct connection to the synthetic squeezing dissipators in Eq. (2): J is the sum of these dissipators, plus a dephasing dissipator. One can easily check that the paired target state |Ψ[r 0 ] q is dark with respect to this sum dissipator. However, it is not unique, which is to be expected: all we have done at this stage is re-written our original problem in a new frame, hence the transformed version of any dark state from Eq. (14) (including |00 ) remains a dark state.
However, in our new frame defined byÛ [r 0 ], the state |Ψ[r 0 ] q is the only dark state of the dissipatorĴ that does not have any single-excitation components, i.e. kets |01 or |10 . This suggests a very simple strategy to make the paired state |Ψ[r 0 ] q a unique dark state. By adding a Hamiltonian that (in the new frame) breaks the degeneracy between the paired and unpaired subspaces, such asĤ we introduce an energy gap µ between the paired states |00 and |11 (which have zero energy) and the unpaired states |10 and |01 . The above extra Hamiltonian dynamics is exactly what we need to pick out a unique entangled steady state. One can directly confirm that for any µ > 0, the target entangled state |Ψ[r 0 ] q is the unique steady state of the dynamics generated by our master equation, which in the new frame has the forṁ Here,ρ =Û [r 0 ]ρÛ −1 [r 0 ] is the transformed-frame density matrix of the two qubits. Note that for large µ Γ, the connection to the synthetic squeezing master equation of Eq. (2) is even more direct. In this limit,Ĥ will causeĴ 1/2 to oscillate at frequencies ±µ, whileĴ Z remains static. In this case, the dynamics is well approximated by dropping cross-terms in the dissipator (as they will be fast oscillating), leading to the approximatioṅ (24) Apart from the last collective dephasing dissipator, the dissipative dynamics is now completely equivalent to the synthetic two-mode squeezing master equation in Eq. (2). Importantly, |Ψ[r 0 ] q is still the unique steady state of this master equation. Moreover, as we discuss in App. B 2, this additional dephasing has little impact on the dynamics in the limit µ Γ. In this limit we can also remove the Hamiltonian term by moving to the interaction frame,ρ = exp(−iĤ t)ρ exp(iĤ t), where the equivalence to Eq. (2) is more obvious, and |Ψ[r 0 ] q still remains the unique steady state.
Finally, having realized a master equation that is strongly analogous to synthetic two-mode squeezing in the new frame, we can move back to the original frame (i.e. via the unitaryÛ † [r 0 ]). The master equation in the original frame has the forṁ whereĴ is the simple collective loss dissipator in Eq. (13), and the lab-frame HamiltonianĤ =Û † [r 0 ]Ĥ Û [r 0 ] has the form with This Hamiltonian describes a two-qubit system subject to Rabi drives at the same frequency ν and amplitude Λ, in the rotating frame determined by the drive frequency. Thus, ±∆ describe the detuning of the drive from each qubit, such that ν is the average of the two qubit frequencies. To be explicit, in the lab frame where the drives are time-dependent, the Hamiltonian has the form We note that the master equation in Eq. (25) was first presented in Ref. [24]. Our derivation above adds to this work by showing the connection to squeezed dissipation, and by providing a clear intuitive explanation for why the scheme stabilizes entanglement. Further, our approach also makes it easy to generalize to the case where the qubits do not couple equally to the collective loss channel (i.e. η = 1 in Eq. (13)), something not previously considered. The asymmetric case is presented in detail in App. B, and the key finding is that dissipative entanglement stabilization is still possible, but the parameters in the required driving Hamiltonian are different: and where ∆, Λ are no longer given by Eq. (27) but are determined via Eqs. (B13) and (B14) in App. B. While the above scheme is attractive as it does not require the ability to modulate the qubit-waveguide coupling, it does have other limitations when compared to the modulated-coupling approach. For the symmetric η = 1 case and a large squeezing parameter r, the Rabidrive amplitude scales as Λ ∼ e 2r (Ω 1 − Ω 1 )/2. Hence, for qubits that are far detuned from one another, one requires extremely large Rabi drive amplitudes. In contrast, the modulated-coupled approach of Sec. III has no analogous limitation.
Finally, we briefly note some of the imperfections of the above scheme. The collective loss dissipator is typically realized by coupling the qubits to a waveguide. When the separation of the qubits along the waveguide is correctly set, only collective loss dissipation is generated as desired. However, if the qubit spacing is not equal to one of the special values, the waveguide will also generate single qubit loss dissipators, as well as a coherent tunneling termĤ wg = J(σ (1) +σ (2) − + h.c.). As discussed in Ref. [24], single qubit loss rapidly degrades the entanglement of the steady state. The coherent tunneling term will also degrade the steady state entanglement, however, one can show numerically the steady state entanglement is much more sensitive to single qubit loss than to the tunneling term. Intuitively, this is because one of the eigenstates ofĤ wg is the α = 1 dark state of the family of Eq. (14) (equivalently the r → ∞ two-qubit squeezed state). The modulated-coupling scheme involving a waveguide discussed in detail in App. A 3 shows an analogous result.

V. TWO-QUBIT ENTANGLEMENT STABILIZATION
As we have discussed, the dynamics of Eq. (2) can be engineered without squeezed vacuum, either using our proposed coupling modulation schemes (see App. A for further details), or the linear driving and correlated loss scheme of Ref. [24] that we have generalized in the previous section (details in App. B). Thus, it is useful to understand the limits of using Eq. (2) for high-quality entangled state stabilization, as it is the generic description of many schemes.
The steady-state of this evolution is the entangled state of Eq. (3), which for r → ∞ tends to a maximally entangled Bell state. Our implementation of Eq. (8) does not require increasing modulation amplitudes to increase r, which would otherwise be an issue given the RWA used to drop fast-oscillating time-dependent terms from Eq. (8). For increasing r, α can be decreased to keep α ± and β± constant. As in the single qubit case, this comes at the cost of reducing the overall decay rate Γ ≡ 4ḡ 2 α 2 /κ. Keeping the modulation amplitudes fixed, more entanglement requires a longer preparation time, a trade-off also observed in similar bosonic schemes [29].
Moreover, saturating the r → ∞ limit leads to a degeneracy in the steady-state of Eq. (2) (see App. C for further details). This degeneracy can be most easily understood by considering dynamics in the r → ∞ limit of Eq. (11) which is unitarily equivalent to Eq. (2). In this form, it is clear that the steady-state degeneracy is due to the evolution conserving total angular momentum of the two qubits (implying the existence of distinct singlet and triplet steady states). In light of this, for r large but not infinite, we expect a slow approach to steady-state due to the weak breaking of total angular momentum conservation. This slow dynamical rate is intrinsic to Eq. (2), as opposed to the slow decay rate discussed in the previous paragraph, which is an effect of our engineering of Eq. (8).

A. Slow Dissipation: Liouvillian Gap and Large-Spin Mean-Field
We now focus on the intrinsically slow dynamics of Eq. (2), which correspond to decay out of the nearly degenerate subspace discussed in App. C. Writingρ q = L(ρ q ), where L corresponds to Eq. (8), we can vectorize this equation as˙ ρ q =L ρ q , where the super-operator ma-trixL describes the dissipative evolution of the qubits. The eigenvalues of this matrix describe the dynamical timescales of the system evolution, and we find that there is one slow rate, corresponding to the eigenvalue with the smallest real part, often referred to as the Liouvillian gap. In the large r limit it scales inversely with r as .
Note that all other eigenvalues monotonically increase with r, scaling as sinh 2 (r) to leading order. This intrinsically slow rate implies that there is a tradeoff between the stabilization speed-limit and the amount of entanglement in steady-state (cf. Eq. (3)). The impact on a specific protocol depends on whether "typical" initial conditions evolve to the steady-state via the process described by this slow rate. We leave this question to future work, but note that prior work has shown that Eq. (2) is quite robust to additional sources of qubit decoherence [23], which ultimately compete with this slow rate and determine the asymptotic entanglement.
This slow rate is a property of the two qubit setup, and does not exist for its bosonic analogue [29], whereσ − 1/2 are replaced by bosonic modesd 1/2 , as the bosonic model has a unique steady-state for balanced dissipators. To understand the transition between these two extremes, we consider a system of two spin-S particles, where S 1/2. The bosonic model is often used as an approximation to such a large spin model via the Holstein-Primakoff (HP) transformation [41], which connects the spin operatorŜ + to a bosonic mode viaŜ + =d † 2S −d †d . In the HP representation, the large spin version of Eq. (8) has dissipatorŝ withn k =d † kd k , and an enhanced dissipative rate Γ S ≡ 16Sḡ 2 α 2 /κ. In the limit S n k the square root factors can be neglected and the bosonic model is obtained.
However, the steady-state of the bosonic model has n k = sinh 2 (r), and the square root factor of Eqs. (33) and (34) cannot be ignored. Treating this factor in the mean-field limit, we obtain a scalar that can be absorbed into an effective asymptotic dissipative rate for the large spin model This rate slows as r increases, similar to what was observed in the two-qubit (S = 1/2) model. Dynamically, Γ MF S will decrease as n k increases during the protocol. This mean-field result is a strong indication that any finite spin model may have a slow rate that results in a trade-off between the speed and quality of entanglement stabilization. Further study, beyond mean-field, is warranted to confirm this result, but is outside the scope of this work.

B. Thermal Dissipative Squeezing
The effect of thermal occupation of the bosonic environment of Eq. (8) is especially relevant in the transmission line version of our scheme, and remains unexplored in prior literature. To study this, we consider the thermal dissipative squeezing master equatioṅ withĴ a andĴ b as before for α − = β − = α cosh(r), α + = β + = α sinh(r), and wheren th is the average photon population of the modes. We characterize the effect of thermal occupation by measuring the concurrence of the two-qubit steady-state as a function of transmission line temperature, with the thermal occupation as a function of temperature given by the Bose-Einstein distribution. For qubits with a frequency around 6 GHz (as is typical in circuit QED) this is shown in Fig. 4. As can be seen, the two-qubit steadystate remains highly entangled (> 0.9 concurrence) up to a temperature of around 70 mK. For comparison, we also plot the purity of the steady-state as a function of temperature. While 70 mK is well above the operating temperature of many cryogenic circuit QED experiments, the relevant quantity is the thermal occupation of the TL near the qubit frequency, which includes thermal noise from classical control lines. If this control line noise is controlled via appropriate filtering [42], high-quality remote entanglement can be generated by our TL scheme.

VI. CONCLUSION
In this paper, we have presented an overlooked connection linking the dynamics of a system of qubits immersed in a squeezed environment, and reservoir engineering using interference of qubit decay and excitation pathways. We have shown that both situations result in the exact same dissipative dynamics, and have proposed schemes which achieve this interference using parametric coupling modulation of the Jaynes-Cummings interaction of cavity-and waveguide-QED. The latter can be used to remotely entangle qubits coupled to a transmission line.
Such modulated coupling is readily available in circuit-QED [32][33][34][35][36], and our scheme could immediately be implemented in such architectures. An alternative implementation of the modulated coupling of our schemes uses a dispersive interaction with classical cavity displacement [14,39]. In either case, the bare-coupling can be very weak, as the amplitude of the modulation can be used to compensate. Adapting to the specifics of a particular setup, and understanding the impact of spurious decoherence sources will be the focus of future work.
The results of our work overturn the long held belief that the unique effects of qubit-squeezed vacuum interaction were a feature of the nonclassical nature of the squeezed environment. In effect, we have traded the nonclassical radiation for a modulated nonlinear coupling. This combination of nonlinearity and classical modulation can be compared to the operation of a parametric amplifier, whose nonlinearity is provided by one or more two-level systems (qubits) coupled to a bosonic mode, such as in a Josephson parametric amplifier [43,44]. In the parametric amplifier it is the bosonic field that is squeezed, while in our scheme, the bosonic mode is never squeezed, but rather the qubit undergoes the dissipative dynamics corresponding to a squeezed environment.
We note that schemes for dissipative Bell state stabilization involving two qubits, one or more coupled cavities, and classical qubit-cavity control have been proposed and implemented in circuit-QED [24,38,[45][46][47][48][49][50]. Our work is the only scheme to make explicit connection to dissipative stabilization using squeezed vacuum, and it remains to be seen if this connection can be made for other proposals that use classical drives or coupling modulation. We conjecture that this connection holds approximately for all schemes with an intrinsically slow rate, as we have shown for the scheme of Ref. [24]. For schemes with no slow rate, such as Ref. [50], understanding where the connection to dissipative squeezing breaks is a topic worthy of future study.

ACKNOWLEDGMENTS
This material is based upon work for which LG was supported by the U.S. Department of Energy, Office of Science, under Award Number DE-SC0019461, and for which AL and AC were supported by the National Science Foundation QLCI for HQAN (NSF award 2016136).
Appendix A: Two Qubit Master Equation Derivations

Two Qubits and Two Cavities
An instructive example is to directly extend the single qubit setup to two qubits, by considering two qubits and cavities coupled via the Hamiltonian with the modulated couplings Both qubits must couple to both cavities, to allow for the interference required by dissipative squeezing. However, the symmetry in the qubit-cavity couplings shown above is not required, and can be compensated for by the modulation amplitudes. Besides the obvious physical implementation of two qubits both coupled to each physical cavity, this setup could also be implemented in a scheme where each qubit couples to only one physical cavity, and the cavities are tunnel coupled. The cavity supermodes would then play the role of the modesâ and b in Eq. (A1). Following the same procedure as for one qubit, we move to the interaction frame for Eq. (A1) and drop all fast-oscillating terms to obtain H =ḡâ † α +σ where we have assumed that |Ω 1 − Ω 2 | ḡ so that all time-dependent terms are fast oscillating and can be neglected. We can eliminate the cavities to obtain an effective master equation for the qubits when κ a , κ b |α ± |ḡ, |β ± |ḡ. In the case where κ a = κ b = κ, this giveṡ as depicted in Fig. 3a) and b). Starting with equal couplingḡ, and setting κ a = κ b = κ is a completely general choice, as differences in the couplings or decay rates can be compensated for by control of the modulation amplitudes.
We note that if instead we modulate with red-sideband tones only on modeâ, and blue sideband tones on modê b, i.e. g 1 (t) =ḡ α − e −i(ωa−Ω1)t + α + e −i(ωa−Ω2)t + h.c., then by a similar procedure we would obtain the dissipators of Eq. (11). Unfortunately, the other two approaches to engineering Eq. (8) discussed in this appendix do not have a similar simple modification that would enable engineering of Eq. (11).

Two Qubits and One Cavity
For a cavity with sufficiently broad linewidth it is also possible to synthetically generate two-qubit dissipative squeezing using our scheme with only one cavity. We start with the Hamiltonian with the modulated coupling whereΩ = (Ω 1 + Ω 2 )/2 is the average qubit frequency. We then move to the interaction frame to obtain the Hamiltonian H =ḡâ † e iδΩt α +σ (1) where δΩ = (Ω 1 − Ω 2 )/2 and Σ 1/2 = (3Ω 1/2 + Ω 2/1 )/2. Ifḡ > δΩ we cannot ignore terms that oscillate at ±δΩ, as they do not average away over timescales of interest. All other terms are manifestly fast oscillating and will average away, except for the terms in the last line at frequencies 2ω a − Σ 1/2 . For these to be fast oscillating, we require that 2ω a − Σ 1/2 ḡ. This simply implies that when designing the system, we must be careful to avoid the "resonance" conditions 2ω a = Σ 1/2 , so as not to activate unwanted qubit evolution.
Having dropped all fast oscillating terms the Hamiltonian bcomeŝ H =ḡâ † e iδΩt α +σ (1) Note thatḡ > δΩ is the opposite operating condition to the two cavity setup, where we require δΩ >ḡ.
Assuming the cavity has a sufficiently broad linewidth, i.e. δΩ κ, then the qubit processes in Eq. (A8) will interact with the environment through the cavity's nonzero density of states at ±δΩ, as depicted in Fig. 3c). If 2δΩ is greater than the effective qubit linewidths, which are at least 4α 2ḡ2 /κ, then the processes at +δΩ and −δΩ will effectively see independent environments. In this case, adiabatically eliminating the cavity would result in a two-qubit master equation of the form of Eq. (8).

Two qubits and a transmission line
We consider two qubits coupled to a long waveguide or transmission line (TL) described by the Hamiltonian We consider the general case where the TL may be an engineered metamaterial, and hence is described by a general band structure describing propagating photonic modes (wavevector k, band index λ) where we assume the dispersion relations satisfy ω k,λ = ω −k,λ . The qubit-TL couplings are described bŷ with x α denoting the position of the α qubit along the TL. f kλ is a wavevector and band dependent coupling coefficient which we assume to satisfy f kλ = f * −kλ . To realize the synthetic squeezing dissipator, the qubit-TL couplings g α (t) are modulated at the sum of the qubit splitting frequencies, and also have a static component: Note that in this description g has units Hz √ m.
After moving into the interaction picture with respect to the system HamiltonianĤ S = 2 α=1 (Ω α /2)σ (α) z and the TL HamiltonianĤ T L , the system-TL interaction can be written aŝ which is of the form with α ∈ {1, 2} indexing the qubits and ω ∈ {±Ω 1 , ±Ω 2 , ±(2Ω 1 + Ω 2 ), ±(Ω 1 + 2Ω 2 )} indexing all of the relevant frequencies. TheÂ α (ω) are e.g.,Â 1 (Ω 1 ) = gη 1σ (1) − andÂ 2 (Ω 1 ) = g 1σ (2) + . WithĤ I,int written in this form and following the master equation derivation of [51], we find the system density matrix master equation to be where we have introduced the bath correlation functions . We next decompose the bath correlation functions as for S αβ (ω) = (Γ αβ (ω) − Γ * βα (ω))/2i and γ αβ (ω) = Γ αβ (ω)+Γ * βα (ω). This decomposition separates the master equation into dissipative terms with correlation function coefficients γ αβ (ω) and Hamiltonian terms with correlation function coefficients S αβ (ω). The master equation can then be recast into the form where the coherent evolution of the system due to the transmission line is given by the "Lamb shift" Hamilto-nianĤ and the dissipation is given by The Lamb shift Hamiltonian contains two kinds of terms: local terms that correspond to energy shifts of each qubit, and more interestingly, a coherent qubit-qubit "pairing" interaction of the form The local qubit energy shifts simply renormalize the bare qubit frequencies, and hence do not play any significant role. In contrast, the pairing interaction can in principle disrupt our entanglement stabilization scheme; we discuss this more later in this section. Before doing this, we first analyze the form of the TL induced dissipation on the qubits. The dissipative terms include dissipative processes whose rates are governed by the TL correlation function γ αβ (ω) for all possible frequencies ω ∈ {±Ω 1 , ±Ω 2 , ±(2Ω 1 + Ω 2 ), ±(Ω 1 + 2Ω 2 )}. We let the TL be at zero temperature so that Thus all negative frequency processes (which correspond to emission of photons from the TL) vanish. While this is necessary to obtain the exact synthetic squeezing dissipation we desire, in Sec. V B we investigate the impact of a thermal environment on entanglement stabilization. The dissipative processes at frequency 2Ω 1 + Ω 2 (Ω 1 + 2Ω 2 ) involve both single qubit loss on qubit 1 (2), due to the static couplings η i , as well as single qubit excitation on qubit 2 (1), due to the modulated couplings i . These dissipative processes destroy the entanglement we seek to generate. Therefore we must engineer the TL to have a vanishing density of states at these frequencies. This could be arranged through standard filtering or bandgap engineering techniques. When the TL has zero density of states at these frequencies, the correlation functions become γ αβ (2Ω 1 + Ω 2 ) = γ αβ (Ω 1 + 2Ω 2 ) = 0. (A24) We are thus left with the dissipative processes at Ω 1 and Ω 2 . The Ω 1 processes take the form where D[Ô]ρ is the usual Lindblad dissipator and D[Ô,P ]ρ =ÔρP † − {P †Ô ,ρ}/2. The Ω 2 process is the same but with the indices 1 ↔ 2. Note that if the correlation functions γ αβ (Ω 1 ) = γ 1 were all equal, this dissipator would reduce to the collective dissipator D Ω1 (ρ(t)) = g 2 γ 1 D[η 1σ (1) The above correlation functions are readily evaluated in general (in the continuum limit 1 L k → dk 2π ) as where the {k j } are all wavevectors for which ω λ (k) = ω and v g,λ (k) = |∂ω λ (k)/∂k| is the group velocity at frequency ω λ (k). Then for the Ω 1 processes, we have where we have dropped the band index λ and absorbed the coupling coefficient f λ (k) into the group velocity.
Here k 1 is the wavevector for which ω λ (k 1 ) = Ω 1 for the correct band λ. We get a similar result for the Ω 2 processes.
Defining the qubit spacing we can rewrite the Ω 1 dissipation as In the ideal case where k 1 l = πn for n ∈ Z, this all reduces to a single dissipator with collective jump operator where the sign is set by the parity of n (+ for even n).
In the more general case where k 1 l = πn, Eq. (A30) tells us the dissipation will both have the desired collective dissipator, but also unwanted uncorrelated single qubit dissipators (loss on qubit 1, excitation on qubit 2). For a small spacing error k 1 l = πn + δ where δ 1, these additional dissipators are given to leading order in δ by As we discuss further below, our scheme can still perform well if these terms are non-zero but have a small amplitude (i.e. the error in the qubit spacing is small). A similar analysis holds for the dissipative processes involving a frequency Ω 2 in the TL. The dissipation reduces to a single collective jump operator when k 2 l = πm for m ∈ Z, with the sign set by the parity of m. Finally, if the qubit frequencies are chosen such that the characteristic wavevectors k 1 , k 2 satisfy k 1 l = πn and where usually |δ j | 10 −3 . Thus Λ may not exactly vanish at the synthetic squeezing conditions, but its value is typically quite suppressed.
To get some sense for how this scheme performs in the face of qubit spacing imperfections, l = l ideal +δl, we consider a simple model system and analyze its performance numerically. To be concrete, consider two superconducting qubits with frequencies Ω 1 = 4 GHz and Ω 2 = 6 GHz and a transmission line with a linear dispersion over the range 0 to 10 GHz, and with a (sharp) cutoff at 10 GHz (we require the transmission line to have zero density of states in a stop-band spanning at least 2Ω 1 + Ω 2 = 14 GHz and Ω 1 + 2Ω 2 = 16 GHz). We choose to space the qubits by a distance or approximately two wavelengths at Ω 1 (three wavelengths at Ω 2 ). The ideal spacing, δl = 0, satisfies the synthetic squeezing conditions k 1 l = 4π and k 2 l = 6π. We letḡ 2 /v g ≡ 1 and set η 1 = η 2 = cosh(r) and 1 = 2 = sinh(r). Then to a very good approximation Λ = sinh(r) cosh(r) (sin(k 1 δl) + sin(k 2 δl)) .
Thus the master equation is given bẏ . We numerically compute the steady state of this master equation as a function of δl and r. In general for δl = 0 the concurrence of the steady state is less than 1 but varies as a function of r. Thus, we numerically maximize the concurrence for a given δl by optimizing over r. The results of this optimization are shown in Fig. 5. The concurrence, state purity, and optimal squeezing parameter are plotted against the spacing error |δl| as a fraction of the wavelength λ 1 = 2π/k 1 . In the plot, we compare the full master equation (solid curves) to the master equation wherein the Hamiltonian is excluded (dashed curves) to show that the single qubit dissipators cause the bulk of the performance degradation. bitrary η > 0. In the main text we sought to provide the intuition for why this scheme works by exploiting the unitary equivalence between the family of the collective loss dissipator dark states Eq. (14) and the family of target entangled states Eq. (15). Even when the loss dissipator is asymmetric, the local unitary equivalence still holds, but the local qubit rotations are no longer of equal angle but opposite sign. As finding the angles when η = 1 is nontrivial, we instead take a different approach. Given that the steady state we seek is simultaneously a dark state of the collective loss operator and a zero-energy eigenstate of the Hamiltonian, we can simply find the Hamiltonian parameters for which a given dark state of Eq. (13) is an eigenstate. We found that for the symmetric case the Hamiltonian is simply local Rabi drives on the qubits given by Eq. (26). We allow for a more general Hamiltonian wherein we allow for asymmetric Rabi drive amplitudes Λ (j) and asymmetric detunings parameterized by . The dark states of the collective of loss Eq. (13) can be expressed as for |α| 2 + (1 + η 2 )|β| 2 = 1. We take α, β > 0 real and positive without loss of generality. Enforcing that |φ (2) (with rates γ1 = γ2 = Γ). For small µ Γ the dissipative gap vanishes as κ slow ∼ µ 2 , and for large µ Γ the gap saturates at a constant which is smaller than the ideal master equation gap by a factor of e −2r . This further exponential suppression of the gap is due to the overall prefactors of e −r in the squeezing dissipatorsĴ 1/2 (c.f. Eqs. (19) and (20)).
parameter r, which is the limiting factor for relaxation to steady-state. The gap is thus the relevant metric for how closely Eq. (24) approximates the dynamics of Eq. (2). To compare, we numerically calculate the dissipative gap of Eq. (B19) as a function of µ/Γ and compare it to the gap of Eq. (2). The results of this are shown in Fig. 6, where the solid curves show the gap for various r and the dashed lines show the gap of Eq. (2). Recall that the gap closes at µ = 0 (no Hamiltonian) because the collective loss dissipator has multiple steady states.
A simple perturbative argument suggests that the gap should open as κ slow ∼ µ 2 which is precisely what is observed for µ Γ. When µ Γ the gap saturates at a value e −2r times smaller than the ideal master equation gap. This factor comes from the prefactors of the jump operatorsĴ 1/2 given in Eqs. (19) and (20). This extra suppression of the dissipative gap, and hence the state stabilization rate, is another important consideration for the practical implementation of this scheme, in addition to those discussed at the end of Sec. IV.