# Stabilizing two-qubit entanglement with engineered synthetic squeezing

L. C. G. Govia<EMAIL_ADDRESS>Quantum Engineering and Computing, Raytheon
BBN Technologies, 10 Moulton St., Cambridge, MA, USA.  Pritzker School of
Molecular Engineering, University of Chicago, Chicago, IL, USA.  A.
Lingenfelter  Pritzker School of Molecular Engineering, University of Chicago,
Chicago, IL, USA.  A. A. Clerk  Pritzker School of Molecular Engineering,
University of Chicago, Chicago, IL, USA.

###### Abstract

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 Gardiner and Zoller (2004); Walls
(1983), have long been of interest, with applications ranging from quantum
metrology and sensing Giovannetti  _et al._ (2004); Caves  _et al._ (1980);
Braginsky and Khalili (1992); LIGO Scientific Collaboration (2013); Iwasawa
_et al._ (2013); Peano  _et al._ (2015); Taylor  _et al._ (2013), to qubit
readout Barzanjeh  _et al._ (2014); Didier  _et al._ (2015a, b); Govia and
Clerk (2017); Eddins  _et al._ (2018) and enhanced quantum gates Puri and
Blais (2016); Royer  _et al._ (2017). 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
Gardiner (1986)

$\displaystyle\dot{\hat{\rho}}_{q}=\gamma\mathcal{D}\left[\cosh(r)\hat{\sigma}_{-}+\sinh(r)\hat{\sigma}_{+}\right]\hat{\rho}_{q},$
(1)

where $\hat{\sigma}_{\pm}$ are standard qubit raising/lowering operators,
$\mathcal{D}[x]\hat{\rho}=x\rho
x^{\dagger}-\left\\{x^{\dagger}x,\hat{\rho}\right\\}/2$ is the usual Lindblad
dissipator, $\gamma$ 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 Gardiner (1986); Murch  _et al._ (2013)
and modifications of the Mollow triplet fluorescence spectrum Carmichael  _et
al._ (1987a, b); Dalton  _et al._ (1999); Toyli  _et al._ (2016).

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 Kraus and Cirac (2004). The
evolution is described by

$\displaystyle\dot{\hat{\rho}}_{q}$
$\displaystyle=\gamma_{1}\mathcal{D}\left[\cosh(r)\hat{\sigma}_{-}^{(1)}+\sinh(r)\hat{\sigma}_{+}^{(2)}\right]\hat{\rho}_{q}$
$\displaystyle+\gamma_{2}\mathcal{D}\left[\sinh(r)\hat{\sigma}_{+}^{(1)}+\cosh(r)\hat{\sigma}_{-}^{(2)}\right]\hat{\rho}_{q},$
(2)

where $\gamma_{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

$\displaystyle\ket{\Psi}_{q}=\frac{1}{\sqrt{\cosh(2r)}}\left(\cosh(r)\ket{00}-\sinh(r)\ket{11}\right).$
(3)

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 two-
qubit 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 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. This second approach
both generalizes and provides intuition for a protocol first introduced in
Ref. Motzoi  _et al._ (2016). 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 Cirac  _et al._
(1993), optomechanics Kronwald  _et al._ (2013) and circuit QED Didier  _et
al._ (2014), and bosonic optomechanical entanglement Woolley and Clerk (2013,
2014). Related master equations have also been studied in the setting of spin
ensembles Parkins  _et al._ (2006); Dalla Torre  _et al._ (2013), 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

Figure 2: Pathways to the dissipative dynamics of Eq. (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 $\omega_{c}\pm\Omega$ interfere such that both qubit decay (red
sideband) and excitation (blue sideband) result in an emitted cavity photon at
frequency $\omega_{c}$.

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:

$\displaystyle\hat{H}=\omega_{c}\hat{a}^{\dagger}\hat{a}+\frac{\Omega}{2}\hat{\sigma}_{z}+g(t)\left(\hat{a}+\hat{a}^{\dagger}\right)\hat{\sigma}_{x}.$
(4)

This is the usual Rabi Hamiltonian describing a cavity of frequency
$\omega_{c}$ coupled to a qubit of frequency $\Omega$, with a time-dependent
coupling $g(t)$. We consider phase-locked double-sideband cosine modulation of
the coupling at frequencies $\omega_{c}\pm\Omega$, described by

$\displaystyle
g(t)=\bar{g}\left(\alpha_{+}e^{-i(\omega_{c}+\Omega)t}+\alpha_{-}e^{-i(\omega_{c}-\Omega)t}\right)+h.c.,$
(5)

with $\alpha_{\rm\pm}$ the complex-valued amplitudes of the modulation. In
circuit QED systems, this parametrically modulated coupling has been directly
demonstrated using a superconducting quantum interference device (SQUID) Chen
_et al._ (2014); McKay  _et al._ (2016); Lu  _et al._ (2017); Zhong  _et al._
(2019); Chang  _et al._ (2020); Noh  _et al._ (2021); Brown  _et al._ (2021),
and alternatively implemented using a dispersive qubit-cavity coupling, a
qubit Rabi drive, and classical cavity displacements Hacohen-Gourgy  _et al._
(2016); Eddins  _et al._ (2018).

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 $\Delta=\omega_{c}-\Omega$, and obtain the Hamiltonian

$\displaystyle\hat{H}^{\prime}$
$\displaystyle=\bar{g}\hat{a}^{\dagger}\left(\alpha_{+}\hat{\sigma}_{+}+\alpha_{-}\hat{\sigma}_{-}\right)+h.c.,$
(6)

where we require that $\absolutevalue{\alpha_{\pm}}\ll 2\Delta/\bar{g}$ for
the RWA. Assuming that the cavity decays at a rate
$\kappa\gg\absolutevalue{\alpha_{\pm}}\bar{g}$ into a zero temperature
environment, adiabatic elimination of the cavity yields a qubit-only master
equation

$\displaystyle\dot{\hat{\rho}}_{q}=\frac{4\bar{g}^{2}}{\kappa}\mathcal{D}\left[\alpha_{+}\hat{\sigma}_{+}+\alpha_{-}\hat{\sigma}_{-}\right]\hat{\rho}_{q}.$
(7)

With the associations $\alpha_{+}=\alpha\sinh(r)$ and
$\alpha_{-}=\alpha\cosh(r)$ we obtain Eq. (1) with decay rate
$\gamma=4\bar{g}^{2}\alpha^{2}/\kappa$ and squeezing parameter $r$ set by
$\absolutevalue{\alpha_{+}/\alpha_{-}}=\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 Lu
_et al._ (2017); Huang  _et al._ (2018), 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.
(1) is a direct consequence of the photon-number pairing in the squeezed
state. To be explicit, if the reservoir ends up in a state with an odd number
of photons $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 $\omega_{c}+\Omega$) with
the red sideband Raman process (driven by the modulation at frequency
$\omega_{c}-\Omega$). 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

Figure 3: Two-qubit interference. a) and b) interference between the qubit
decay and excitation processes mediated by the coupling modulation that
results in the combined dissipators $\hat{J}_{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
$\Omega_{1}$ ($\Omega_{2}$).

The above ideas can easily be extended to two qubits (frequencies $\Omega_{1}$
and $\Omega_{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 equation

$\displaystyle\dot{\hat{\rho}}_{q}=\frac{4\bar{g}^{2}}{\kappa}\Big{(}\mathcal{D}\left[\hat{\mathcal{J}}_{a}\right]+\mathcal{D}\left[\hat{\mathcal{J}}_{b}\right]\Big{)}\hat{\rho}_{q},$
(8)

with

$\displaystyle\hat{\mathcal{J}}_{a}=\alpha_{-}\hat{\sigma}^{(1)}_{-}+\alpha_{+}\hat{\sigma}^{(2)}_{+},\leavevmode\nobreak\
\leavevmode\nobreak\
\hat{\mathcal{J}}_{b}=\beta_{+}\hat{\sigma}^{(1)}_{+}+\beta_{-}\hat{\sigma}^{(2)}_{-},$
(9)

describing the independent qubit dissipative channels mediated by the bosonic
modes shown in Fig. 3a) and 3b). Here $\bar{g}$ and $\alpha_{\pm}$
($\beta_{\pm}$) define the qubit-mode coupling modulation for mode $a$ ($b$),
in direct analogy to Eq. (5). By setting $\alpha_{-}=\beta_{-}=\alpha\cosh(r)$
and $\alpha_{+}=\beta_{+}=\alpha\sinh(r)$ we obtain Eq. (2) with
$\gamma_{1}=\gamma_{2}=4\bar{g}^{2}\alpha^{2}/\kappa$. 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

$\displaystyle\hat{H}_{\rm
coup}=\left[g_{1}(t)\left(\hat{a}+\hat{a}^{\dagger}\right)+g_{2}(t)\left(\hat{b}+\hat{b}^{\dagger}\right)\right]\left(\hat{\sigma}^{(1)}_{x}+\hat{\sigma}^{(2)}_{x}\right).$
(10)

The coupling $g_{1}(t)$ is modulated such that the $a$ cavity mediates the
interference process of Fig. 3a) described by the dissipator
$\hat{\mathcal{J}}_{a}$, and similarly $g_{2}(t)$ and the $b$ cavity mediate
Fig. 3b) and $\hat{\mathcal{J}}_{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 dissipators

$\displaystyle\hat{\mathcal{J}}^{f}_{a}=\alpha_{-}\hat{\sigma}^{(1)}_{-}+\alpha_{+}\hat{\sigma}^{(2)}_{-},\leavevmode\nobreak\
\leavevmode\nobreak\
\hat{\mathcal{J}}^{f}_{b}=\beta_{+}\hat{\sigma}^{(1)}_{+}+\beta_{-}\hat{\sigma}^{(2)}_{+},$
(11)

which leads to Eq. (2) with the state of qubit 2 flipped.

The second setup uses a single cavity with a broad linewidth ($\kappa$) 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

$\displaystyle\hat{H}_{\rm
coup}=g(t)\left(\hat{a}+\hat{a}^{\dagger}\right)\left(\hat{\sigma}^{(1)}_{x}+\hat{\sigma}^{(2)}_{x}\right).$
(12)

In this setup, as a result of the modulated coupling, one pair of sideband
processes results in cavity photon emission at frequency
$\omega_{a}+\delta\Omega$, and the other pair results in emission at
$\omega_{a}-\delta\Omega$, with $\delta\Omega=(\Omega_{1}-\Omega_{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 sum-frequency
drives a Raman process where qubit 1 is excited and a photon of frequency
$\Omega_{2}$ is emitted into the TL. Due to the static coupling, qubit 2 also
decays into the TL at frequency $\Omega_{2}$. These two processes interfere,
as do the analogous processes at emission frequency $\Omega_{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 Zhong  _et al._
(2019); Chang  _et al._ (2020), 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
qubit-waveguide 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. Motzoi
_et al._ (2016).

Consider two qubits which are coupled to a common source of loss at a rate
$\Gamma$. The dissipation is described by the master equation
$\dot{\hat{\rho}}_{q}=\Gamma\mathcal{D}\left[\hat{\mathcal{J}}\right]\hat{\rho}_{q}$
with the collective loss operator

$\hat{\mathcal{J}}=\hat{\sigma}_{-}^{(1)}+\eta\hat{\sigma}_{-}^{(2)},$ (13)

where $\eta$ allows for an asymmetry between the qubit-bath couplings. For
clarity we take $\eta=1$ in what follows; as we show later, analogous results
hold for $\eta\neq 1$. $\hat{\mathcal{J}}$ has a two-parameter family of dark
states of the form

$\ket{\Phi[\alpha,\phi]}=\sqrt{1-\alpha^{2}}\ket{00}+\frac{e^{i\phi}\alpha}{\sqrt{2}}(\ket{01}-\ket{10}),$
(14)

where $0\leq\alpha\leq 1$. These states are entangled for any $\alpha>0$,
whereas for $\alpha=0$ we have the trivial dark state $\ket{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 $\alpha>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 $\phi=0$; in this case we have the
explicit mapping

$\displaystyle\ket{\Psi[\alpha]}_{q}$
$\displaystyle\equiv\hat{U}[\alpha]\ket{\Phi[\alpha]}$
$\displaystyle=\frac{1}{\sqrt{\cosh(2r)}}\left(\cosh(r)\ket{00}-\sinh(r)\ket{11}\right),$
(15)

where $r[\alpha]$ is determined via $\alpha=\sqrt{\tanh(2r)}$. Crucially, the
required unitary is local, i.e.
$\hat{U}[\alpha]=\hat{U}^{(1)}[\alpha]\otimes\hat{U}^{(2)}[\alpha]$. It
corresponds to opposite local rotations of each qubit about the y-axis

$\displaystyle\hat{U}^{(j)}[\alpha]=\exp(\pm
i\frac{\theta}{2}\hat{\sigma}_{y}^{(j)}),$ (16)

where qubit 1 (2) takes the upper (lower) sign. The rotation angle $\theta$ is
a function of $\alpha$ (or equivalently $r$), and is given by

$\cos\theta=\sqrt{\frac{1-\alpha^{2}}{1+\alpha^{2}}}=e^{-2r}.$ (17)

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 $\hat{U}$
and $\ket{\Psi}_{q}$ to be functions of $r$ rather than $\alpha$. The next
step in our mapping to synthetic two-mode squeezing is to examine the form of
the collective loss dissipator $\hat{\mathcal{J}}$ in the new frame defined by
$\hat{U}[r_{0}]$. We have

$\hat{\mathcal{J}}^{\prime}\equiv\hat{U}[r_{0}]\hat{\mathcal{J}}\hat{U}^{\dagger}[r_{0}]=\hat{\mathcal{J}}^{\prime}_{1}+\hat{\mathcal{J}}^{\prime}_{2}+\hat{\mathcal{J}}^{\prime}_{Z}.$
(18)

Here, $\hat{\mathcal{J}}_{1}^{\prime}$ and $\hat{\mathcal{J}}_{2}^{\prime}$
are precisely the two squeezing dissipators needed for synthetic squeezing
(c.f. Eq. 2))

$\displaystyle\hat{\mathcal{J}}_{1}^{\prime}$
$\displaystyle=e^{-r_{0}}\left(\cosh(r_{0})\hat{\sigma}_{-}^{(1)}-\sinh(r_{0})\hat{\sigma}_{+}^{(2)}\right),$
(19) $\displaystyle\hat{\mathcal{J}}_{2}^{\prime}$
$\displaystyle=e^{-r_{0}}\left(\cosh(r_{0})\hat{\sigma}_{-}^{(2)}-\sinh(r_{0})\hat{\sigma}_{+}^{(1)}\right).$
(20)

Further, $\hat{\mathcal{J}}_{Z}^{\prime}$ is a correlated dephasing term

$\hat{\mathcal{J}}_{Z}^{\prime}=\frac{1}{2}\sqrt{1-e^{-4r_{0}}}\left(\hat{\sigma}_{z}^{(1)}-\hat{\sigma}_{z}^{(2)}\right).$
(21)

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): $\mathcal{J}^{\prime}$ is the sum of these
dissipators, plus a dephasing dissipator. One can easily check that the paired
target state $\ket{\Psi[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 $\ket{00}$)
remains a dark state.

However, in our new frame defined by $\hat{U}[r_{0}]$, the state
$\ket{\Psi[r_{0}]}_{q}$ is the only dark state of the dissipator
$\hat{\mathcal{J}}^{\prime}$ that does not have any single-excitation
components, i.e. kets $\ket{01}$ or $\ket{10}$. This suggests a very simple
strategy to make the paired state $\ket{\Psi[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

$\hat{H}^{\prime}=\frac{\mu}{2}\left(\hat{\sigma}_{z}^{(1)}-\hat{\sigma}_{z}^{(2)}\right),$
(22)

we introduce an energy gap $\mu$ between the paired states $\ket{00}$ and
$\ket{11}$ (which have zero energy) and the unpaired states $\ket{10}$ and
$\ket{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 $\mu>0$,
the target entangled state $\ket{\Psi[r_{0}]}_{q}$ is the unique steady state
of the dynamics generated by our master equation, which in the new frame has
the form

$\dot{\hat{\rho}}^{\prime}=-i[\hat{H}^{\prime},\hat{\rho}^{\prime}]+\Gamma\mathcal{D}\left[\hat{\mathcal{J}}^{\prime}_{1}+\hat{\mathcal{J}}^{\prime}_{2}+\hat{\mathcal{J}}^{\prime}_{Z}\right]\hat{\rho}^{\prime}.$
(23)

Here, $\hat{\rho}^{\prime}=\hat{U}[r_{0}]\hat{\rho}\hat{U}^{-1}[r_{0}]$ is the
transformed-frame density matrix of the two qubits.

Note that for large $\mu\gg\Gamma$, the connection to the synthetic squeezing
master equation of Eq. (2) is even more direct. In this limit,
$\hat{H}^{\prime}$ will cause $\hat{\mathcal{J}}^{\prime}_{1/2}$ to oscillate
at frequencies $\pm\mu$, while $\hat{\mathcal{J}}^{\prime}_{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
approximation

$\dot{\hat{\rho}}^{\prime}\simeq-i[\hat{H}^{\prime},\hat{\rho}^{\prime}]+\Gamma\left(\mathcal{D}\left[\hat{\mathcal{J}}^{\prime}_{1}\right]+\mathcal{D}\left[\hat{\mathcal{J}}^{\prime}_{2}\right]+\mathcal{D}\left[\hat{\mathcal{J}}^{\prime}_{Z}\right]\right)\hat{\rho}^{\prime}.$
(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, $\ket{\Psi[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
$\mu\gg\Gamma$. In this limit we can also remove the Hamiltonian term by
moving to the interaction frame,
$\hat{\rho}^{\prime\prime}=\exp\small(-i\hat{H}^{\prime}t\small)\hat{\rho}^{\prime}\exp\small(i\hat{H}^{\prime}t\small)$,
where the equivalence to Eq. (2) is more obvious, and $\ket{\Psi[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 $\hat{U}^{\dagger}[r_{0}]$). The master
equation in the original frame has the form

$\dot{\hat{\rho}}=-i[\hat{H},\hat{\rho}]+\Gamma\mathcal{D}\left[\hat{\mathcal{J}}\right]\hat{\rho},$
(25)

where $\hat{\mathcal{J}}$ is the simple collective loss dissipator in Eq.
(13), and the lab-frame Hamiltonian
$\hat{H}=\hat{U}^{\dagger}[r_{0}]\hat{H}^{\prime}\hat{U}[r_{0}]$ has the form

$\hat{H}=\frac{\Delta}{2}\left(\hat{\sigma}_{z}^{(1)}-\hat{\sigma}_{z}^{(2)}\right)+\Lambda\left(\hat{\sigma}_{x}^{(1)}+\hat{\sigma}_{x}^{(2)}\right).$
(26)

with

$\displaystyle\Delta=\mu
e^{-2r_{0}},\,\,\,\,\,\Lambda=(\mu/2)\sqrt{1-e^{-4r_{0}}}.$ (27)

This Hamiltonian describes a two-qubit system subject to Rabi drives at the
same frequency $\nu$ and amplitude $\Lambda$, in the rotating frame determined
by the drive frequency. Thus, $\pm\Delta$ describe the detuning of the drive
from each qubit, such that $\nu$ 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

$\hat{H}_{\rm
lab}=\frac{1}{2}\left(\Omega_{1}\hat{\sigma}_{z}^{(1)}+\Omega_{2}\hat{\sigma}_{z}^{(2)}\right)+\Lambda\left[e^{i\nu
t}\left(\hat{\sigma}_{-}^{(1)}+\hat{\sigma}_{-}^{(2)}\right)+\textrm{h.c.}\right]$
(28)

with $\nu=(\Omega_{1}+\Omega_{2})/2$ and $\Delta=\Omega_{1}-\Omega_{2}$.

We note that the master equation in Eq. (25) was first presented in Ref.
Motzoi  _et al._ (2016). 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. $\eta\neq 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:

$\hat{H}=\frac{\Delta+\epsilon}{2}\hat{\sigma}_{z}^{(1)}-\frac{\Delta-\epsilon}{2}\hat{\sigma}_{z}^{(2)}+\Lambda\left(\eta\hat{\sigma}_{x}^{(1)}+\hat{\sigma}_{x}^{(2)}\right),$
(29)

where

$\epsilon=\frac{\Lambda^{2}}{\Delta}\left(1-\eta^{2}\right),$ (30)

and where $\Delta,\leavevmode\nobreak\ \Lambda$ are no longer given by Eq.
(27) but are determined via Eqs. (96) and (97) 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 $\eta=1$ case
and a large squeezing parameter $r$, the Rabi-drive amplitude scales as
$\Lambda\sim e^{2r}(\Omega_{1}-\Omega_{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 $\hat{H}_{\rm wg}=J(\hat{\sigma}_{+}^{(1)}\hat{\sigma}_{-}^{(2)}+{\rm
h.c.})$. As discussed in Ref. Motzoi  _et al._ (2016), 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 $\hat{H}_{\rm wg}$ is the $\alpha=1$ dark state of the family
of Eq. (14) (equivalently the $r\to\infty$ 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. Motzoi  _et al._ (2016) 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\rightarrow\infty$ 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$,
$\alpha$ can be decreased to keep $\alpha_{\pm}$ and $\beta{\pm}$ constant. As
in the single qubit case, this comes at the cost of reducing the overall decay
rate $\Gamma\equiv 4\bar{g}^{2}\alpha^{2}/\kappa$. Keeping the modulation
amplitudes fixed, more entanglement requires a longer preparation time, a
trade-off also observed in similar bosonic schemes Woolley and Clerk (2014).

Moreover, saturating the $r\rightarrow\infty$ 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\rightarrow\infty$ limit of Eq. (11)

$\displaystyle\dot{\hat{\rho}}_{q}$
$\displaystyle=\gamma\left(\mathcal{D}\left[\hat{\sigma}_{1}^{-}+\hat{\sigma}_{2}^{-}\right]+\mathcal{D}\left[\hat{\sigma}_{1}^{+}+\hat{\sigma}_{2}^{+}\right]\right)\hat{\rho}_{q},$
(31)

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).

###  V.1 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
$\dot{\hat{\rho}}_{q}=\mathcal{L}(\hat{\rho}_{q})$, where $\mathcal{L}$
corresponds to Eq. (8), we can vectorize this equation as
$\dot{\vec{\rho}}_{q}=\hat{\mathcal{L}}\vec{\rho}_{q}$, where the super-
operator matrix $\hat{\mathcal{L}}$ 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

$\displaystyle\kappa_{\rm
slow}=\frac{\Gamma}{3\sinh^{2}(r)}+\mathcal{O}\left(\frac{1}{\sinh^{4}(r)}\right).$
(32)

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 trade-off 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 Kraus and Cirac (2004), 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 Woolley and Clerk (2014), where $\hat{\sigma}_{1/2}^{-}$
are replaced by bosonic modes $\hat{d}_{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\gg\leavevmode\nobreak\ 1/2$. The bosonic model is often used as an
approximation to such a large spin model via the Holstein-Primakoff (HP)
transformation Holstein and Primakoff (1940), which connects the spin operator
$\hat{S}^{+}$ to a bosonic mode via
$\hat{S}^{+}=\hat{d}^{\dagger}\sqrt{2S-\hat{d}^{\dagger}\hat{d}}$. In the HP
representation, the large spin version of Eq. (8) has dissipators

$\displaystyle\hat{\mathcal{J}}_{a}=\cosh(r)\sqrt{1-\frac{\hat{n}_{1}}{2S}}\leavevmode\nobreak\
\hat{d}_{1}+\sinh(r)\hat{d}_{2}^{\dagger}\sqrt{1-\frac{\hat{n}_{2}}{2S}},$
(33)
$\displaystyle\hat{\mathcal{J}}_{b}=\cosh(r)\sqrt{1-\frac{\hat{n}_{2}}{2S}}\leavevmode\nobreak\
\hat{d}_{2}+\sinh(r)\hat{d}_{1}^{\dagger}\sqrt{1-\frac{\hat{n}_{1}}{2S}},$
(34)

with $\hat{n}_{k}=\hat{d}_{k}^{\dagger}\hat{d}_{k}$, and an enhanced
dissipative rate $\Gamma_{S}\equiv 16S\bar{g}^{2}\alpha^{2}/\kappa$. In the
limit $S\gg\left<\hat{n}_{k}\right>$ the square root factors can be neglected
and the bosonic model is obtained.

However, the steady-state of the bosonic model has
$\left<\hat{n}_{k}\right>=\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

$\displaystyle\Gamma^{\rm
MF}_{S}=\frac{16S\bar{g}^{2}\alpha^{2}}{\kappa}\left(1-\frac{\sinh^{2}(r)}{2S}\right).$
(35)

This rate slows as $r$ increases, similar to what was observed in the two-
qubit ($S=1/2$) model. Dynamically, $\Gamma^{\rm MF}_{S}$ will decrease as
$\left<\hat{n}_{k}\right>$ 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.

###  V.2 Thermal Dissipative Squeezing

Figure 4: Thermal dissipative squeezing. Concurrence (black solid line) and
purity (red dashed line) of the steady-state of thermal dissipative squeezing,
cf. Eq. (36), as a function of the temperature of the transmission line. The
inset shows the concurrence as a function of the squeezing parameter $r$ for
zero temperature, and the dashed horizontal line indicates the squeezed
parameter value used in the main plot.

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 equation

$\displaystyle\dot{\hat{\rho}}_{q}=\leavevmode\nobreak\ $
$\displaystyle\Gamma\Big{(}\left(1+\bar{n}_{\rm
th}\right)\mathcal{D}\left[\hat{\mathcal{J}}_{a}\right]+\bar{n}_{\rm
th}\mathcal{D}\left[\hat{\mathcal{J}}^{\dagger}_{a}\right]$
$\displaystyle+\left(1+\bar{n}_{\rm
th}\right)\mathcal{D}\left[\hat{\mathcal{J}}_{b}\right]+\bar{n}_{\rm
th}\mathcal{D}\left[\hat{\mathcal{J}}^{\dagger}_{b}\right]\Big{)}\hat{\rho}_{q},$
(36)

with $\hat{\mathcal{J}}_{a}$ and $\hat{\mathcal{J}}_{b}$ as before for
$\alpha_{-}=\beta_{-}=\alpha\cosh(r)$, $\alpha_{+}=\beta_{+}=\alpha\sinh(r)$,
and where $\bar{n}_{\rm 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 steady-
state 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 Wang  _et al._ (2019), 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 Chen  _et al._
(2014); McKay  _et al._ (2016); Lu  _et al._ (2017); Zhong  _et al._ (2019);
Chang  _et al._ (2020), 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
Hacohen-Gourgy  _et al._ (2016); Eddins  _et al._ (2018). 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 Yurke  _et al._ (1989); Castellanos-Beltran  _et al._ (2008). 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 Leghtas  _et al._ (2013); Shankar
_et al._ (2013); Aron  _et al._ (2014); Kimchi-Schwartz  _et al._ (2016);
Motzoi  _et al._ (2015, 2016); Doucet  _et al._ (2020); Brown  _et al._
(2021). 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. Motzoi  _et al._ (2016). For schemes with no slow rate, such as
Ref. Doucet  _et al._ (2020), understanding where the connection to
dissipative squeezing breaks is a topic worthy of future study.

###### Acknowledgements.

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).

## References

  * Gardiner and Zoller (2004) C. W. Gardiner and P. Zoller,  _Quantum Noise_ (Springer, Heidelberg, 2004).
  * Walls (1983) D. F. Walls, Nature 306, 141 (1983).
  * Giovannetti  _et al._ (2004) V. Giovannetti, S. Lloyd, and L. Maccone, Science 306, 1330 (2004).
  * Caves  _et al._ (1980) C. M. Caves, K. S. Thorne, R. W. P. Drever, V. D. Sandberg, and M. Zimmermann, Rev. Mod. Phys. 52, 341 (1980).
  * Braginsky and Khalili (1992) V. B. Braginsky and F. Y. Khalili,  _Quantum Measurement_ (Cambridge University Press, New York, NY, USA, 1992).
  * LIGO Scientific Collaboration (2013) LIGO Scientific Collaboration, Nature Photonics 7, 613 (2013).
  * Iwasawa  _et al._ (2013) K. Iwasawa, K. Makino, H. Yonezawa, M. Tsang, A. Davidovic, E. Huntington, and A. Furusawa, Phys. Rev. Lett. 111, 163602 (2013).
  * Peano  _et al._ (2015) V. Peano, H. G. L. Schwefel, C. Marquardt, and F. Marquardt, Phys. Rev. Lett. 115, 243603 (2015).
  * Taylor  _et al._ (2013) M. A. Taylor, J. Janousek, V. Daria, J. Knittel, B. Hage, BachorHans-A., and W. P. Bowen, Nat Photon 7, 229 (2013).
  * Barzanjeh  _et al._ (2014) S. Barzanjeh, D. P. DiVincenzo, and B. M. Terhal, Physical Review B 90, 134515 (2014).
  * Didier  _et al._ (2015a) N. Didier, A. Kamal, W. D. Oliver, A. Blais, and A. A. Clerk, Phys. Rev. Lett. 115, 093604 (2015a).
  * Didier  _et al._ (2015b) N. Didier, J. Bourassa, and A. Blais, Phys. Rev. Lett. 115, 203601 (2015b).
  * Govia and Clerk (2017) L. C. G. Govia and A. A. Clerk, New Journal of Physics 19, 023044 (2017).
  * Eddins  _et al._ (2018) A. Eddins, S. Schreppler, D. M. Toyli, L. S. Martin, S. Hacohen-Gourgy, L. C. G. Govia, H. Ribeiro, A. A. Clerk, and I. Siddiqi, Phys. Rev. Lett. 120, 040505 (2018).
  * Puri and Blais (2016) S. Puri and A. Blais, Phys. Rev. Lett. 116, 180501 (2016).
  * Royer  _et al._ (2017) B. Royer, A. L. Grimsmo, N. Didier, and A. Blais, Quantum 1, 11 (2017).
  * Gardiner (1986) C. W. Gardiner, Phys. Rev. Lett. 56, 1917 (1986).
  * Murch  _et al._ (2013) K. W. Murch, S. J. Weber, K. M. Beck, E. Ginossar, and I. Siddiqi, Nature 499, 62 (2013).
  * Carmichael  _et al._ (1987a) H. J. Carmichael, A. S. Lane, and D. F. Walls, Phys. Rev. Lett. 58, 2539 (1987a).
  * Carmichael  _et al._ (1987b) H. Carmichael, A. Lane, and D. Walls, J. Mod. Opt. 34, 821 (1987b).
  * Dalton  _et al._ (1999) B. J. Dalton, Z. Ficek, and S. Swain, J. Mod. Opt. 46, 379 (1999).
  * Toyli  _et al._ (2016) D. M. Toyli, A. W. Eddins, S. Boutin, S. Puri, D. Hover, V. Bolkhovsky, W. D. Oliver, A. Blais, and I. Siddiqi, Phys. Rev. X 6, 031004 (2016).
  * Kraus and Cirac (2004) B. Kraus and J. I. Cirac, Phys. Rev. Lett. 92, 013602 (2004).
  * Motzoi  _et al._ (2016) F. Motzoi, E. Halperin, X. Wang, K. B. Whaley, and S. Schirmer, Phys. Rev. A 94, 032313 (2016).
  * Cirac  _et al._ (1993) J. I. Cirac, A. S. Parkins, R. Blatt, and P. Zoller, Phys. Rev. Lett. 70, 556 (1993).
  * Kronwald  _et al._ (2013) A. Kronwald, F. Marquardt, and A. A. Clerk, Phys. Rev. A 88, 063833 (2013).
  * Didier  _et al._ (2014) N. Didier, F. Qassemi, and A. Blais, Phys. Rev. A 89, 013820 (2014).
  * Woolley and Clerk (2013) M. J. Woolley and A. A. Clerk, Phys. Rev. A 87, 063846 (2013).
  * Woolley and Clerk (2014) M. J. Woolley and A. A. Clerk, Phys. Rev. A 89, 063805 (2014).
  * Parkins  _et al._ (2006) A. S. Parkins, E. Solano, and J. I. Cirac, Phys. Rev. Lett. 96, 053602 (2006).
  * Dalla Torre  _et al._ (2013) E. G. Dalla Torre, J. Otterbach, E. Demler, V. Vuletic, and M. D. Lukin, Phys. Rev. Lett. 110, 120402 (2013).
  * Chen  _et al._ (2014) Y. Chen, C. Neill, P. Roushan, N. Leung, M. Fang, R. Barends, J. Kelly, B. Campbell, Z. Chen, B. Chiaro, A. Dunsworth, E. Jeffrey, A. Megrant, J. Y. Mutus, P. J. J. O’Malley, C. M. Quintana, D. Sank, A. Vainsencher, J. Wenner, T. C. White, M. R. Geller, A. N. Cleland, and J. M. Martinis, Phys. Rev. Lett. 113, 220502 (2014).
  * McKay  _et al._ (2016) D. C. McKay, S. Filipp, A. Mezzacapo, E. Magesan, J. M. Chow, and J. M. Gambetta, Phys. Rev. Applied 6, 064007 (2016).
  * Lu  _et al._ (2017) Y. Lu, S. Chakram, N. Leung, N. Earnest, R. K. Naik, Z. Huang, P. Groszkowski, E. Kapit, J. Koch, and D. I. Schuster, Phys. Rev. Lett. 119, 150502 (2017).
  * Zhong  _et al._ (2019) Y. P. Zhong, H. S. Chang, K. J. Satzinger, M. H. Chou, A. Bienfait, C. R. Conner, É. Dumur, J. Grebel, G. A. Peairs, R. G. Povey, D. I. Schuster, and A. N. Cleland, Nat. Phys. 15, 741 (2019).
  * Chang  _et al._ (2020) H.-S. Chang, Y. P. Zhong, A. Bienfait, M.-H. Chou, C. R. Conner, E. Dumur, J. Grebel, G. A. Peairs, R. G. Povey, K. J. Satzinger, and A. N. Cleland, Phys. Rev. Lett. 124, 240502 (2020).
  * Noh  _et al._ (2021) T. Noh, Z. Xiao, K. Cicak, X. Y. Jin, E. Doucet, J. Teufel, J. Aumentado, L. C. G. Govia, L. Ranzani, A. Kamal, and R. W. Simmonds, “Strong parametric dispersive shifts in a statically decoupled multi-qubit cavity qed system,” (2021), arXiv:2103.09277 .
  * Brown  _et al._ (2021) T. Brown, E. Doucet, D. Ristè, G. Ribeill, K. Cicak, J. Aumentado, R. Simmonds, L. Govia, A. Kamal, and L. Ranzani, “Trade off-free entanglement stabilization in a superconducting qutrit-qubit system,” (2021), arXiv:2107.13579 .
  * Hacohen-Gourgy  _et al._ (2016) S. Hacohen-Gourgy, L. S. Martin, E. Flurin, V. V. Ramasesh, K. B. Whaley, and I. Siddiqi, Nature 538, 491 (2016).
  * Huang  _et al._ (2018) Z. Huang, Y. Lu, E. Kapit, D. I. Schuster, and J. Koch, Phys. Rev. A 97, 062345 (2018).
  * Holstein and Primakoff (1940) T. Holstein and H. Primakoff, Phys. Rev. 58, 1098 (1940).
  * Wang  _et al._ (2019) Z. Wang, S. Shankar, Z. Minev, P. Campagne-Ibarcq, A. Narla, and M. Devoret, Phys. Rev. Applied 11, 014031 (2019).
  * Yurke  _et al._ (1989) B. Yurke, L. R. Corruccini, P. G. Kaminsky, L. W. Rupp, A. D. Smith, A. H. Silver, R. W. Simon, and E. A. Whittaker, Phys. Rev. A 39, 2519 (1989).
  * Castellanos-Beltran  _et al._ (2008) M. A. Castellanos-Beltran, K. D. Irwin, G. C. Hilton, L. R. Vale, and K. W. Lehnert, Nat. Phys 4, 929 (2008).
  * Leghtas  _et al._ (2013) Z. Leghtas, U. Vool, S. Shankar, M. Hatridge, S. M. Girvin, M. H. Devoret, and M. Mirrahimi, Phys. Rev. A 88, 023849 (2013).
  * Shankar  _et al._ (2013) S. Shankar, M. Hatridge, Z. Leghtas, K. M. Sliwa, A. Narla, U. Vool, S. M. Girvin, L. Frunzio, M. Mirrahimi, and M. H. Devoret, Nature 504, 419 (2013).
  * Aron  _et al._ (2014) C. Aron, M. Kulkarni, and H. E. Türeci, Phys. Rev. A 90, 062305 (2014).
  * Kimchi-Schwartz  _et al._ (2016) M. E. Kimchi-Schwartz, L. Martin, E. Flurin, C. Aron, M. Kulkarni, H. E. Tureci, and I. Siddiqi, Phys. Rev. Lett. 116, 240503 (2016).
  * Motzoi  _et al._ (2015) F. Motzoi, K. B. Whaley, and M. Sarovar, Phys. Rev. A 92, 032308 (2015).
  * Doucet  _et al._ (2020) E. Doucet, F. Reiter, L. Ranzani, and A. Kamal, Phys. Rev. Research 2, 023370 (2020).
  * Breuer and Petruccione (2002) H.-P. Breuer and F. Petruccione,  _The Theory of Open Quantum Systems_ (Oxford University Press, Oxford, 2002).

##  Appendix A Two Qubit Master Equation Derivations

###  A.1 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

$\displaystyle\hat{H}$
$\displaystyle=\omega_{a}\hat{a}^{\dagger}\hat{a}+\frac{\Omega_{1}}{2}\hat{\sigma}^{(1)}_{z}+g_{1}(t)\left(\hat{a}+\hat{a}^{\dagger}\right)\left(\hat{\sigma}^{(1)}_{x}+\hat{\sigma}^{(2)}_{x}\right)$
$\displaystyle+\omega_{b}\hat{b}^{\dagger}\hat{b}+\frac{\Omega_{2}}{2}\hat{\sigma}^{(2)}_{z}+g_{2}(t)\left(\hat{b}+\hat{b}^{\dagger}\right)\left(\hat{\sigma}^{(1)}_{x}+\hat{\sigma}^{(2)}_{x}\right),$
(37)

with the modulated couplings

$\displaystyle
g_{1}(t)=\bar{g}\left(\alpha_{-}e^{-i\left(\omega_{a}-\Omega_{1}\right)t}+\alpha_{+}e^{-i\left(\omega_{a}+\Omega_{2}\right)t}\right)+h.c.,$
$\displaystyle
g_{2}(t)=\bar{g}\left(\beta_{-}e^{-i\left(\omega_{b}-\Omega_{2}\right)t}+\beta_{+}e^{-i\left(\omega_{b}+\Omega_{1}\right)t}\right)+h.c..$

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 $\hat{a}$ and $\hat{b}$ in Eq. (37).

Following the same procedure as for one qubit, we move to the interaction
frame for Eq. (37) and drop all fast-oscillating terms to obtain

$\displaystyle\hat{H}^{\prime}$
$\displaystyle=\bar{g}\hat{a}^{\dagger}\left(\alpha_{+}\hat{\sigma}^{(2)}_{+}+\alpha_{-}\hat{\sigma}^{(1)}_{-}\right)+h.c.$
$\displaystyle+\bar{g}\hat{b}^{\dagger}\left(\beta_{+}\hat{\sigma}^{(1)}_{+}+\beta_{-}\hat{\sigma}^{(2)}_{-}\right)+h.c..$
(38)

where we have assumed that $\absolutevalue{\Omega_{1}-\Omega_{2}}\gg\bar{g}$
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
$\kappa_{a},\kappa_{b}\gg\absolutevalue{\alpha_{\pm}}\bar{g},\absolutevalue{\beta_{\pm}}\bar{g}$.
In the case where $\kappa_{a}=\kappa_{b}=\kappa$, this gives

$\displaystyle\dot{\hat{\rho}}_{q}=\frac{4\bar{g}^{2}}{\kappa}\Big{(}\mathcal{D}\left[\hat{\mathcal{J}}_{a}\right]+\mathcal{D}\left[\mathcal{J}_{b}\right]\Big{)}\hat{\rho}_{q},$
(39)

with

$\displaystyle\hat{\mathcal{J}}_{a}=\alpha_{-}\hat{\sigma}^{(1)}_{-}+\alpha_{+}\hat{\sigma}^{(2)}_{+},$
(40)
$\displaystyle\hat{\mathcal{J}}_{b}=\beta_{+}\hat{\sigma}^{(1)}_{+}+\beta_{-}\hat{\sigma}^{(2)}_{-},$
(41)

as depicted in Fig. 3a) and b). Starting with equal coupling $\bar{g}$, and
setting $\kappa_{a}=\kappa_{b}=\kappa$ 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
$\hat{a}$, and blue sideband tones on mode $\hat{b}$, i.e.

$\displaystyle
g_{1}(t)=\bar{g}\left(\alpha_{-}e^{-i\left(\omega_{a}-\Omega_{1}\right)t}+\alpha_{+}e^{-i\left(\omega_{a}-\Omega_{2}\right)t}\right)+h.c.,$
$\displaystyle
g_{2}(t)=\bar{g}\left(\beta_{-}e^{-i\left(\omega_{b}+\Omega_{2}\right)t}+\beta_{+}e^{-i\left(\omega_{b}+\Omega_{1}\right)t}\right)+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).

###  A.2 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

$\displaystyle\hat{H}$
$\displaystyle=\omega_{a}\hat{a}^{\dagger}\hat{a}+\sum_{j=1}^{2}\left(\frac{\Omega_{j}}{2}\hat{\sigma}^{(j)}_{z}+g(t)\left(\hat{a}+\hat{a}^{\dagger}\right)\hat{\sigma}^{(j)}_{x}\right),$
(42)

with the modulated coupling

$\displaystyle
g(t)=\bar{g}\left(\alpha_{-}e^{-i\left(\omega_{a}-\bar{\Omega}\right)t}+\alpha_{+}e^{-i\left(\omega_{a}+\bar{\Omega}\right)t}\right)+h.c.,$

where $\bar{\Omega}=(\Omega_{1}+\Omega_{2})/2$ is the average qubit frequency.

We then move to the interaction frame to obtain the Hamiltonian

$\displaystyle\hat{H}^{\prime}=\bar{g}\hat{a}^{\dagger}\Big{[}$ $\displaystyle
e^{i\delta\Omega
t}\left(\alpha_{+}\hat{\sigma}_{+}^{(1)}+\alpha_{-}\hat{\sigma}_{-}^{(2)}\right)+e^{-i\delta\Omega
t}\left(\alpha_{+}\hat{\sigma}_{+}^{(2)}+\alpha_{-}\hat{\sigma}_{-}^{(1)}\right)$
$\displaystyle+e^{i\Sigma_{1}t}\alpha_{-}\hat{\sigma}_{+}^{(1)}+e^{-i\Sigma_{1}t}\alpha_{+}\hat{\sigma}_{-}^{(1)}$
$\displaystyle+e^{i\Sigma_{2}t}\alpha_{-}\hat{\sigma}_{+}^{(2)}+e^{-i\Sigma_{2}t}\alpha_{+}\hat{\sigma}_{-}^{(2)}$
$\displaystyle+e^{i(2\omega_{a}+\delta\Omega)t}\left(\alpha^{*}_{+}\hat{\sigma}_{-}^{(2)}+\alpha^{*}_{-}\hat{\sigma}_{+}^{(1)}\right)$
$\displaystyle+e^{i(2\omega_{a}-\delta\Omega)t}\left(\alpha^{*}_{+}\hat{\sigma}_{-}^{(1)}+\alpha_{-}^{*}\hat{\sigma}_{+}^{(2)}\right)$
$\displaystyle+e^{i(2\omega_{a}+\Sigma_{1})t}\alpha^{*}_{+}\hat{\sigma}_{+}^{(1)}+e^{i(2\omega_{a}+\Sigma_{2})t}\alpha^{*}_{+}\hat{\sigma}_{+}^{(2)}$
$\displaystyle\left.+e^{i(2\omega_{a}-\Sigma_{1})t}\alpha^{*}_{-}\hat{\sigma}_{-}^{(1)}+e^{i(2\omega_{a}-\Sigma_{2})t}\alpha^{*}_{-}\hat{\sigma}_{-}^{(2)}\right]+h.c.,$
(43)

where $\delta\Omega=(\Omega_{1}-\Omega_{2})/2$ and
$\Sigma_{1/2}=(3\Omega_{1/2}+\Omega_{2/1})/2$. If $\bar{g}>\delta\Omega$ we
cannot ignore terms that oscillate at $\pm\delta\Omega$, 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\omega_{a}-\Sigma_{1/2}$. For these to be fast oscillating, we
require that $2\omega_{a}-\Sigma_{1/2}\gg\bar{g}$. This simply implies that
when designing the system, we must be careful to avoid the “resonance”
conditions $2\omega_{a}=\Sigma_{1/2}$, so as not to activate unwanted qubit
evolution.

Having dropped all fast oscillating terms the Hamiltonian bcomes

$\displaystyle\hat{H}^{\prime}=\bar{g}\hat{a}^{\dagger}\left[e^{i\delta\Omega
t}\left(\alpha_{+}\hat{\sigma}_{+}^{(1)}+\alpha_{-}\hat{\sigma}_{-}^{(2)}\right)\right.$
$\displaystyle\left.+e^{-i\delta\Omega
t}\left(\alpha_{+}\hat{\sigma}_{+}^{(2)}+\alpha_{-}\hat{\sigma}_{-}^{(1)}\right)\right]+h.c..$
(44)

Note that $\bar{g}>\delta\Omega$ is the opposite operating condition to the
two cavity setup, where we require $\delta\Omega>\bar{g}$.

Assuming the cavity has a sufficiently broad linewidth, i.e.
$\delta\Omega\ll\kappa$, then the qubit processes in Eq. (44) will interact
with the environment through the cavity’s nonzero density of states at
$\pm\delta\Omega$, as depicted in Fig. 3c). If $2\delta\Omega$ is greater than
the effective qubit linewidths, which are at least
$4\alpha^{2}\bar{g}^{2}/\kappa$, then the processes at $+\delta\Omega$ and
$-\delta\Omega$ 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).

###  A.3 Two qubits and a transmission line

We consider two qubits coupled to a long waveguide or transmission line (TL)
described by the Hamiltonian

$\hat{H}=\sum_{\alpha=1}^{2}\frac{\Omega_{\alpha}}{2}\hat{\sigma}_{z}^{(\alpha)}+\hat{H}_{TL}+\hat{H}_{I}.$
(45)

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 $\lambda$)

$\hat{H}_{TL}=\sum_{k\lambda}\omega_{k\lambda}\hat{b}_{k\lambda}^{\dagger}\hat{b}_{k\lambda},$
(46)

where we assume the dispersion relations satisfy
$\omega_{k,\lambda}=\omega_{-k,\lambda}$. The qubit-TL couplings are described
by

$\hat{H}_{I}=\sum_{\alpha=1}^{2}g_{\alpha}(t)\hat{\sigma}_{x}^{(\alpha)}\hat{B}_{\alpha},$
(47)

where

$\hat{B}_{\alpha}=\sum_{k\lambda}\left(\frac{1}{\sqrt{L}}e^{ikx_{\alpha}}f_{k\lambda}\hat{b}_{k\lambda}^{\dagger}+\mathrm{h.c.}\right),$
(48)

with $x_{\alpha}$ denoting the position of the $\alpha$ qubit along the TL.
$f_{k\lambda}$ is a wavevector and band dependent coupling coefficient which
we assume to satisfy $f_{k\lambda}=f_{-k\lambda}^{*}$. To realize the
synthetic squeezing dissipator, the qubit-TL couplings $g_{\alpha}(t)$ are
modulated at the sum of the qubit splitting frequencies, and also have a
static component:

$g_{\alpha}(t)=\overline{g}\left(\eta_{\alpha}+2\epsilon_{\alpha}\cos\left[\left(\Omega_{1}+\Omega_{2}\right)t\right]\right).$
(49)

Note that in this description $\overline{g}$ has units
$\mathrm{Hz}\sqrt{\mathrm{m}}$.

After moving into the interaction picture with respect to the system
Hamiltonian
$\hat{H}_{S}=\sum_{\alpha=1}^{2}(\Omega_{\alpha}/2)\hat{\sigma}_{z}^{(\alpha)}$
and the TL Hamiltonian $\hat{H}_{TL}$, the system-TL interaction can be
written as

$\displaystyle\hat{H}_{I,\mathrm{int}}/\overline{g}$
$\displaystyle=\left[\eta_{1}e^{i\Omega_{1}t}+\epsilon_{1}e^{-i\Omega_{2}t}+\epsilon_{1}e^{i\left(2\Omega_{1}+\Omega_{2}\right)t}\right]\hat{\sigma}_{+}^{(1)}\hat{B}_{1}(t)$
$\displaystyle+\,\left[\eta_{2}e^{i\Omega_{2}t}+\epsilon_{2}e^{-i\Omega_{1}t}+\epsilon_{2}e^{i\left(\Omega_{1}+2\Omega_{2}\right)t}\right]\hat{\sigma}_{+}^{(2)}\hat{B}_{2}(t)$
$\displaystyle+\,\mathrm{h.c.}$ (50)

which is of the form

$\hat{H}_{I,\mathrm{int}}=\sum_{\alpha,\omega}e^{-i\omega
t}\hat{A}_{\alpha}(\omega)\hat{B}_{\alpha}(t)$ (51)

with $\alpha\in\\{1,2\\}$ indexing the qubits and
$\omega\in\\{\pm\Omega_{1},\pm\Omega_{2},\pm(2\Omega_{1}+\Omega_{2}),\pm(\Omega_{1}+2\Omega_{2})\\}$
indexing all of the relevant frequencies. The $\hat{A}_{\alpha}(\omega)$ are
e.g., $\hat{A}_{1}(\Omega_{1})=\overline{g}\eta_{1}\hat{\sigma}_{-}^{(1)}$ and
$\hat{A}_{2}(\Omega_{1})=\overline{g}\epsilon_{1}\hat{\sigma}_{+}^{(2)}$. With
$\hat{H}_{I,\mathrm{int}}$ written in this form and following the master
equation derivation of Breuer and Petruccione (2002), we find the system
density matrix master equation to be

$\displaystyle\frac{d}{dt}\hat{\rho}(t)=\sum_{\alpha,\beta}\sum_{\omega}\Gamma_{\alpha\beta}(\omega)\left(\hat{A}_{\beta}(\omega)\hat{\rho}(t)\hat{A}_{\alpha}^{\dagger}(\omega)\right.$
$\displaystyle\left.-\>\hat{A}_{\alpha}^{\dagger}(\omega)\hat{A}_{\beta}(\omega)\hat{\rho}(t)\right)+\mathrm{h.c.}$
(52)

where we have introduced the bath correlation functions

$\Gamma_{\alpha\beta}(\omega)\equiv\int_{0}^{\infty}dse^{i\omega
s}\langle\hat{B}_{\alpha}^{\dagger}(s)\hat{B}_{\beta}(0)\rangle.$ (53)

Here
$\langle\hat{B}_{\alpha}^{\dagger}(s)\hat{B}_{\beta}(0)\rangle=\mathrm{Tr}_{TL}\\{\hat{B}_{\alpha}^{\dagger}(s)\hat{B}_{\beta}(0)\hat{\rho}_{TL}\\}$
and we have assumed the TL state $\hat{\rho}_{TL}$ is stationary so
$\langle\hat{B}_{\alpha}^{\dagger}(t)\hat{B}_{\beta}(t-s)\rangle=\langle\hat{B}_{\alpha}^{\dagger}(s)\hat{B}_{\beta}(0)\rangle$.

We next decompose the bath correlation functions as

$\Gamma_{\alpha\beta}(\omega)=\frac{1}{2}\gamma_{\alpha\beta}(\omega)+iS_{\alpha\beta}(\omega),$
(54)

for
$S_{\alpha\beta}(\omega)=(\Gamma_{\alpha\beta}(\omega)-\Gamma_{\beta\alpha}^{*}(\omega))/2i$
and
$\gamma_{\alpha\beta}(\omega)=\Gamma_{\alpha\beta}(\omega)+\Gamma_{\beta\alpha}^{*}(\omega)$.
This decomposition separates the master equation into dissipative terms with
correlation function coefficients $\gamma_{\alpha\beta}(\omega)$ and
Hamiltonian terms with correlation function coefficients
$S_{\alpha\beta}(\omega)$. The master equation can then be recast into the
form

$\frac{d}{dt}\hat{\rho}(t)=-i[\hat{H}_{LS},\hat{\rho}(t)]+D(\hat{\rho}(t)),$
(55)

where the coherent evolution of the system due to the transmission line is
given by the “Lamb shift” Hamiltonian

$\hat{H}_{LS}=\sum_{\omega}\sum_{\alpha,\beta}S_{\alpha\beta}(\omega)\hat{A}_{\alpha}^{\dagger}(\omega)\hat{A}_{\beta}(\omega),$
(56)

and the dissipation is given by

$\displaystyle
D(\hat{\rho}(t))=\sum_{\omega}\sum_{\alpha,\beta}\gamma_{\alpha\beta}(\omega)\left(\hat{A}_{\beta}(\omega)\hat{\rho}(t)\hat{A}_{\alpha}^{\dagger}(\omega)\right.$
$\displaystyle\left.-\>\frac{1}{2}\\{\hat{A}_{\alpha}^{\dagger}(\omega)\hat{A}_{\beta}(\omega),\hat{\rho}(t)\\}\right).$
(57)

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

$\hat{H}_{\mathrm{qb-
qb}}=\Lambda\hat{\sigma}_{+}^{(1)}\hat{\sigma}_{+}^{(2)}+\mathrm{h.c.}$ (58)

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 $\gamma_{\alpha\beta}(\omega)$ for all possible
frequencies
$\omega\in\\{\pm\Omega_{1},\pm\Omega_{2},\pm(2\Omega_{1}+\Omega_{2}),\pm(\Omega_{1}+2\Omega_{2})\\}$.
We let the TL be at zero temperature so that

$\gamma_{\alpha\beta}(\omega<0)=0.$ (59)

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.2 we investigate the impact of a
thermal environment on entanglement stabilization.

The dissipative processes at frequency $2\Omega_{1}+\Omega_{2}$
($\Omega_{1}+2\Omega_{2}$) involve both single qubit loss on qubit 1 (2), due
to the static couplings $\eta_{i}$, as well as single qubit excitation on
qubit 2 (1), due to the modulated couplings $\epsilon_{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

$\gamma_{\alpha\beta}(2\Omega_{1}+\Omega_{2})=\gamma_{\alpha\beta}(\Omega_{1}+2\Omega_{2})=0.$
(60)

We are thus left with the dissipative processes at $\Omega_{1}$ and
$\Omega_{2}$. The $\Omega_{1}$ processes take the form

$\displaystyle D_{\Omega_{1}}(\hat{\rho}(t))/\overline{g}^{2}$
$\displaystyle=\eta_{1}^{2}\gamma_{11}(\Omega_{1})\mathcal{D}[\hat{\sigma}_{-}^{(1)}]\hat{\rho}(t)$
$\displaystyle+\,\eta_{1}\epsilon_{2}\gamma_{12}(\Omega_{1})\mathcal{D}[\hat{\sigma}_{+}^{(2)},\hat{\sigma}_{-}^{(1)}]\hat{\rho}(t)$
$\displaystyle+\,\eta_{1}\epsilon_{2}\gamma_{21}(\Omega_{1})\mathcal{D}[\hat{\sigma}_{-}^{(1)},\hat{\sigma}_{+}^{(2)}]\hat{\rho}(t)$
$\displaystyle+\,\epsilon_{2}^{2}\gamma_{22}(\Omega_{1})\mathcal{D}[\hat{\sigma}_{+}^{(2)}]\hat{\rho}(t),$
(61)

where $\mathcal{D}[\hat{O}]\hat{\rho}$ is the usual Lindblad dissipator and
$\mathcal{D}[\hat{O},\hat{P}]\hat{\rho}=\hat{O}\hat{\rho}\hat{P}^{\dagger}-\\{\hat{P}^{\dagger}\hat{O},\hat{\rho}\\}/2$.
The $\Omega_{2}$ process is the same but with the indices $1\leftrightarrow
2$. Note that if the correlation functions
$\gamma_{\alpha\beta}(\Omega_{1})=\gamma_{1}$ were all equal, this dissipator
would reduce to the collective dissipator

$D_{\Omega_{1}}(\hat{\rho}(t))=\overline{g}^{2}\gamma_{1}\mathcal{D}[\eta_{1}\hat{\sigma}_{-}^{(1)}+\epsilon_{2}\hat{\sigma}_{+}^{(2)}]\hat{\rho}(t).$
(62)

The above correlation functions are readily evaluated in general (in the
continuum limit $\frac{1}{L}\sum_{k}\to\int\frac{dk}{2\pi}$) as

$\gamma_{\alpha\beta}(\omega)=\sum_{\lambda}\sum_{j}\frac{|f_{\lambda}(k_{j})|^{2}}{v_{g,\lambda}(k_{j})}e^{ik_{j}(x_{\beta}-x_{\alpha})},$
(63)

where the $\\{k_{j}\\}$ are all wavevectors for which
$\omega_{\lambda}(k)=\omega$ and
$v_{g,\lambda}(k)=|\partial\omega_{\lambda}(k)/\partial k|$ is the group
velocity at frequency $\omega_{\lambda}(k)$. Then for the $\Omega_{1}$
processes, we have

$\gamma_{\alpha\beta}(\Omega_{1})=\frac{2}{v_{g}(k_{1})}\cos\left(k_{1}(x_{\beta}-x_{\alpha})\right),$
(64)

where we have dropped the band index $\lambda$ and absorbed the coupling
coefficient $f_{\lambda}(k)$ into the group velocity. Here $k_{1}$ is the
wavevector for which $\omega_{\lambda}(k_{1})=\Omega_{1}$ for the correct band
$\lambda$. We get a similar result for the $\Omega_{2}$ processes.

Defining the qubit spacing

$l\equiv x_{2}-x_{1},$ (65)

we can rewrite the $\Omega_{1}$ dissipation as

$\displaystyle
D_{\Omega_{1}}(\hat{\rho}(t))=\frac{2\overline{g}^{2}}{v_{g}(k_{1})}\cos(k_{1}l)\mathcal{D}[\eta_{1}\hat{\sigma}_{-}^{(1)}+\epsilon_{2}\hat{\sigma}_{+}^{(2)}]\hat{\rho}(t)$
$\displaystyle+\>\frac{2\overline{g}^{2}}{v_{g}(k_{1})}\left(1-\cos(k_{1}l)\right)\left(\eta_{1}^{2}\mathcal{D}[\hat{\sigma}_{-}^{(1)}]+\epsilon_{2}^{2}\mathcal{D}[\hat{\sigma}_{+}^{(2)}]\right)\hat{\rho}(t).$
(66)

In the ideal case where $k_{1}l=\pi n$ for $n\in\mathbb{Z}$, this all reduces
to a single dissipator with collective jump operator

$\hat{\mathcal{J}}_{1}=\eta_{1}\hat{\sigma}_{-}^{(1)}\pm\epsilon_{2}\hat{\sigma}_{+}^{(2)},$
(67)

where the sign is set by the parity of $n$ ($+$ for even $n$). In the more
general case where $k_{1}l\neq\pi n$, Eq. (66) 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=\pi n+\delta$ where $\delta\ll 1$, these
additional dissipators are given to leading order in $\delta$ by

$\displaystyle D_{\Omega_{1},{\rm
bad}}(\hat{\rho}(t))=\delta^{2}\frac{2\bar{g}^{2}}{v_{g}(k_{1})}\left(\eta_{1}^{2}\mathcal{D}[\hat{\sigma}_{-}^{(1)}]+\epsilon_{2}^{2}\mathcal{D}[\hat{\sigma}_{+}^{(2)}]\right)\hat{\rho}(t).$
(68)

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
$\Omega_{2}$ in the TL. The dissipation reduces to a single collective jump
operator

$\hat{\mathcal{J}}_{2}=\epsilon_{1}\hat{\sigma}_{+}^{(1)}\pm\eta_{2}\hat{\sigma}_{-}^{(2)},$
(69)

when $k_{2}l=\pi m$ for $m\in\mathbb{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=\pi n$ and $k_{2}l=\pi m$ for
some spacing $l$ and integers $n$ and $m$, the final master equation is

$\frac{d}{dt}\hat{\rho}(t)=-i[\hat{H}_{LS},\hat{\rho}]+\Gamma_{1}\mathcal{D}[\hat{\mathcal{J}}_{1}]\hat{\rho}+\Gamma_{2}\mathcal{D}[\hat{\mathcal{J}}_{2}]\hat{\rho}.$
(70)

Identifying $\eta_{1}=\alpha_{-}$, $\eta_{2}=\beta_{-}$,
$\epsilon_{1}=\beta_{+}$, and $\epsilon_{2}=\alpha_{+}$, we arrive at the
desired master equation Eq. (8) _except_ with the additional Lamb shift
Hamiltonian. Note that if the qubit spacing is such that $k_{1}l=\pi n$ and
$k_{2}l=\pi m$ for $m$, $n$ of opposite parity, we must flip the modulation
phase for one of the qubit-waveguide couplings Eq. (49), e.g.,
$\epsilon_{2}\mapsto-\epsilon_{2}$. The signs in the dissipators
$\hat{\mathcal{J}}_{1}$, $\hat{\mathcal{J}}_{2}$ must be the same to stabilize
the desired entangled state.

We now consider the Lamb shift Hamiltonian $\hat{H}_{LS}$. It takes the
general form

$\hat{H}_{LS}=\sum_{\alpha=1}^{2}\Omega_{\alpha,LS}\hat{\sigma}_{z}^{(\alpha)}+\left(\Lambda\hat{\sigma}_{+}^{(1)}\hat{\sigma}_{+}^{(2)}+\mathrm{h.c.}\right),$
(71)

where $\Omega_{\alpha,LS}$ are the inconsequential single qubit frequency
shifts, and the pairing interaction amplitude is

$\displaystyle\Lambda/\overline{g}^{2}=\eta_{1}\epsilon_{2}\left(S_{12}(\Omega_{1})+S_{21}(-\Omega_{1})\right)$
$\displaystyle+\>\epsilon_{1}\eta_{2}\left(S_{21}(\Omega_{2})+S_{12}(-\Omega_{2})\right).$
(72)

We show below that $\Lambda$ is real, thus the eigenstates of the pairing
interaction are paired states:

$\ket{\psi_{\pm}}=\frac{1}{\sqrt{2}}\left(\ket{00}\pm\ket{11}\right).$ (73)

Recall that synthetic squeezing generates the paired state given by Eq. (3),
where the relative phase between $\ket{00}$ and $\ket{11}$ is $\pm 1$. Thus,
one of the eigenstates of the pairing Hamiltonian is always the $r\to\infty$
limit of the synthetic squeezing steady state, so that as $r$ is made larger
the Hamiltonian will induce less disruption to our entanglement scheme. For
any $r<\infty$ it will reduce the purity of the steady state, so we now
evaluate $\Lambda$ to determine the magnitude of this impact.

Unlike the dissipation correlation functions $\gamma_{\alpha\beta}(\omega)$
which are readily evaluated with minimal specification of the TL band
structure, the $S_{\alpha\beta}(\omega)$ cannot be evaluated without
specifying the band structure of the TL or the $k$-dependent coefficients
$f_{k\lambda}$ of the system-TL coupling. Nevertheless, we can make some
general statements about $S_{\alpha\beta}(\omega)$, and thus about $\Lambda$,
given the assumptions we have already made about the TL. The
$S_{\alpha\beta}(\omega)$ are given by

$S_{\alpha\beta}(\omega)=\sum_{\lambda}\,{\rm
pv}\\!\int_{-k_{c}}^{k_{c}}\frac{dk}{2\pi}\frac{e^{ik(x_{\beta}-x_{\alpha})}}{\omega-\omega_{\lambda}(k)}|f_{\lambda}(k)|^{2},$
(74)

where pv denotes the principal value and $k_{c}$ is the band edge cutoff. At
the outset, we had assumed reciprocity of the TL which implies
$\omega_{\lambda}(k)=\omega_{\lambda}(-k)$ and
$f_{\lambda}(k)=f_{\lambda}(-k)$. Thus the integral can be rewritten as

$S_{\alpha\beta}(\omega)=2\sum_{\lambda}\,{\rm
pv}\\!\int_{0}^{k_{c}}\frac{dk}{2\pi}\frac{\cos(kl_{\alpha\beta})}{\omega-\omega_{\lambda}(k)}|f_{\lambda}(k)|^{2},$
(75)

where $l_{\alpha\beta}\equiv x_{\beta}-x_{\alpha}$. Therefore we see that
$S_{\alpha\beta}(\omega)$ are real, thus so is $\Lambda$.

We can say little more about the correlation functions
$S_{\alpha\alpha}(\omega)$ in general without more details of the TL, as these
correlation functions are highly dependent on the band structure of the TL.
Fortunately the pairing amplitude $\Lambda$ contains only
$S_{\alpha\neq\beta}(\omega)$ factors which have the oscillating term
$\cos(kl)$ in their integrals. If we let the qubit spacing be large,
$k_{j}l\gg 1$, then the support of the integrals will be in a narrow region
around the pole $\omega-\omega_{\lambda}(k)$. This implies a few things:
first, the negative frequency correlation functions are negligible
$|S_{\alpha\beta}(-\Omega_{j})|\ll|S_{\beta\alpha}(+\Omega_{j})|$ since they
never encounter a pole. Second, the contributions to the positive frequency
correlation functions by all bands other than that which contains $\Omega_{j}$
are similarly negligible. Finally, even in the band that contains $\Omega_{j}$
the integral is not particularly sensitive to either the band edge nor the
dispersion relation since the rapid oscillation and $1/k$ behavior drastically
reduce the integral’s support away from the pole.

With the above general observations, we evaluate a rough estimate of $\Lambda$
for a typical TL. Suppose the TL has a low frequency band ($\lambda=1)$ for
which $\omega(k_{j})=\Omega_{j}$ (dropping the band index). For concreteness
let’s consider $S_{12}(\Omega_{1})$. We linearize the dispersion relation
about the pole $\omega(k)=\Omega_{1}$ as
$\omega(k)\approx\omega(k_{1})+v_{g}(k_{1})(k-k_{1})$ (for a typical TL far
from the cutoff, this is a reasonable approximation). Then
$S_{12}(\Omega_{1})$ is

$S_{12}(\Omega_{1})\approx\frac{2}{v_{g}(k_{1})}\,{\rm
pv}\\!\int_{0}^{k_{c}}\frac{dk}{2\pi}\frac{\cos(kl)}{k_{1}-k}|f(k)|^{2}.$ (76)

Finally, we assume that $f(k)\approx{\rm const}$ in the region around $k_{1}$
and absorb it into the group velocity. Since the integrand oscillates rapidly
and quickly dies off away from $k=k_{1}$, we extend the integration bounds to
$\pm\infty$, yielding

$S_{12}(\Omega_{1})\approx\frac{1}{v_{g}(k_{1})}\sin(k_{1}l).$ (77)

A similar result holds for $S_{21}(\Omega_{2})$. Therefore, we find the
pairing amplitude $\Lambda$ to be

$\Lambda/\bar{g}^{2}\approx\frac{\eta_{1}\epsilon_{2}}{v_{g}(k_{1})}\sin(k_{1}l)+\frac{\eta_{2}\epsilon_{1}}{v_{g}(k_{2})}\sin(k_{2}l),$
(78)

which vanishes for $k_{j}l=\pi n_{j}$ (integers $n_{j}$), exactly the
conditions needed to realize synthetic squeezing.

In general, $\Lambda$ won’t vanish exactly at the synthetic squeezing
conditions. However, numerical evaluation of $\Lambda$ for various $\omega(k)$
(examples considered include $\omega=v_{g}|k|$, $\omega=(v_{g}/a)|\sin ka|$,
and $\omega=v_{g}(|k|-\alpha k^{2})$) and for moderate qubit spacing
$k_{j}l\sim 10\pi$ shows that, for spacing near the synthetic squeezing
conditions, typically

$\Lambda/\bar{g}^{2}\approx\frac{\eta_{1}\epsilon_{2}}{v_{g}(k_{1})}(\sin(k_{1}l)+\delta_{1})+\frac{\eta_{2}\epsilon_{1}}{v_{g}(k_{2})}(\sin(k_{2}l)+\delta_{2})$
(79)

where usually $|\delta_{j}|\simeq 10^{-3}$. Thus $\Lambda$ 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_{\rm ideal}+\delta l$, we consider a simple model system
and analyze its performance numerically. To be concrete, consider two
superconducting qubits with frequencies $\Omega_{1}=4$ GHz and $\Omega_{2}=6$
GHz and a transmission line with a linear dispersion

$\omega(k)=v_{g}|k|,$ (80)

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\Omega_{1}+\Omega_{2}=14$ GHz and $\Omega_{1}+2\Omega_{2}=16$ GHz).
We choose to space the qubits by a distance

$l=4\pi/k_{1}+\delta l,$ (81)

or approximately two wavelengths at $\Omega_{1}$ (three wavelengths at
$\Omega_{2}$). The ideal spacing, $\delta l=0$, satisfies the synthetic
squeezing conditions $k_{1}l=4\pi$ and $k_{2}l=6\pi$. We let
$\bar{g}^{2}/v_{g}\equiv 1$ and set $\eta_{1}=\eta_{2}=\cosh(r)$ and
$\epsilon_{1}=\epsilon_{2}=\sinh(r)$. Then to a very good approximation

$\Lambda=\sinh(r)\cosh(r)\left(\sin(k_{1}\delta l)+\sin(k_{2}\delta
l)\right).$ (82)

Thus the master equation is given by

$\displaystyle\dot{\hat{\rho}}=-i[\hat{H},\hat{\rho}]+\sum_{j}2\cos(k_{j}\delta
l)\mathcal{D}[\hat{\mathcal{J}}_{j}]\hat{\rho}$ (83)
$\displaystyle+4\sin^{2}\left(\frac{k_{1}\delta
1}{2}\right)\left(\cosh^{2}(r)\mathcal{D}[\hat{\sigma}_{-}^{(1)}]+\sinh^{2}(r)\mathcal{D}[\hat{\sigma}_{+}^{(2)}]\right)\hat{\rho}$
$\displaystyle+4\sin^{2}\left(\frac{k_{2}\delta
l}{2}\right)\left(\cosh^{2}(r)\mathcal{D}[\hat{\sigma}_{-}^{(2)}]+\sinh^{2}(r)\mathcal{D}[\hat{\sigma}_{+}^{(1)}]\right)\hat{\rho},$

with $\hat{H}=\Lambda(\hat{\sigma}_{-}^{(1)}\hat{\sigma}_{-}^{(2)}+{\rm
h.c.})$ and
$\hat{\mathcal{J}}_{1/2}=\cosh(r)\hat{\sigma}_{-}^{(1/2)}+\sinh(r)\hat{\sigma}_{+}^{(2/1)}$.

We numerically compute the steady state of this master equation as a function
of $\delta l$ and $r$. In general for $\delta l\neq 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 $\delta 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 $|\delta l|$ as a fraction of the wavelength $\lambda_{1}=2\pi/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.

Figure 5: Performance of transmission line scheme with qubit spacing errors.
The purity (red) and optimal concurrence (black) of the steady state of Eq.
(83) as a function of qubit spacing error $\delta l$ (measured as a fraction
of the wavelength $\lambda_{1}$). The concurrence is optimized over the
squeezing parameter $r$, and the optimal squeezing (blue) is plotted against
the right axis. The optimal squeezing diverges as $\delta l\to 0$. The solid
curves show the results of the full master equation and the dashed curves show
the results of the master equation _without_ the Hamiltonian. The degradation
of the entanglement is driven mostly by the single qubit dissipators.

##  Appendix B Collective loss plus local qubit driving scheme: additional
details

###  B.1 Asymmetric waveguide couplings

We consider the generalization of the scheme presented in Sec. IV to the
general collective loss dissipator with arbitrary $\eta>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 $\eta\neq 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

$\hat{H}=\frac{\Delta+\epsilon}{2}\hat{\sigma}_{z}^{(1)}-\frac{\Delta-\epsilon}{2}\hat{\sigma}_{z}^{(2)}+\Lambda^{(1)}\hat{\sigma}_{x}^{(1)}+\Lambda^{(2)}\hat{\sigma}_{x}^{(2)},$
(84)

wherein we allow for asymmetric Rabi drive amplitudes $\Lambda^{(j)}$ and
asymmetric detunings parameterized by $\epsilon$.

The dark states of the collective of loss Eq. (13) can be expressed as

$\ket{\phi}=\alpha\ket{00}+\beta\left(\ket{01}-\eta\ket{10}\right),$ (85)

for $|\alpha|^{2}+(1+\eta^{2})|\beta|^{2}=1$. We take $\alpha,\beta>0$ real
and positive without loss of generality. Enforcing that $\ket{\phi}$ is an
eigenstate of $\hat{H}$, i.e. $\hat{H}\ket{\phi}=E\ket{\phi}$, imposes the
following constraints on the Hamiltonian parameters

$\displaystyle\Lambda^{(1)}=\eta\Lambda^{(2)},$ (86)
$\displaystyle\epsilon=\frac{\Lambda^{2}}{\Delta}\left(1-\eta^{2}\right),$
(87)

where in the latter we have defined $\Lambda\equiv\Lambda^{(2)}$. Imposing
these conditions yields Eq. (29) as the final Hamiltonian.

The unique steady state expressed in terms of the Hamiltonian parameters is
given by

$\ket{\Phi}=\frac{1}{\sqrt{\Delta^{2}+\Lambda^{2}(1+\eta^{2})}}\left(\Delta\ket{00}+\Lambda\left(\ket{01}-\eta\ket{10}\right)\right).$
(88)

For symmetric loss, there is a frame in which this state is an entangled state
of the form of Eq. (3), and we now show this is also the case for asymmetric
loss. As with the symmetric case, the correct frame will be defined by the
eigenbasis of the Hamiltonian of Eq. (29). This Hamiltonian is diagonalized by
rotating about each qubit’s $y$-axis through the angles

$\displaystyle\theta_{1}$
$\displaystyle=\arctan\left(\frac{2\eta\Lambda}{\Delta+\epsilon}\right),$ (89)
$\displaystyle\theta_{2}$
$\displaystyle=\arctan\left(\frac{-2\Lambda}{\Delta-\epsilon}\right),$ (90)

with these angles lying on the principal branch of the arctangent.

In the new frame, Eq. (29) transforms to

$\hat{H}^{\prime}=\frac{\mu}{2}\left(\hat{\tau}_{z}^{(1)}-\hat{\tau}_{z}^{(2)}\right),$
(91)

for

$\mu=\sqrt{(\Delta-\epsilon)^{2}+4\Lambda^{2}}.$ (92)

Making this transformation on the state $\ket{\Phi}$ yields our target state

$\ket{\Phi}^{\prime}=\frac{1}{\sqrt{\beta_{-}^{2}+\beta_{+}^{2}}}\left(\beta_{-}\ket{00}^{\prime}-\beta_{+}\ket{11}^{\prime}\right),$
(93)

where $\beta_{-}>\beta_{+}$, and they are given by

$\displaystyle\beta_{-}$ $\displaystyle=\mu+\Delta+\epsilon,$ (94)
$\displaystyle\beta_{+}$ $\displaystyle=\eta(\mu-\Delta+\epsilon).$ (95)

$\ket{\Phi}^{\prime}$ can be made to look explicitly like a squeezed state by
introducing the squeezing parameter $r={\rm artanh}\,\beta_{+}/\beta_{-}$. For
the symmetric case ($\eta=1$) this expression reduces to
$r=(1/4)\ln(1+4\Lambda^{2}/\Delta^{2})$ which is a rewriting of Eq. (27). In
the general case, one cannot easily find the expressions for $\Delta$ and
$\Lambda$ in terms of $r$, $\mu$, and $\eta$. Instead one must solve the
following equations for $\Delta$ and $\Lambda$ in terms of the other
parameters (recalling that $\epsilon$ is defined in terms of $\Delta$ and
$\Lambda$, c.f. Eq. (87))

$\displaystyle\Delta$
$\displaystyle=-\frac{e^{2r}-1-\eta\left(1+e^{2r}\right)}{e^{2r}-1+\eta\left(1+e^{2r}\right)}\left(\mu+\epsilon\right),$
(96) $\displaystyle\Lambda$
$\displaystyle=\frac{1}{2}\sqrt{\mu^{2}-(\Delta-\epsilon)^{2}}.$ (97)

Finally in the Rabi eigenbasis, the collective loss dissipator becomes

$\displaystyle\hat{\mathcal{J}}^{\prime}$
$\displaystyle=\hat{\mathcal{J}}^{\prime}_{Z}+\hat{\mathcal{J}}^{\prime}_{1}+\hat{\mathcal{J}}^{\prime}_{2},$
(98) $\displaystyle\hat{\mathcal{J}}^{\prime}_{Z}$
$\displaystyle=\frac{\eta\Lambda}{\mu}\left(\hat{\tau}_{z}^{(1)}-\hat{\tau}_{z}^{(2)}\right),$
(99) $\displaystyle\hat{\mathcal{J}}^{\prime}_{1}$
$\displaystyle=\frac{1}{2\mu}\left(\beta_{-}\hat{\tau}_{-}^{(1)}-\beta_{+}\hat{\tau}_{+}^{(2)}\right),$
(100) $\displaystyle\hat{\mathcal{J}}^{\prime}_{2}$
$\displaystyle=\frac{(\mu-\epsilon)^{2}-\Delta^{2}}{8\eta\mu\Lambda^{2}}\left(\beta_{-}\hat{\tau}_{-}^{(2)}-\beta_{+}\hat{\tau}_{+}^{(1)}\right).$
(101)

This has precisely the same operator form as $\hat{\mathcal{J}}^{\prime}$
given in the main text for the symmetric case. One difference that appears in
the asymmetric case is that the coefficients of
$\hat{\mathcal{J}}^{\prime}_{1}$ and $\hat{\mathcal{J}}^{\prime}_{2}$ differ.
Only when $\eta=1$ do the coefficients coincide. This has no effect on the
steady state as the dissipation rates do not change the form of the steady
state.

###  B.2 Comparing dynamics against the ideal two-qubit squeezed dissipation
master equation

In Sec. IV we showed that with collective loss and local Rabi drives on two
qubits, there is a frame in which the steady state is a two-qubit entangled
state and that the master equation in this frame is given by

$\dot{\hat{\rho}}^{\prime}=-i[\hat{H}^{\prime},\hat{\rho}^{\prime}]+\Gamma\mathcal{D}\left[\hat{\mathcal{J}}^{\prime}_{1}+\hat{\mathcal{J}}^{\prime}_{2}+\hat{\mathcal{J}}^{\prime}_{Z}\right]\hat{\rho}^{\prime},$
(102)

with the Hamiltonian given by Eq. (22) and the dissipator terms given by Eqs.
(19)-(21). We further showed that when the Hamiltonian energy $\mu\gg\Gamma$
(the collective loss rate), the three terms in the dissipator can be split
into three separate dissipators, as cross-terms between them become fast-
oscillating in the interaction frame defined by $\hat{H}^{\prime}$.

The resulting master equation, Eq. (24), is nearly identical to the dynamics
of the ideal two-qubit squeezing master equation, Eq. (2) (with
$\gamma_{1}=\gamma_{2}=\Gamma$). Moving to the interaction frame of
$\hat{H}^{\prime}$ in Eq. (24) does not affect the dissipators, so that the
only remaining difference between the two master equations is the addition of
the correlated dephasing $\hat{\mathcal{J}}_{Z}^{\prime}$ to Eq. (24). We now
show that this additional dissipation has virtually no qualitative or
quantitative effect on the stabilization of the entangled two-qubit state, and
thus that Eq. (102) and Eq. (2) are essentially the same when $\mu\gg\Gamma$.

As we discussed in Sec. V, Eq. (2) has a dissipative gap $\kappa_{\rm slow}$
which is exponentially small in the squeezing 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. (102) as a
function of $\mu/\Gamma$ 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
$\mu=0$ (no Hamiltonian) because the collective loss dissipator has multiple
steady states.

A simple perturbative argument suggests that the gap should open as
$\kappa_{\rm slow}\sim\mu^{2}$ which is precisely what is observed for
$\mu\ll\Gamma$. When $\mu\gg\Gamma$ 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 $\hat{\mathcal{J}}^{\prime}_{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.

Figure 6: Comparison of dynamics between collective loss plus local driving
and ideal two-qubit dissipative squeezing. The dissipative gap $\kappa_{\rm
slow}$ of Eq. (102) as a function of the Hamiltonian energy $\mu$, computed
numerically for various squeezing parameters $r$. The dashed lines show the
dissipative gap of the ideal two-qubit squeezing master equation, Eq. (2)
(with rates $\gamma_{1}=\gamma_{2}=\Gamma$). For small $\mu\ll\Gamma$ the
dissipative gap vanishes as $\kappa_{\rm slow}\sim\mu^{2}$, and for large
$\mu\gg\Gamma$ 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 $\hat{\mathcal{J}}_{1/2}^{\prime}$ (c.f. Eqs. (19) and
(20)).

##  Appendix C Balanced Two-Qubit Dissipators

For a setup described by Eq. (2), only in the infinite squeezing limit,
$r\rightarrow\infty$, is the steady-state a maximally entangled state.
However, for our system we have the freedom to choose the modulation
amplitudes, and in particular could choose
$\alpha_{-}=\beta_{-}=\alpha_{+}=\beta_{+}=\sqrt{\bar{m}}$. In this case the
master equation becomes

$\displaystyle\dot{\hat{\rho}}_{q}=\frac{4\bar{m}\bar{g}^{2}}{\kappa}\left(\mathcal{D}\left[\hat{\sigma}^{(1)}_{-}+\hat{\sigma}^{(2)}_{+}\right]+\mathcal{D}\left[\hat{\sigma}^{(1)}_{+}+\hat{\sigma}^{(2)}_{-}\right]\right)\hat{\rho}_{q}.$
(103)

As this is the $r\rightarrow\infty$ continuation of Eq. (2), one might be
tempted to believe that its steady-state is the maximally entangled Bell-state
$\ket{\Phi^{-}}=(\ket{00}-\ket{11})/\sqrt{2}$, regardless of the coherent
amplitude $\sqrt{\bar{m}}$, which now only controls the speed of the
dissipative dynamics.

However, while the state $\ket{\Phi^{-}}$ is indeed a steady-state of the
system, it is not the unique steady-state, and in fact any two-qubit state of
the form

$\displaystyle\hat{\rho}_{\rm
SS}=\nu\outerproduct{\Phi^{-}}{\Phi^{-}}+\frac{1}{4}(1-\nu)\hat{\mathbb{I}},$
(104)

is a steady-state of Eq. (103), where $-1/3\leq\nu\leq 1$ to ensure
$\hat{\rho}_{\rm SS}$ is a valid density matrix. This one parameter family of
states includes the maximally entangled state $\ket{\Phi^{-}}$ ($\nu=1$), but
also the maximally mixed state ($\nu=0$). As such, for the purposes of
generating steady-state qubit entanglement, it is important to apply
unbalanced drives to preserve the uniqueness of the steady-state, and avoid
the degeneracy described above.