# Dissipative stabilization of entangled cat states using a driven Bose-
Hubbard dimer

M. Mamaev  Department of Physics, McGill University, Montréal, Québec, Canada.
L. C. G. Govia  Institute for Molecular Engineering, University of Chicago,
5640 S. Ellis Ave., Chicago, IL 60637  A. A. Clerk  Institute for Molecular
Engineering, University of Chicago, 5640 S. Ellis Ave., Chicago, IL 60637

###### Abstract

We analyze a modified Bose-Hubbard model, where two cavities having on-site
Kerr interactions are subject to two-photon driving and correlated
dissipation. We derive an exact solution for the steady state of this
interacting driven-dissipative system, and use it show that the system permits
the preparation and stabilization of pure entangled non-Gaussian states, so-
called entangled cat states. Unlike previous proposals for dissipative
stabilization of such states, our approach requires only a linear coupling to
a single engineered reservoir (as opposed to nonlinear couplings to two or
more reservoirs). Our scheme is within the reach of state-of-the-art
experiments in circuit QED.

##  1 Introduction

To harness the full potential of continuous variable quantum computing
[Braunstein and van Loock(2005)], it is necessary to prepare and ideally
stabilize non-trivial entangled states of multiple bosonic modes. An ideal
approach here is reservoir engineering [Poyatos  _et al._(1996)Poyatos, Cirac,
and Zoller], where controlled interactions with tailored dissipation are used
to accomplish both these goals. The simplest kind of entangled bosonic states
are the two-mode squeezed states, and there exist many closely-related
reservoir engineering protocols which allow their stabilization (see e.g.
[Parkins  _et al._(2006)Parkins, Solano, and Cirac, Krauter  _et
al._(2011)Krauter, Muschik, Jensen, Wasilewski, Petersen, Cirac, and Polzik,
Muschik  _et al._(2011)Muschik, Polzik, and Cirac, Woolley and Clerk(2014)]).

While useful in certain applications, two-mode squeezed states are Gaussian,
and hence describable by a positive-definite Wigner function. The ability to
generate and stabilize truly non-Gaussian entangled states (which exhibit
negativity in their Wigner functions) would provide an even more powerful
potential resource. A canonical example of such states are two-mode entangled
Schrödinger cat states [Haroche and Raimond(2006)]. These states are integral
to bosonic codes for quantum computing [Vlastakis  _et al._(2013)Vlastakis,
Kirchmair, Leghtas, Nigg, Frunzio, Girvin, Mirrahimi, Devoret, and Schoelkopf,
Mirrahimi  _et al._(2014)Mirrahimi, Leghtas, Albert, Touzard, Schoelkopf,
Jiang, and Devoret, Wang  _et al._(2016)Wang, Gao, Reinhold, Heeres, Ofek,
Chou, Axline, Reagor, Blumoff, Sliwa, Frunzio, Girvin, Jiang, Mirrahimi,
Devoret, and Schoelkopf], and are equivalent to so-called “entangled coherent
states”, which have applications in many aspects of quantum information
[Sanders(2012)].

Dissipative stabilization of non-Gaussian states is in principle possible, but
requires more complicated setups than in the Gaussian case; in particular,
nonlinearity is now an essential resource. This complexity is present even
when trying to reservoir engineer cat states in a single mode: for example, a
recently implemented protocol requires one to engineer a two-photon loss
operator [Mirrahimi  _et al._(2014)Mirrahimi, Leghtas, Albert, Touzard,
Schoelkopf, Jiang, and Devoret, Leghtas  _et al._(2015)Leghtas, Touzard, Pop,
Kou, Vlastakis, Petrenko, Sliwa, Narla, Shankar, Hatridge, Reagor, Frunzio,
Schoelkopf, Mirrahimi, and Devoret]. There exist protocols for stabilizing
entangled cat states [Sarlette  _et al._(2012)Sarlette, Leghtas, Brune,
Raimond, and Rouchon, Arenz  _et al._(2013)Arenz, Cormick, Vitali, and Morigi]
based on passing a stream of few-level atoms through a pair of cavities. These
approaches effectively make use of either a complex nonlinear dissipator or
highly nonlinear cavity driving.

Figure 1: Left panel: Schematic of the system: two bosonic modes $\hat{a}$,
$\hat{b}$ are subjected to two-photon driving (amplitude $\lambda$), and have
on-site Kerr-interactions which are repulsive (attractive) for mode $\hat{a}$
(mode $\hat{b}$). The two modes couple to a single engineered Markovian
reservoir (e.g. a waveguide) which mediates the dissipative equivalent of a
tunneling interaction. For appropriately matched parameters (see main text),
the system has a unique pure steady state corresponding to an entangled
Schrödinger cat state. This state is an equal superposition of each cavity
being in an even-parity cat state, and each cavity being in an odd-parity cat
state. This is depicted schematically in the right panel, using Wigner
functions for the cat states involved in the superpositions.

In this paper, we present a surprisingly simple scheme for stabilizing a two-
mode entangled cat-state. Our setup is sketched in Fig. 1, and resembles the
well-studied dissipative Bose-Hubbard dimer model (see e.g. [Pudlik  _et
al._(2013)Pudlik, Hennig, Witthaut, and Campbell, Grujic  _et
al._(2013)Grujic, Clark, Jaksch, and Angelakis, Cao  _et al._(2016)Cao,
Mahmud, and Hafezi, Casteels and Ciuti(2017), Casteels and Wouters(2017), Noh
and Angelakis(2017)]), where two driven bosonic modes with on-site Kerr
interactions are coupled via coherent tunneling. We make some crucial
modifications to this standard setup. In our system, both cavities are subject
to coherent two-photon (i.e. parametric) driving, as opposed to the more
standard linear driving. Further, there is no direct tunneling between the
cavities. Instead, both cavities couple linearly to a common dissipative
environment (the engineered reservoir), which mediates the dissipative
equivalent of a tunnel interaction. This dissipative coupling could be
generated in many different ways, including by simply coupling both cavities
passively to a low-loss waveguide [Chang  _et al._(2012)Chang, Jiang,
Gorshkov, and Kimble, Metelmann and Clerk(2015)]. Finally, we require one
cavity to have an attractive interaction, while the other has a repulsive
interaction.

Despite being an interacting open system, we are able to analytically solve
for the stationary state, something that is not possible for the standard
driven-dissipative Bose Hubbard dimer [Cao  _et al._(2016)Cao, Mahmud, and
Hafezi]. Using our exact solution, we identify a simple parameter regime where
the steady state is a pure entangled cat state. Our scheme is also easily
modified so that it stabilizes a two-dimensional manifold spanned by two even-
parity cat states; this allows a potential two-cavity generalization of the
cat-state encoded qubit introduced in Ref. [Mirrahimi  _et
al._(2014)Mirrahimi, Leghtas, Albert, Touzard, Schoelkopf, Jiang, and
Devoret]. We stress that all the key ingredients required for our scheme are
available in state of the art circuit QED setups (namely parametric drives,
see e.g. [Yamamoto  _et al._(2008)Yamamoto, Inomata, Watanabe, Matsuba,
Miyazaki, Oliver, Nakamura, and Tsai, Krantz  _et al._(2016)Krantz, Bengtsson,
Simoen, Gustavsson, Shumeiko, Oliver, Wilson, Delsing, and Bylander, Mutus
_et al._(2014)Mutus, White, Barends, Chen, Chen, Chiaro, Dunsworth, Jeffrey,
Kelly, Megrant, Neill, O’Malley, Roushan, Sank, Vainsencher, Wenner,
Sundqvist, Cleland, and Martinis], strong on-site Kerr interactions [Kirchmair
_et al._(2013)Kirchmair, Vlastakis, Leghtas, Nigg, Paik, Ginossar, Mirrahimi,
Frunzio, Girvin, and Schoelkopf] and engineered dissipation [Sliwa  _et
al._(2015)Sliwa, Hatridge, Narla, Shankar, Frunzio, Schoelkopf, and Devoret,
Lecocq  _et al._(2017)Lecocq, Ranzani, Peterson, Cicak, Simmonds, Teufel, and
Aumentado, Lecocq  _et al._(2016)Lecocq, Clark, Simmonds, Aumentado, and
Teufel]). While the relaxation rate associated with our stabilization scheme
can become very slow for large photon number entangled cat states, we show
explicitly that our scheme is still effective in state-of-the-art systems and
resilient to low levels of internal cavity loss.

##  2 Single-mode cat state stabilization

###  2.1 Dissipative stabilization of a single-mode squeezed state

We introduce our scheme by first reviewing the well-understood reservoir
engineering protocol for stabilizing a single-mode squeezed state [Cirac  _et
al._(1993)Cirac, Parkins, Blatt, and Zoller, Rabl  _et al._(2004)Rabl,
Shnirman, and Zoller, Kronwald  _et al._(2013)Kronwald, Marquardt, and Clerk].
This protocol was recently realized experimentally both in optomechanics
[Wollman  _et al._(2015)Wollman, Lei, Weinstein, Suh, Kronwald, Marquardt,
Clerk, and Schwab, Lecocq  _et al._(2015)Lecocq, Clark, Simmonds, Aumentado,
and Teufel, Pirkkalainen  _et al._(2015)Pirkkalainen, Damskägg, Brandt,
Massel, and Sillanpää] as well as in a trapped ion system [Kienzler  _et
al._(2015)Kienzler, Lo, Keitch, de Clercq, Leupold, Lindenfelser, Marinelli,
Negnevitsky, and Home]. The protocol is based on the observation that a
squeezed state is the vacuum of a bosonic Bogoliubov operator
$\hat{\beta}[r]=\cosh(r)\hat{d}+\sinh(r)\hat{d}^{\dagger}$, where $\hat{d}$ is
the annihilation operator for the mode of interest. Hence, by simply cooling
the mode $\hat{\beta}[r]$, one can stabilize a squeezed state.

The required dissipative cooling dynamics is obtained most simply from a
Lindblad master equation for the reduced density matrix $\hat{\rho}$ of the
mode:

$\frac{d}{dt}\hat{\rho}=\Gamma\left(\hat{z}\hat{\rho}\hat{z}^{\dagger}-\frac{1}{2}\\{\hat{z}^{\dagger}\hat{z},\hat{\rho}\\}\right)\equiv\Gamma\mathcal{L}[\hat{z}]\hat{\rho},$
(1)

with a jump operator $\hat{z}$ given by

$\hat{z}\equiv\mu\hat{d}+\nu\hat{d}^{\dagger}.$ (2)

$\Gamma$ represents the rate of the dissipative interaction (i.e. the cooling
rate), and the parameters $\mu,\nu$ are set to

$\mu=\cosh r,\,\,\,\nu=\sinh r,$ (3)

making the jump operator $\hat{z}$ equal to the Bogoliubov-mode annihilation
operator $\hat{\beta}[r]$. The unique steady state of this master equation is
then the dark state of the jump operator $\hat{z}$, which is nothing else than
the vacuum squeezed state with squeeze parameter $r$.

###  2.2 Dissipative stabilization of a single-mode cat state

To modify this scheme to stabilize a non-Gaussian state, we now introduce
nonlinearity into the jump operator $\hat{z}$. One simple way to do this is to
make the parameter $\mu$ a photon-number dependent operator. We take
$\mu=\mu_{1}\hat{d}^{\dagger}\hat{d}$, and consider a master equation of the
form in Eq. (1) with jump operator

$\hat{z}=\left(\mu_{1}\hat{d}^{\dagger}\hat{d}\right)\hat{d}+\nu\hat{d}^{\dagger}.$
(4)

The steady states of this master equation correspond to the dark states of the
jump operator $\hat{z}$. Despite the nonlinearity, these dark states are
easily found: they correspond to a two-dimensional subspace spanned by even
and odd Schrödinger cat states. Defining even/odd cat-states in terms of
coherent states $\ket{\alpha}$ as

$\Ket{\mathcal{C}_{\pm}(\alpha)}=\frac{\ket{\alpha}\pm\ket{-\alpha}}{\sqrt{2(1\pm
e^{-2|\alpha|^{2}})}},$ (5)

we find that our jump operator has dark states

$\displaystyle\hat{z}\ket{\mathcal{C}_{\pm}(\sqrt{2}\overline{\alpha})}=0,$
(6)

where the amplitude $\overline{\alpha}$ is given by

$\overline{\alpha}=i\sqrt{\nu/(2\mu_{1})}.$ (7)

For some applications, the ability to stabilize a manifold of dark-states can
be extremely useful, see e.g. [Mirrahimi  _et al._(2014)Mirrahimi, Leghtas,
Albert, Touzard, Schoelkopf, Jiang, and Devoret]. However, for the simplest
state preparation applications, it is ideal to have a unique steady state. We
would thus like to break the symmetry between the even and odd cat-states
above, and stabilize, e.g., only the even parity cat. To do this, we will
simply interpolate between the jump operator in Eq. (2) (which stabilizes an
even-parity squeezed state) and the nonlinear jump operator in Eq. (4). We
thus consider a modified jump operator

$\hat{z}=\left(\mu_{0}+\mu_{1}\hat{d}^{\dagger}\hat{d}\right)\hat{d}+\nu\hat{d}^{\dagger}\equiv\hat{z}_{\rm
cat}.$ (8)

As we show explicitly in Appendix A, for $\mu_{0}\neq 0$, the only possible
dark states of this dissipator are even-parity states. We thus find a unique
dark state (and hence pure steady state of Eq. (1)):

$\hat{z}_{\rm cat}\Ket{\tilde{C}[\mu_{0},\mu_{1},\nu]}=0.$ (9)

One can obtain an analytic expression for this state, we discuss this further
in Sec. 4 (c.f. Eq. (22)). As could already be expected from the results
above, we find that in the limit of small $\mu_{0}$, this state asymptotically
approaches an even-parity cat state:

$\lim_{\mu_{0}/\mu_{1}\rightarrow
0}\left|\Braket{\tilde{C}[\mu_{0},\mu_{1},\nu]}{\mathcal{C}_{+}[\sqrt{2}\overline{\alpha}]}\right|\rightarrow
1.$ (10)

Even for small but finite $\mu_{0}/\mu_{1}$,
$\ket{\tilde{C}[\mu_{0},\mu_{1},\nu]}$ has an extremely high fidelity with the
ideal even parity cat-state (c.f. Fig. 2). We have thus shown how (in
principle) the usual reservoir engineering recipe for stabilizing a single-
mode squeezed state can be extended to stabilizing an even-parity cat state.

###  2.3 Auxiliary cavity method for realizing a nonlinear dissipator

We now consider how one might realize the nonlinear jump operator in Eq. (8).
The usual recipe is to introduce a second highly damped mode $\hat{c}$ which
couples nonlinearly to the principle mode $\hat{d}$ of interest. To that end,
we consider a two-mode system whose non-dissipative dynamics is described by
the Hamiltonian

$\displaystyle\hat{H}_{1}$
$\displaystyle=\Lambda\left(\hat{c}^{\dagger}\hat{z}_{\rm cat}+\hat{z}_{\rm
cat}^{\dagger}\hat{c}\right)$
$\displaystyle=\Lambda\hat{c}^{\dagger}\left[\mu_{0}\hat{d}+\mu_{1}(\hat{d}^{\dagger}\hat{d})\hat{d}+\nu\hat{d}^{\dagger}\right]+h.c.$
(11)

The quadratic terms in this Hamiltonian correspond to a simple beam-splitter
interaction ($\propto\mu_{0}$) and two-mode squeezing interaction
($\propto\nu$). The remaining nonlinear term ($\propto\mu_{1}$) describes
photon-number dependent tunneling.

We now add damping of the $\hat{c}$ mode (loss at rate $\gamma$), such that
the full evolution of the two modes is given by

$\frac{d}{dt}\hat{\rho}=-i[\hat{H_{1}},\hat{\rho}]+\gamma\mathcal{L}[\hat{c}]\hat{\rho}.$
(12)

In the limit where the $\hat{c}$ mode is heavily damped (i.e.
$\gamma\gg\Lambda$), standard adiabatic elimination yields a master equation
of the form in Eq. (1) with jump operator $\hat{z}_{\rm cat}$ and with
$\Gamma\propto\Lambda^{2}/\gamma$. In this limit, our two-mode system has a
unique steady state analogous to that in Eq. (9):

$\displaystyle\ket{\psi_{\rm ss}}$
$\displaystyle=\ket{0}_{c}\ket{\tilde{C}[\mu_{0},\mu_{1},\nu]}_{d}.$ (13)

As discussed earlier, the state $\ket{\tilde{C}[\mu_{0},\mu_{1},\nu]}_{d}$
asymptotically approaches an even cat-state in the limit
$\mu_{0}/\mu_{1}\rightarrow 0$.

We next make a crucial observation: even if one is far from the large-damping
limit, the above state is still the steady state of the full master equation
in Eq. (12). This follows directly from the fact that $\hat{H}_{1}$ and
$\hat{c}$ both annihilate this state. Thus, despite the nonlinearity in our
two-mode system, we have found its steady state (analytically, c.f. Eq. (22))
for all values of the parameters $\gamma,\Lambda,\mu_{0},\mu_{1}$ and $\nu$.

##  3 Stabilization of two-mode entangled cat states

The previous section described a setup for stabilizing a single-mode cat state
in a cavity, using a nonlinear interaction with a damped auxiliary cavity. We
will now connect this result to the stated goal of this paper: stabilization
of non-Gaussian entangled states using a system which resembles a driven-
dissipative Bose-Hubbard dimer model. Surprisingly, we will find that the
resulting setup is much simpler than that described in Eq. (11).

We again start with our two modes $\hat{c}$ and $\hat{d}$, with the
dissipative dynamics described in the previous section. We now imagine that
these modes are delocalized modes, describing equal superpositions of two
distinct localized cavity modes $\hat{a}$ and $\hat{b}$:

$\hat{a}=\frac{1}{\sqrt{2}}(\hat{c}+\hat{d}),\>\>\hat{b}=\frac{1}{\sqrt{2}}(\hat{c}-\hat{d}).$
(14)

The dissipative dynamics described in Eq. (12) stabilizes the state
$\ket{\psi_{\rm ss}}$ defined in Eq. (13). Transforming this state to the
“localized” $\hat{a}$, $\hat{b}$ mode basis is equivalent to passing the state
through a beam splitter. It is well known that two-mode entangled coherent
states can be generated by sending single-mode cat states into the arms of a
$50/50$ beam-splitter [Asbóth  _et al._(2004)Asbóth, Adam, Koniorczyk, and
Janszky]. We thus find that in the localized basis (and for
$\mu_{0}\ll\mu_{1}$), our steady state has the form of a two-mode entangled
cat state:

$\ket{\psi_{\mathrm{ss}}}\approx\frac{1}{\sqrt{2}}\left[\ket{\mathcal{C}_{+}(\overline{\alpha})}_{a}\ket{\mathcal{C_{+}}(\overline{\alpha})}_{b}-\ket{\mathcal{C}_{-}(\overline{\alpha})}_{a}\ket{\mathcal{C_{-}}(\overline{\alpha})}_{b}\right].$
(15)

To determine the feasibility of our scheme, we need to express the Hamiltonian
$\hat{H}_{1}$ in terms of the localized mode operators $\hat{a},\ \hat{b}$.
Before doing this, we first modify $\hat{H}_{1}$ to make it more symmetric,
but without changing the steady state:

$\displaystyle\hat{H}_{2}$
$\displaystyle=\hat{H}_{1}+\left(\Lambda\mu_{1}\hat{c}^{\dagger}\left(\hat{c}^{\dagger}\hat{c}\right)\hat{d}+h.c.\right)$
$\displaystyle=\Lambda\hat{c}^{\dagger}\left[\mu_{0}\hat{d}+\mu_{1}(\hat{c}^{\dagger}\hat{c}+\hat{d}^{\dagger}\hat{d})\hat{d}+\nu\hat{d}^{\dagger}\right]+h.c.$
(16)

We again consider the master equation

$\frac{d}{dt}\hat{\rho}=-i[\hat{H_{2}},\hat{\rho}]+\gamma\mathcal{L}[\hat{c}]\hat{\rho}.$
(17)

It is easy to see that the state $\ket{\psi_{\rm ss}}$ (c.f. Eq. (13))
continues to be the steady state of this master equation; this follows from
the fact that the term we added to the Hamiltonian also annihilates this
state.

We now consider the form of the nonlinear Hamiltonian $\hat{H}_{2}$ in the
localized mode basis. Defining for convenience

$\displaystyle\Delta=\Lambda\mu_{0},\>\>\lambda=\Lambda\nu/2,\>\>U=\Lambda\mu_{1},$
(18)

we obtain

$\displaystyle\frac{d}{dt}\hat{\rho}=-i[\hat{H}_{2},\hat{\rho}]+\gamma\mathcal{L}[\hat{a}+\hat{b}]\hat{\rho},$
(19)

with

$\displaystyle\hat{H}_{2}$
$\displaystyle=\Delta\hat{a}^{\dagger}\hat{a}+\lambda(\hat{a}^{\dagger}\hat{a}^{\dagger}+h.c.)+U\hat{a}^{\dagger}\hat{a}^{\dagger}\hat{a}\hat{a}$
$\displaystyle-\Delta\hat{b}^{\dagger}\hat{b}-\lambda(\hat{b}^{\dagger}\hat{b}^{\dagger}+h.c.)-U\hat{b}^{\dagger}\hat{b}^{\dagger}\hat{b}\hat{b}.$
(20)

We see that the required Hamiltonian has a remarkably simple form,
corresponding to the schematic in Fig. 1. It describes two uncoupled nonlinear
driven cavities whose Hamiltonians differ only by a sign. Each cavity has a
Kerr nonlinearity (strength $U,-U$) and is subject to a parametric drive
(drive detunings $\pm\Delta$, drive amplitudes $\pm\lambda$). The only
interaction between the two cavities comes from the non-local (but linear)
dissipator in the master equation, arising from a linear interaction with a
common reservoir. Such dissipators have been discussed extensively and even
realized experimentally (e.g. [Fang  _et al._(2017)Fang, Luo, Metelmann,
Matheny, Marquardt, Clerk, and Painter]) in many different contexts. For the
simple case where each cavity has the same resonance frequency, it can be
generated by coupling both cavities to a transmission line with an
appropriately chosen distance between the cavities [Chang  _et
al._(2012)Chang, Jiang, Gorshkov, and Kimble, Metelmann and Clerk(2015)].
Further details about implementation issues (including how to realize a
repulsive Kerr interaction in circuit QED) are discussed in Appendix C.

Eqs. (19)-(20) thus completely describe our dissipative scheme for the
preparation and stabilization of entangled non-Gaussian states. Note that the
emergence of cat states and their adiabatic preparation in a single
parametrically-driven Kerr cavity is the subject of several recent works
[Goto(2016a), Goto(2016b), Minganti  _et al._(2016)Minganti, Bartolo, Lolli,
Casteels, and Ciuti, Puri  _et al._(2017)Puri, Boutin, and Blais]; our scheme
leverages similar resources in two cavities (along with engineered
dissipation) to now stabilize an entangled cat state. Note that the amplitude
$\overline{\alpha}$ of our cat state (c.f. Eq. (7)) in terms of the energy
scales in the Hamiltonian is now $\overline{\alpha}=i\sqrt{\lambda/U}$.

##  4 Exact solution for the steady state

We now explore in more detail our analytic expression for the pure steady
state of our system. For simplicity, we work in the non-local $\hat{c},\
\hat{d}$ mode basis. In this basis the steady state $\ket{\psi_{\rm ss}}$ has
the $\hat{c}$ mode in vacuum, while the $\hat{d}$ mode has a definite even
photon number parity (c.f. Eq. (9)). By solving a simple recursion relation
(see Appendix A), the even-parity state of the $\hat{d}$ mode is given by:

$\Ket{\tilde{C}[\mu_{0},\mu_{1},\nu]}_{d}=\mathcal{N}\sum_{n=0}^{\infty}\tilde{c}_{2n}\ket{2n}_{d},$
(21)
$\tilde{c}_{2n}=\left(\frac{-\nu}{2\mu_{1}}\right)^{n}\sqrt{\frac{(2n-1)!!}{(2n)!!}}\frac{\Gamma\left(\frac{1}{2}+\frac{\mu_{0}}{2\mu_{1}}\right)}{\Gamma\left(n+\frac{1}{2}+\frac{\mu_{0}}{2\mu_{1}}\right)}.$
(22)

Here $\ket{2n}_{d}$ denotes a Fock state with $2n$ photons, $\Gamma(z)$ is the
Gamma function, and $\mathcal{N}$ is an overall normalization constant. Our
ability to easily find an analytic form relies on the symmetry of the problem
and the existence of a pure steady state; it is much simpler than more general
methods based on the positive-P function [Drummond and Walls(1980), Minganti
_et al._(2016)Minganti, Bartolo, Lolli, Casteels, and Ciuti, Cao  _et
al._(2016)Cao, Mahmud, and Hafezi]. Our solution method is reminiscent of that
used in Ref. [Stannigel  _et al._(2012)Stannigel, Rabl, and Zoller], which
studied dissipative entanglement in cascaded quantum systems; in contrast to
that work, our approach does not require cascaded interactions, or even any
breaking of time-reversal symmetry.

Figure 2: (a-c) Wigner functions of the steady state (c.f. Eq. (13)) for the
non-local mode $\hat{d}$ at snapshots of $\Delta/U$. This is equivalent to a
cut through the two-mode Wigner function for the local modes, defined by the
combined quadratures $(\hat{x}_{a}-\hat{x}_{b})/\sqrt{2}$ and
$(\hat{p}_{a}-\hat{p}_{b})/\sqrt{2}$, where $\hat{x}_{k}$ and $\hat{p}_{k}$
are the usual quadrature operators for mode $k$. As $\Delta/U$ increases, the
evolution from an even-parity cat state to a squeezed state is clearly
visible. (d) Fidelity of the exact steady state (c.f. Eq. (21)) with both the
squeezed state (dashed lines) obtained for $\Delta/U\rightarrow\infty$, and
the even cat-state (solid lines) obtained for $\Delta/U\rightarrow 0$ (c.f.
Eq. (15)). The fidelity is defined by
$F\left(\ket{\psi_{1}},\ket{\psi_{2}}\right)=\left|\bra{\psi_{1}}\psi_{2}\rangle\right|^{2}$.

The analytic expression for $\ket{\psi_{\rm ss}}$ allows us to confirm the
expected parameter dependence of the state. For $\mu_{0}\gg\mu_{1}$,
nonlinearity plays a minor role, and the state tends to a Gaussian squeezed
state. In contrast, for $\mu_{0}\ll\mu_{1}$, the nonlinearity is crucial to
the existence of a steady state, and $\ket{\psi_{\rm ss}}$ tends to the
desired even cat-state. Fig. 2 uses the exact expression of the steady state
$\ket{\tilde{C}}_{d}$ to show this evolution as a function of the ratio
$\Delta/U=\mu_{0}/\mu_{1}$. The upper panel depicts the Wigner function for
the $\hat{d}$ mode at three specific values of $\Delta/U$. As described in the
figure caption, this describes a cut through the two-mode Wigner function
describing the state in the localized $\hat{a},\ \hat{b}$ basis. One clearly
sees the cross-over from an even cat-state to a squeezed state as $\Delta$ is
increased. This figure also demonstrates that even for modest values of
$\Delta/U$, the fidelity with the desired even-parity cat state can be
extremely high.

##  5 Effect of imperfections

While our ideal system always possesses an entangled pure steady state,
imperfections in any realistic implementation will cause deviations from the
desired dynamics. The most important imperfection will come in the form of
unwanted internal losses in each cavity. To assess the impact of such losses,
we add single-photon loss (at a rate $\kappa$ in each cavity) to our master
equation:

$\frac{d}{dt}\hat{\rho}=-i[\hat{H}_{2},\hat{\rho}]+\gamma\mathcal{L}[\hat{a}+\hat{b}]\rho+\kappa\mathcal{L}[\hat{a}]\hat{\rho}+\kappa\mathcal{L}[\hat{b}]\hat{\rho}.$
(23)

With the inclusion of single particle loss, the steady state will no longer be
pure, and we cannot find an analytic closed form solution. We instead
numerically simulate the master equation. At a heuristic level, the
sensitivity to weak internal loss will depend on the typical timescale
required by the ideal system to reach its steady state. As discussed in detail
in Appendix B, this timescale can be extremely long, scaling as
$e^{4|\overline{\alpha}|^{2}}$. As a result, low internal loss levels
(relative to $U$, $\lambda$, and $\Delta$) are required for good performance.

Fig. 3 plots the fidelity of the numerically-obtained steady state with both
the ideal, loss-free steady state $\ket{\psi_{\rm ss}}$ (c.f. Eq. (13)), as
well as with the ideal two-mode entangled cat-state of Eq. (15), as a function
of the internal loss rate $\kappa$. We find that for $\Delta/U=1$, we get
$F\approx 0.95$ with the exact solution and $F\approx 0.86$ with the ideal cat
state for a loss rate $\kappa=10^{-3}U$. This value of loss is compatible with
state-of-the-art experiments in circuit QED. An attractive Kerr nonlinearity
in a 3D superconducting microwave cavity at the level of $U\sim$ 100-1000 kHz
has been realized [Kirchmair  _et al._(2013)Kirchmair, Vlastakis, Leghtas,
Nigg, Paik, Ginossar, Mirrahimi, Frunzio, Girvin, and Schoelkopf], while in
state-of-the-art 3D cavities internal leakage rates are on the order of 10-100
Hz [Reagor  _et al._(2013)Reagor, Paik, Catelani, Sun, Axline, Holland, Pop,
Masluk, Brecht, Frunzio, Devoret, Glazman, and Schoelkopf, Reagor  _et
al._(2016)Reagor, Pfaff, Axline, Heeres, Ofek, Sliwa, Holland, Wang, Blumoff,
Chou, Hatridge, Frunzio, Devoret, Jiang, and Schoelkopf].

Note that for the results in Fig. (3), we have taken the coupling $\gamma$ to
the engineered dissipation to be $\gamma=\sqrt{2}\Delta$. As discussed in more
detail in Appendix B, this choice minimizes the timescale required to prepare
the steady state, and thus helps minimize the effects of internal loss.
However, the fidelity is not extremely sensitive to this precise tuning of the
value of $\gamma$.

Figure 3: Effect of intrinsic cavity losses $\kappa$ (identical for both
cavities, c.f. Eq. (23)) on the steady state, for various values of the
parametric drive detuning parameter $\Delta$. Solid lines are the fidelity of
the numerically-obtained steady state including loss, with the ideal loss-free
steady state of Eq. (21). Dashed lines are the fidelity with the two-mode
entangled cat state of Eq. (15) (which coincides with the ideal steady state
when $\Delta/U\rightarrow 0$). For all curves the parametric drive strength is
set to $\lambda/U=2$ (thus setting the amplitude
$\overline{\alpha}=i\sqrt{2}$), and the coupling rate $\gamma$ to the
correlated dissipation is taken to be $\gamma=\sqrt{2}\Delta$ to maximize the
speed of internal system dynamics (see Appendix B). The fidelity is defined by
$F(\rho_{1},\rho_{2})=\left|\text{Tr}\left[\sqrt{\sqrt{\rho_{1}}\rho_{2}\sqrt{\rho_{2}}}\right]\right|^{2}$.

A qualitative description of the effect of internal loss can be found in
Appendix D, where we show the Wigner functions of an example steady state with
and without internal loss. Internal loss causes a reduction in the magnitude
and area of Wigner function negativity, indicating a reduction in quantum
correlations between the two modes. For realistic internal loss rates this
reduction is minimal, and the state remains entangled.

Finally, while the ideal system requires a matching of the Hamiltonian
parameters of the two cavities, the desired behaviour is still found for small
mismatches. For the same parameter choices as in Fig. 3, a drive or
nonlinearity strength mismatch by 5% or less still produces a steady state
with $F>0.95$. Further details can be found in Appendix E.

##  6 Conclusion

We have a presented a relatively simple scheme that uses a modified version of
a driven-dissipative Bose Hubbard dimer to stabilize a two-cavity entangled
cat state. The scheme uses detuned parametric drives applied to each cavity,
and couples them only through a linear dissipative tunneling interaction that
is mediated by an external reservoir. While we have focused on the case where
there is a unique steady state, by setting the parametric drive detunings
$\Delta=0$, one would instead stabilize a two dimensional manifold spanned by
the entangled cat states (the states $\ket{2}$ and $\ket{1}$ defined in Eq.
(27)).

##  7 Acknowledgements

We acknowledge R. P. Tiwari and J. M. Torres for fruitful discussions. This
work was supported by NSERC, and by the AFOSR MURI FA9550-15-1-0029.

## APPENDICES

##  Appendix A Analytic expression for the steady state

We derive here the exact steady state of the system described by the
Hamiltonian $\hat{H}_{2}$ in Eq. (16) and the master equation in Eq. (17).
Motivated by the fact that only the $\hat{c}$ mode experiences any damping, we
make an ansatz that there exists a pure-state steady state in which the
$\hat{c}$ mode is in vacuum, i.e.

$\ket{\psi}=\ket{0}_{c}\otimes\sum_{n=0}^{\infty}\alpha_{n}\ket{n}_{d}.$ (24)

Here, $\alpha_{n}$ are arbitrary normalized coefficients and $\ket{n}_{d}$ are
Fock states.

The above state is trivially a dark-state of the jump operator $\hat{c}$
appearing in the dissipator of the master equation in Eq. (17). To be a steady
state, it must then necessarily also be an eigenstate of the Hamiltonian.
Acting on $\ket{\psi}$ with $\hat{H}_{2}$ yields,

$\displaystyle\hat{H}_{2}\ket{\psi}=\ket{1}_{c}\otimes\sum_{n=0}^{\infty}$
$\displaystyle\alpha_{n}\big{(}\left[\Delta+U(n-1)\right]\sqrt{n}\ket{n-1}_{d}$
$\displaystyle+2\lambda\sqrt{n+1}\ket{n+1}_{d}\big{)}.$ (25)

where we have expressed $\mu_{0},\mu_{1}$ and $\nu$ in terms of $\Delta,U$ and
$\lambda$ using Eq. (18).

The RHS of Eq. (25) has the $\hat{c}$ mode in the state $\ket{1}$, while in
$\ket{\psi}$ it is in the state $\ket{0}$. Thus, the only way that
$\ket{\psi}$ can be an eigenstate of $\hat{H}_{2}$ is if it has zero energy.
We thus need $\hat{H}_{2}$ to annihilate $\ket{\psi}$, and therefore require
that the coefficient of each $\ket{n}_{d}$ state in Eq. (25) vanish. This
leads to the recursion relation

$\left(\Delta+Un\right)\sqrt{n+1}\alpha_{n+1}=\begin{cases}0&n=0,\\\
-2\lambda\sqrt{n}\alpha_{n-1}&n\geq 1.\end{cases}$ (26)

The above recursion relation does not mix the even and odd-parity subspaces.
Note that for $\Delta\neq 0$, the $n=0$ case of Eq. (22) forces
$\alpha_{1}=0$, and as a result, forces all odd $\alpha_{j}$ to be zero. The
solution thus has only even-parity Fock states. Solving Eq. (26) for the
remaining even parity coefficients directly gives the solution presented in
Eq. (22), with $\alpha_{2n}\equiv\tilde{c}_{n}$.

Note also that if $\Delta=0$, then the $n=0$ case of Eq. (22) is trivially
satisfied regardless of the value of $\alpha_{1}$, and one can find both even
and odd parity steady states. Thus, as discussed in the main text, a non-zero
detuning $\Delta$ is crucial in order to obtain a unique steady state.

While the above discussion has identified a pure steady state for our system
when $\Delta\neq 0$, the question remains whether this is a unique steady
state. In the case where $U=0$, the system is linear and one can exactly solve
the system. The linear system is stable for $\lambda\leq\Delta$, and in this
regime the unique steady state is exactly that given by Eq. (22). For non-zero
$U$, we have no rigorous analytic argument precluding the existence of
multiple steady states (i.e. in addition to the state we describe analytically
here). We have however investigated the model numerically over a wide range of
parameters, and find that as long as $\Delta$ is non-zero, the steady state
described by Eq. (21) is indeed the unique steady state.

##  Appendix B Dynamical timescales

We consider in this appendix the parameter dependence of the slow rates of the
master equation in Eq. (19) that determine the typical time needed to prepare
our non-Gaussian entangled steady state. As discussed, this timescale sets the
sensitivity to unwanted perturbations such as internal loss. Note first that
when the drive detuning $\Delta=0$, each cavity has exact zero-energy
eigenstates $\ket{\mathcal{C}_{\pm}(\overline{\alpha})}$ [Goto(2016a),
Goto(2016b), Puri  _et al._(2017)Puri, Boutin, and Blais]. For the interesting
limit of a large cat state amplitude $\overline{\alpha}$, the slowest
timescales in our system correspond to dynamics within the four dimensional
subspace spanned by these cat states. Entangled cat states form an orthonormal
basis for this subspace:

$\displaystyle\ket{\tilde{2}}$
$\displaystyle=\frac{1}{\sqrt{2}}\left(\ket{\mathcal{C}_{+}(\overline{\alpha})}\ket{\mathcal{C}_{+}(\overline{\alpha})}-\ket{\mathcal{C}_{-}(\overline{\alpha})}\ket{\mathcal{C}_{-}(\overline{\alpha})}\right),$
$\displaystyle\ket{\tilde{4}}$
$\displaystyle=\frac{1}{\sqrt{2}}\left(\ket{\mathcal{C}_{+}(\overline{\alpha})}\ket{\mathcal{C}_{+}(\overline{\alpha})}+\ket{\mathcal{C}_{-}(\overline{\alpha})}\ket{\mathcal{C}_{-}(\overline{\alpha})}\right),$
$\displaystyle\ket{1}$
$\displaystyle=\frac{1}{\sqrt{2}}\left(\ket{\mathcal{C}_{+}(\overline{\alpha})}\ket{\mathcal{C}_{-}(\overline{\alpha})}-\ket{\mathcal{C}_{-}(\overline{\alpha})}\ket{\mathcal{C}_{+}(\overline{\alpha})}\right),$
$\displaystyle\ket{3}$
$\displaystyle=\frac{1}{\sqrt{2}}\left(\ket{\mathcal{C}_{+}(\overline{\alpha})}\ket{\mathcal{C}_{-}(\overline{\alpha})}+\ket{\mathcal{C}_{-}(\overline{\alpha})}\ket{\mathcal{C}_{+}(\overline{\alpha})}\right).$
(27)

The states $\ket{\tilde{2}},\ket{\tilde{4}}$ have even total photon number
parity, while the states $\ket{1},\ket{3}$ have an odd total photon number
parity.

It is useful to make a slight rotation of the even parity states:

$\displaystyle\ket{2}$
$\displaystyle\equiv\cos\theta\ket{\tilde{2}}+\sin\theta\ket{\tilde{4}},$ (28)
$\displaystyle\ket{4}$
$\displaystyle\equiv-\sin\theta\ket{\tilde{2}}+\cos\theta\ket{\tilde{4}}.$
(29)

where for large $|\overline{\alpha}|$ we have

$\sin\theta\simeq 2e^{-2|\overline{\alpha}|^{2}}.$ (30)

Note that in the large $|\overline{\alpha}|$ limit,
$\ket{2}\rightarrow\ket{\tilde{2}}$, which is our expected entangled-cat
steady state.

The states $\ket{1},\ket{2},\ket{3},\ket{4}$ will form a useful basis for our
cat state subspace. One finds that there are two distinct dark states in this
subspace: they correspond to the basis states $\ket{1}$ and $\ket{2}$. If we
now restrict our master equation to this cat-state subspace, we find the
effective dynamics shown schematically in Fig. 4. The detuning term in
$\hat{H}_{2}$ becomes a coherent tunneling term of amplitude $\tilde{\Delta}$
between the odd parity states $\ket{1}$ and $\ket{3}$, whereas the dissipator
gives rise to the three depicted incoherent transitions (with rates
$\Gamma_{3,4},\Gamma_{4,3},\Gamma_{2,3}$). One sees clearly from this
schematic that, as expected, the dark state $\ket{2}$ will be our eventual
steady state.

For large $|\overline{\alpha}|$, one finds that $\tilde{\Delta}$ and
$\Gamma_{2,3}$ become exponentially small:

$\displaystyle\tilde{\Delta}$
$\displaystyle\simeq-4|\overline{\alpha}|^{2}e^{-2|\overline{\alpha}|^{2}}\Delta,$
(31) $\displaystyle\Gamma_{2,3}$ $\displaystyle\simeq
16|\overline{\alpha}|^{2}e^{-4|\overline{\alpha}|^{2}}\gamma.$ (32)

The exponential suppression of these rates corresponds to the extremely weak
overlap of the coherent states $\ket{\overline{\alpha}}$ and
$\ket{-\overline{\alpha}}$. The above processes will then form the bottleneck
in having populations relax to the state $\ket{2}$, and will determine the
slow timescales characterizing state preparation.

The superoperator describing this reduced master equation can be readily
diagonalized, yielding the relevant dynamical rates. In the large
$\overline{\alpha}$ limit, we find that the rate-limiting slow relaxation
between the dark states $\ket{1}$ and $\ket{2}$ is described by

$\Gamma_{\rm
rel}=\frac{1}{2}\left(2\Gamma_{\Delta}+\Gamma_{2,3}-\sqrt{\left(2\Gamma_{\Delta}\right)^{2}+\left(\Gamma_{2,3}\right)^{2}}\right),$
(33)

where we have defined the effective incoherent tunneling rate between the
state $\ket{1}$ and the $\\{\ket{3}$, $\ket{4}\\}$ bright-state manifold as:

$\Gamma_{\Delta}=\frac{4\tilde{\Delta}^{2}}{\Gamma_{2,3}+\Gamma_{4,3}}\simeq
16|\overline{\alpha}|^{2}e^{-4|\overline{\alpha}|^{2}}\frac{\Delta^{2}}{\gamma}.$
(34)

The rate $\Gamma_{\rm rel}$ describes an effective two step process where
population first incoherently tunnels from $\ket{1}$ to the bright state
manifold (rate $\Gamma_{\Delta}$), and then relaxes to the steady state
$\ket{2}$ (rate $\Gamma_{2,3}$). Note that $\Gamma_{\rm
rel}\neq\left(\Gamma_{\Delta}^{-1}+\Gamma_{2,3}^{-1}\right)^{-1}$, as the
incoherent tunneling process described by $\Gamma_{\Delta}$ can also move
population the wrong way, i.e. from $\ket{3}$ back to the state $\ket{1}$.

It is instructive to consider the behaviour of $\Gamma_{\rm rel}$ for small
and large $\Delta/\gamma$:

$\Gamma_{\rm
rel}\simeq\begin{cases}8|\overline{\alpha}|^{2}e^{-4|\overline{\alpha}|^{2}}\gamma&\text{
if }\Delta\gg\gamma,\\\
16|\overline{\alpha}|^{2}e^{-4|\overline{\alpha}|^{2}}\frac{\Delta^{2}}{\gamma}&\text{if
}\Delta\ll\gamma.\end{cases}$ (35)

As expected, the relaxation rate vanishes in both these extreme limits, as the
relaxation process from $\ket{1}$ to $\ket{2}$ involves both the
$\tilde{\Delta}\propto\Delta$ tunneling process and the dissipative processes
$\Gamma_{2,3}\propto\gamma$. One finds that for fixed $\Delta$ and
$\overline{\alpha}$, $\Gamma_{\rm rel}$ is maximized when the damping rate
$\gamma$ is chosen to be $\gamma=\sqrt{2}\Delta$. This motivates the choice
used in Fig. 3.

Figure 4: Schematic showing the effective slow dynamics in the four
dimensional manifold spanned by even and odd cat-states in each cavity. There
are two dark states of the dissipator in this subspace, denoted $\ket{2}$ and
$\ket{1}$, with $\ket{2}$ the eventual steady state. The detuning $\Delta$ of
the parametric drives induces tunneling between the odd parity states
$\ket{1}$ and $\ket{3}$, while dissipation induces the indicated transitions
between states of different parity. Exponentially small rates are indicates in
dark red.

##  Appendix C Implementation issues

###  C.1 Parametric drives

In the lab frame, the Hamiltonian for two parametrically driven cavities is

$\displaystyle\hat{H}_{\rm Lab}$
$\displaystyle=\omega_{a}\hat{a}^{\dagger}\hat{a}+\lambda(\hat{a}^{\dagger}\hat{a}^{\dagger}e^{-i\omega_{p_{a}}t}+h.c.)$
$\displaystyle+\omega_{b}\hat{b}^{\dagger}\hat{b}-\lambda(\hat{b}^{\dagger}\hat{b}^{\dagger}e^{-i\omega_{p_{b}}t}+h.c.),$
(36)

where $\omega_{a/b}$ are the cavity resonance frequencies and
$\omega_{p_{a/b}}$ are the frequencies of the parametric drives on each
cavity. To obtain the quadratic part of the Hamiltonian in Eq. (20), we go to
a frame rotating at $\omega_{p_{a}}/2$ for cavity $\hat{a}$ and
$\omega_{p_{b}}/2$ for cavity $\hat{b}$, and choose the parametric drive
frequencies such that
$\Delta=\omega_{a}-\omega_{p_{a}}/2=\omega_{p_{b}}/2-\omega_{b}$.

It is most convenient to consider a situation where the parametric drives are
at the same frequency, i.e. $\omega_{p_{a}}=\omega_{p_{b}}\equiv\omega_{p}$,
as this simplifies phase locking of the two drives, which is required to
maintain the relative sign difference in the two drive amplitudes. In this
case, the cavities must have different frequencies, and the parametric drive
frequency is simply taken to be the average cavity frequency:
$2\omega_{p}=\omega_{a}+\omega_{b}$. This results in the detuning parameter
$\Delta$ being given by $\Delta=(\omega_{a}-\omega_{b})/2$.

Choosing identical pump frequencies also has a strong advantage when
considering the correlated dissipator needed for our scheme, as now, this
dissipator can be implemented in a completely passive manner. One can for
example simply couple the two cavities to a waveguide, as discussed in the
next subsection. For this choice of pump frequencies, we can go to an
interaction picture at $\omega_{p}/2$ for the waveguide coupled to the
cavities, with the result that the cavity-waveguide interaction is time
independent for both cavities. As discussed in the next subsection, one can
then derive the needed $\hat{a}+\hat{b}$ dissipator by arranging the distance
between the cavities appropriately.

If instead the cavities are not pumped at the same parametric drive frequency,
the simple waveguide scheme for generating the needed dissipator fails, as we
can no longer define a global interaction frame where both cavity-waveguide
interactions are time-independent. In this case, a correlated dissipator can
still be achieved using frequency conversion via engineered dissipation with
active elements (see e.g. [Sliwa  _et al._(2015)Sliwa, Hatridge, Narla,
Shankar, Frunzio, Schoelkopf, and Devoret, Lecocq  _et al._(2017)Lecocq,
Ranzani, Peterson, Cicak, Simmonds, Teufel, and Aumentado, Lecocq  _et
al._(2016)Lecocq, Clark, Simmonds, Aumentado, and Teufel]).

###  C.2 Non-local dissipator

As discussed in the main text, the non-local dissipator in Eq. (19)
(corresponding to correlated single photon loss) can be generated by coupling
the two cavities to a one-dimensional transmission line or waveguide. A full
derivation of this is provided in Ref. [Chang  _et al._(2012)Chang, Jiang,
Gorshkov, and Kimble] and in Appendix B of Ref. [Metelmann and Clerk(2015)].
In general, coupling two cavities to a waveguide and then eliminating the
waveguide will generate not only a dissipator of the form in Eq. (19), but
also a direct, Hamiltonian tunneling interaction between the cavities. The
relative strength of these two kinds of couplings is a function of the
distance $d$ between the attachment points of the cavities to the waveguide.
For $d=m\tilde{\lambda}$ (where $m$ is an integer, and $\tilde{\lambda}$ is
the wavelength of a waveguide excitation at frequency $\omega_{p}/2$), the
induced Hamiltonian tunneling interaction vanishes, and the only effect is the
desired non-local dissipator with the correct phase. Note that this result
requires that the relevant dissipative dynamics is slow compared to the
propagation time $d/v_{g}$, where $v_{g}$ is the waveguide group velocity. The
slow dynamics of our system (see Appendix B) means that this condition can
easily be satisfied.

###  C.3 Positive $U$

While the use of Josephson junctions in superconducting circuits to induce
effective attractive Kerr nonlinearities is by now routine (see e.g.
[Kirchmair  _et al._(2013)Kirchmair, Vlastakis, Leghtas, Nigg, Paik, Ginossar,
Mirrahimi, Frunzio, Girvin, and Schoelkopf]), it is also possible to generate
the repulsive interaction required in our scheme. One option is to use a qubit
based on an inductively shunted Josephson junction (such as a flux qubit or
fluxonium [Manucharyan  _et al._(2009)Manucharyan, Koch, Glazman, and
Devoret]). The effective Hamiltonian describing such a system has the form
[Girvin(2014)]

$\hat{H}_{\rm
qb}=E_{C}\hat{n}^{2}+\frac{1}{2}E_{L}\hat{\phi}^{2}-E_{J}\cos\left(\hat{\phi}+\phi_{g}\right),$
(37)

where $E_{C}$ is the capacitive charging energy, $E_{L}$ is the inductive
energy associated with the shunt inductance, and $E_{J}$ is the Josephson
energy. $\hat{\phi}$ ($\hat{n}$) is the junction phase (charge) operator. The
parameter $\phi_{g}$ describes a flux bias of the system, and is proportional
to the externally applied magnetic flux. If one applies half a flux quantum,
$\phi_{g}=\pi$, then one effectively changes the sign of the Josephson
potential. Expanding the cosine potential in this case yields a quartic term
with a positive coefficient, corresponding to a repulsive Kerr interaction. If
one now weakly and non-resonantly couples a cavity to this qubit circuit (i.e.
the dispersive regime), the cavity will inherit this repulsive nonlinearity.

An alternate approach would be to use a standard transmon qubit [Koch  _et
al._(2007)Koch, Yu, Gambetta, Houck, Schuster, Majer, Blais, Devoret, Girvin,
and Schoelkopf], but work in the so-called straddling regime, where the cavity
frequency is between the frequencies for the ground to first excited state and
first excited to second excited state transitions of the transmon. In the
dispersive regime, one can treat the interaction between the transmon and the
cavity perturbatively, and standard time-independent perturbation theory to
fourth order [Zhu  _et al._(2013)Zhu, Ferguson, Manucharyan, and Koch] yields
an effective self-Kerr interaction for the cavity given by

$\displaystyle\hat{H}_{\rm Kerr}=U_{\rm
eff}\hat{a}^{\dagger}\hat{a}^{\dagger}\hat{a}\hat{a},\ U_{\rm
eff}=\frac{g^{4}U_{\rm qb}}{2\delta^{3}(\delta+\frac{U_{\rm qb}}{2})}.$ (38)

Here $\delta=\omega_{01}-\omega_{c}$ is the detuning, and $g$ the coupling
between the transmon’s ground to first exited state transition and the cavity,
and $U_{\rm qb}$ is the anharmonicity of the transmon. For a transmon $U_{\rm
qb}\approx-E_{C}$, i.e. $U_{\rm qb}<0$. If $\delta<0$, or $\delta>0$ and
$\delta+U_{\rm qb}/2>0$, then $U_{\rm eff}<0$ and the cavity nonlinearity is
attractive. However, if $\delta>0$ and $\delta+U_{\rm qb}/2<0$, which lies in
the straddling regime (defined by $\delta+U_{\rm qb}<0$), then $U_{\rm eff}>0$
and the cavity will inherit a repulsive nonlinearity from the transmon.

Other, more recent circuit designs also allow for strong repulsive Kerr
nonlinearities, and even the ability to tune the strength of the nonlinearity
_in situ_ [Frattini  _et al._(2017)Frattini, Vool, Shankar, Narla, Sliwa, and
Devoret, Zhang  _et al._(2017)Zhang, Huang, Gershenson, and Bell].

##  Appendix D Wigner Function Negativity

While the fidelity is a good quantitative measure of the success of our scheme
in the presence of internal loss, it would also be interesting to understand
how the quantum correlations in the steady state diminish as a result of
internal loss. To that end, we compare the Wigner functions of the steady
state with no internal loss, and that with finite internal loss rate, where
the lossy steady state (panel (b)) has $90\%$ fidelity with the ideal steady
state.

The relevant information is contained in the negativity of the Wigner
function, and while the negativity is reduced by internal loss, even for a
$10\%$ fidelity loss there is still significant negativity. Further, the
interference fringes of the Wigner function of the non-local mode $\hat{d}$
correspond to the entanglement between the local modes $\hat{a}$ and
$\hat{b}$. As can be seen, the fringes are still clearly visible in the lossy
case, indicating that entanglement between the local modes remains.

Figure 5: Wigner function of the steady state of the non-local mode $\hat{d}$
for (a) zero internal loss (c.f. Eq. (13)) and (b) finite internal loss
$\kappa/U=2\times 10^{-4}$ (identical for both cavities, c.f. Eq. (23)). We
use system parameters of $\lambda/U=2$, $\Delta/U=0.5$, $\gamma/U=1/\sqrt{2}$,
which for (b) corresponds to a fidelity of $F\approx 0.9$, as seen in Fig. 3.
The minimum Wigner function negativity is $W=-0.21$ for zero internal loss in
panel (a), and $W=-0.16$ for finite internal loss in panel (b).

##  Appendix E Parameter mismatch

The pure steady state of our system requires matching of parameters between
two Kerr cavities. In this appendix, we quantify the loss of fidelity from
having the second cavity imperfectly matched with the first. To do so, we
numerically calculate the steady state of the Hamiltonian

$\displaystyle\hat{H}_{\mathrm{mis}}$
$\displaystyle=\Delta\left[\hat{a}^{\dagger}\hat{a}-(1+\delta\Delta)\hat{b}^{\dagger}\hat{b}\right]$
(39)
$\displaystyle+\lambda\left[\hat{a}^{\dagger}\hat{a}^{\dagger}-(1+\delta\lambda)\hat{b}^{\dagger}\hat{b}^{\dagger}+h.c.\right]$
$\displaystyle+U\left[\hat{a}^{\dagger}\hat{a}^{\dagger}\hat{a}\hat{a}-(1+\delta
U)\hat{b}^{\dagger}\hat{b}^{\dagger}\hat{b}\hat{b}\right].$

where $\delta\Delta$, $\delta\lambda$, and $\delta U$ are the dimensionless
mismatch in the cavity detuning, parametric drive strength, and self-Kerr
nonlinearity respectively. The original system is recovered for
$\delta\Delta=\delta\lambda=\delta U=0$. Figure 6 shows the fidelity of the
numerically-computed steady state with the ideal steady state
$\ket{0}_{c}\otimes\Ket{\tilde{C}}_{d}$ of Eq. (21). We find that $F>0.95$ is
maintained for mismatches as large as $5\%$. Having imperfectly matched
cavities thus poses no significant concern for the scheme.

Figure 6: Fidelity $F$ for our system with mismatched (a) cavity detuning
$\Delta$, (b) parametric drive $\lambda$ and (c) self-Kerr nonlinearity $U$,
as described by the modified Hamiltonian in Eq. (39). We compare to the exact
steady state $\ket{0}_{c}\otimes\Ket{\tilde{C}}_{d}$ of Eq. (21). The mismatch
parameters $\delta\Delta$, $\delta\lambda$, $\delta U$ are dimensionless
relative differences between the two cavities. For all curves the parametric
drive strength $\lambda/U=2$, the detuning $\Delta/U=1$, and the correlated
dissipation rate $\gamma/U=2$. As can be seen, $F>0.95$ can be obtained for up
to $5\%$ relative mismatch.

## References

  * [Braunstein and van Loock(2005)] S. L. Braunstein and P. van Loock, Rev. Mod. Phys. 77, 513 (2005).
  * [Poyatos  _et al._(1996)Poyatos, Cirac, and Zoller] J. F. Poyatos, J. I. Cirac, and P. Zoller, Phys. Rev. Lett. 77, 4728 (1996).
  * [Parkins  _et al._(2006)Parkins, Solano, and Cirac] A. S. Parkins, E. Solano, and J. I. Cirac, Phys. Rev. Lett. 96, 053602 (2006).
  * [Krauter  _et al._(2011)Krauter, Muschik, Jensen, Wasilewski, Petersen, Cirac, and Polzik] H. Krauter, C. A. Muschik, K. Jensen, W. Wasilewski, J. M. Petersen, J. I. Cirac, and E. S. Polzik, Phys. Rev. Lett. 107, 080503 (2011).
  * [Muschik  _et al._(2011)Muschik, Polzik, and Cirac] C. A. Muschik, E. S. Polzik, and J. I. Cirac, Phys. Rev. A 83, 052312 (2011).
  * [Woolley and Clerk(2014)] M. J. Woolley and A. A. Clerk, Phys. Rev. A 89, 063805 (2014).
  * [Haroche and Raimond(2006)] S. Haroche and J.-M. Raimond,  _Exploring the Quantum: Atoms, Cavities and Photons_ (Oxford University Press, Oxford, 2006).
  * [Vlastakis  _et al._(2013)Vlastakis, Kirchmair, Leghtas, Nigg, Frunzio, Girvin, Mirrahimi, Devoret, and Schoelkopf] B. Vlastakis, G. Kirchmair, Z. Leghtas, S. E. Nigg, L. Frunzio, S. M. Girvin, M. Mirrahimi, M. H. Devoret, and R. J. Schoelkopf, Science 342, 607 (2013).
  * [Mirrahimi  _et al._(2014)Mirrahimi, Leghtas, Albert, Touzard, Schoelkopf, Jiang, and Devoret] M. Mirrahimi, Z. Leghtas, V. V. Albert, S. Touzard, R. J. Schoelkopf, L. Jiang, and M. H. Devoret, New J. Phys. 16, 045014 (2014).
  * [Wang  _et al._(2016)Wang, Gao, Reinhold, Heeres, Ofek, Chou, Axline, Reagor, Blumoff, Sliwa, Frunzio, Girvin, Jiang, Mirrahimi, Devoret, and Schoelkopf] C. Wang, Y. Y. Gao, P. Reinhold, R. W. Heeres, N. Ofek, K. Chou, C. Axline, M. Reagor, J. Blumoff, K. M. Sliwa, L. Frunzio, S. M. Girvin, L. Jiang, M. Mirrahimi, M. H. Devoret, and R. J. Schoelkopf, Science 352, 1087 (2016).
  * [Sanders(2012)] B. C. Sanders, J. Phys. A: Math. Theor. 45, 244002 (2012).
  * [Leghtas  _et al._(2015)Leghtas, Touzard, Pop, Kou, Vlastakis, Petrenko, Sliwa, Narla, Shankar, Hatridge, Reagor, Frunzio, Schoelkopf, Mirrahimi, and Devoret] Z. Leghtas, S. Touzard, I. M. Pop, A. Kou, B. Vlastakis, A. Petrenko, K. M. Sliwa, A. Narla, S. Shankar, M. J. Hatridge, M. Reagor, L. Frunzio, R. J. Schoelkopf, M. Mirrahimi, and M. H. Devoret, Science 347, 853 (2015).
  * [Sarlette  _et al._(2012)Sarlette, Leghtas, Brune, Raimond, and Rouchon] A. Sarlette, Z. Leghtas, M. Brune, J. M. Raimond, and P. Rouchon, Phys. Rev. A 86, 012114 (2012).
  * [Arenz  _et al._(2013)Arenz, Cormick, Vitali, and Morigi] C. Arenz, C. Cormick, D. Vitali, and G. Morigi, J. Phys. B: At. Mol. Opt. Phys. 46, 224001 (2013).
  * [Pudlik  _et al._(2013)Pudlik, Hennig, Witthaut, and Campbell] T. Pudlik, H. Hennig, D. Witthaut, and D. K. Campbell, Phys. Rev. A 88, 063606 (2013).
  * [Grujic  _et al._(2013)Grujic, Clark, Jaksch, and Angelakis] T. Grujic, S. R. Clark, D. Jaksch, and D. G. Angelakis, Phys. Rev. A 87, 053846 (2013).
  * [Cao  _et al._(2016)Cao, Mahmud, and Hafezi] B. Cao, K. W. Mahmud, and M. Hafezi, Phys. Rev. A 94, 063805 (2016).
  * [Casteels and Ciuti(2017)] W. Casteels and C. Ciuti, Phys. Rev. A 95, 013812 (2017).
  * [Casteels and Wouters(2017)] W. Casteels and M. Wouters, Phys. Rev. A 95, 043833 (2017).
  * [Noh and Angelakis(2017)] C. Noh and D. G. Angelakis, Rep. Prog. Phys 80, 016401 (2017).
  * [Chang  _et al._(2012)Chang, Jiang, Gorshkov, and Kimble] D. E. Chang, L. Jiang, A. V. Gorshkov, and H. J. Kimble, New J. Phys. 14, 063003 (2012).
  * [Metelmann and Clerk(2015)] A. Metelmann and A. A. Clerk, Phys. Rev. X 5, 021025 (2015).
  * [Yamamoto  _et al._(2008)Yamamoto, Inomata, Watanabe, Matsuba, Miyazaki, Oliver, Nakamura, and Tsai] T. Yamamoto, K. Inomata, M. Watanabe, K. Matsuba, T. Miyazaki, W. D. Oliver, Y. Nakamura, and J. S. Tsai, Appl. Phys. Lett. 93, 042510 (2008).
  * [Krantz  _et al._(2016)Krantz, Bengtsson, Simoen, Gustavsson, Shumeiko, Oliver, Wilson, Delsing, and Bylander] P. Krantz, A. Bengtsson, M. Simoen, S. Gustavsson, V. Shumeiko, W. D. Oliver, C. M. Wilson, P. Delsing, and J. Bylander, Nat. Comm. 7, 11417 EP (2016).
  * [Mutus  _et al._(2014)Mutus, White, Barends, Chen, Chen, Chiaro, Dunsworth, Jeffrey, Kelly, Megrant, Neill, O’Malley, Roushan, Sank, Vainsencher, Wenner, Sundqvist, Cleland, and Martinis] J. Y. Mutus, T. C. White, R. Barends, Y. Chen, Z. Chen, B. Chiaro, A. Dunsworth, E. Jeffrey, J. Kelly, A. Megrant, C. Neill, P. J. J. O’Malley, P. Roushan, D. Sank, A. Vainsencher, J. Wenner, K. M. Sundqvist, A. N. Cleland, and J. M. Martinis, Appl. Phys. Lett. 104, 263513 (2014).
  * [Kirchmair  _et al._(2013)Kirchmair, Vlastakis, Leghtas, Nigg, Paik, Ginossar, Mirrahimi, Frunzio, Girvin, and Schoelkopf] G. Kirchmair, B. Vlastakis, Z. Leghtas, S. E. Nigg, H. Paik, E. Ginossar, M. Mirrahimi, L. Frunzio, S. M. Girvin, and R. J. Schoelkopf, Nature 495, 205 EP (2013).
  * [Sliwa  _et al._(2015)Sliwa, Hatridge, Narla, Shankar, Frunzio, Schoelkopf, and Devoret] K. M. Sliwa, M. Hatridge, A. Narla, S. Shankar, L. Frunzio, R. J. Schoelkopf, and M. H. Devoret, Phys. Rev. X 5, 041020 (2015).
  * [Lecocq  _et al._(2017)Lecocq, Ranzani, Peterson, Cicak, Simmonds, Teufel, and Aumentado] F. Lecocq, L. Ranzani, G. A. Peterson, K. Cicak, R. W. Simmonds, J. D. Teufel, and J. Aumentado, Phys. Rev. Applied 7, 024028 (2017).
  * [Lecocq  _et al._(2016)Lecocq, Clark, Simmonds, Aumentado, and Teufel] F. Lecocq, J. B. Clark, R. W. Simmonds, J. Aumentado, and J. D. Teufel, Phys. Rev. Lett. 116, 043601 (2016).
  * [Cirac  _et al._(1993)Cirac, Parkins, Blatt, and Zoller] J. I. Cirac, A. S. Parkins, R. Blatt, and P. Zoller, Phys. Rev. Lett. 70, 556 (1993).
  * [Rabl  _et al._(2004)Rabl, Shnirman, and Zoller] P. Rabl, A. Shnirman, and P. Zoller, Phys. Rev. B 70, 205304 (2004).
  * [Kronwald  _et al._(2013)Kronwald, Marquardt, and Clerk] A. Kronwald, F. Marquardt, and A. A. Clerk, Phys. Rev. A 88, 063833 (2013).
  * [Wollman  _et al._(2015)Wollman, Lei, Weinstein, Suh, Kronwald, Marquardt, Clerk, and Schwab] E. E. Wollman, C. U. Lei, A. J. Weinstein, J. Suh, A. Kronwald, F. Marquardt, A. A. Clerk, and K. C. Schwab, Science 349, 952 (2015).
  * [Lecocq  _et al._(2015)Lecocq, Clark, Simmonds, Aumentado, and Teufel] F. Lecocq, J. B. Clark, R. W. Simmonds, J. Aumentado, and J. D. Teufel, Phys. Rev. X 5, 041037 (2015).
  * [Pirkkalainen  _et al._(2015)Pirkkalainen, Damskägg, Brandt, Massel, and Sillanpää] J.-M. Pirkkalainen, E. Damskägg, M. Brandt, F. Massel, and M. A. Sillanpää, Phys. Rev. Lett. 115, 243601 (2015).
  * [Kienzler  _et al._(2015)Kienzler, Lo, Keitch, de Clercq, Leupold, Lindenfelser, Marinelli, Negnevitsky, and Home] D. Kienzler, H.-Y. Lo, B. Keitch, L. de Clercq, F. Leupold, F. Lindenfelser, M. Marinelli, V. Negnevitsky, and J. P. Home, Science 347, 53 (2015).
  * [Asbóth  _et al._(2004)Asbóth, Adam, Koniorczyk, and Janszky] J. K. Asbóth, P. Adam, M. Koniorczyk, and J. Janszky, Eur. Phys. J. B 30, 403 (2004).
  * [Fang  _et al._(2017)Fang, Luo, Metelmann, Matheny, Marquardt, Clerk, and Painter] K. Fang, J. Luo, A. Metelmann, M. H. Matheny, F. Marquardt, A. A. Clerk, and O. Painter, Nature Physics 13, 465 EP (2017).
  * [Goto(2016a)] H. Goto, Phys. Rev. A 93, 050301 (2016a).
  * [Goto(2016b)] H. Goto, Sci. Rep. 6, 21686 (2016b).
  * [Minganti  _et al._(2016)Minganti, Bartolo, Lolli, Casteels, and Ciuti] F. Minganti, N. Bartolo, J. Lolli, W. Casteels, and C. Ciuti, Sci. Rep. 6, 26987 (2016).
  * [Puri  _et al._(2017)Puri, Boutin, and Blais] S. Puri, S. Boutin, and A. Blais, npj Quant. Inf. 3, 18 (2017).
  * [Drummond and Walls(1980)] P. D. Drummond and D. F. Walls, J. Phys. A: Math. Gen. 13, 725 (1980).
  * [Stannigel  _et al._(2012)Stannigel, Rabl, and Zoller] K. Stannigel, P. Rabl, and P. Zoller, New J. Phys. 14, 063014 (2012).
  * [Reagor  _et al._(2013)Reagor, Paik, Catelani, Sun, Axline, Holland, Pop, Masluk, Brecht, Frunzio, Devoret, Glazman, and Schoelkopf] M. Reagor, H. Paik, G. Catelani, L. Sun, C. Axline, E. Holland, I. M. Pop, N. A. Masluk, T. Brecht, L. Frunzio, M. H. Devoret, L. Glazman, and R. J. Schoelkopf, Appl. Phys. Lett. 102, 192604 (2013).
  * [Reagor  _et al._(2016)Reagor, Pfaff, Axline, Heeres, Ofek, Sliwa, Holland, Wang, Blumoff, Chou, Hatridge, Frunzio, Devoret, Jiang, and Schoelkopf] M. Reagor, W. Pfaff, C. Axline, R. W. Heeres, N. Ofek, K. Sliwa, E. Holland, C. Wang, J. Blumoff, K. Chou, M. J. Hatridge, L. Frunzio, M. H. Devoret, L. Jiang, and R. J. Schoelkopf, Phys. Rev. B 94, 014506 (2016).
  * [Manucharyan  _et al._(2009)Manucharyan, Koch, Glazman, and Devoret] V. E. Manucharyan, J. Koch, L. I. Glazman, and M. H. Devoret, Science 326, 113 (2009).
  * [Girvin(2014)] S. M. Girvin, in  _Quantum Machines: Measurement and Control of Engineered Quantum Systems: Lecture Notes of the Les Houches Summer School_, edited by M. Devoret, R. J. Schoelkopf, B. Huard, and L. F. Cugliandolo (Oxford University Press, 2014) Chap. 3.
  * [Koch  _et al._(2007)Koch, Yu, Gambetta, Houck, Schuster, Majer, Blais, Devoret, Girvin, and Schoelkopf] J. Koch, T. M. Yu, J. Gambetta, A. A. Houck, D. I. Schuster, J. Majer, A. Blais, M. H. Devoret, S. M. Girvin, and R. J. Schoelkopf, Phys. Rev. A 76, 042319 (2007).
  * [Zhu  _et al._(2013)Zhu, Ferguson, Manucharyan, and Koch] G. Zhu, D. G. Ferguson, V. E. Manucharyan, and J. Koch, Phys. Rev. B 87, 024510 (2013).
  * [Frattini  _et al._(2017)Frattini, Vool, Shankar, Narla, Sliwa, and Devoret] N. E. Frattini, U. Vool, S. Shankar, A. Narla, K. M. Sliwa, and M. H. Devoret, Appl. Phys. Lett. 110, 222603 (2017).
  * [Zhang  _et al._(2017)Zhang, Huang, Gershenson, and Bell] W. Zhang, W. Huang, M. E. Gershenson, and M. T. Bell, Phys. Rev. Applied 8, 051001 (2017).


* [Roberts et al., 2017] Roberts, K., Gururaj, A. E., Chen, X., Pournejati, S., Hersh, W. R., Demner-Fushman, D., Ohno-Machado, L., Cohen, T., and Xu, H. (2017).  Information retrieval for biomedical datasets: the 2016 bioCADDIE dataset retrieval challenge.  Database, 2017.
  * [SiriJodha Khalsa, 22z1] SiriJodha Khalsa, Peter Cotroneo, M. W. (2018, doi:10.17632/7j43z6n22z.1).  A survey of current practices in data search services.  Technical report, Research Data Alliance Data (RDA) Discovery Paradigms Interest Group.
  * [Taylor et al., 1411] Taylor, C. F., Field, D., Sansone, S.-A., Aerts, J., Apweiler, R., Ashburner, M., Ball, C. A., Binz, P.-A., Bogue, M., Booth, T., Brazma, A., Brinkman, R. R., Michael Clark, A., Deutsch, E. W., Fiehn, O., Fostel, J., Ghazal, P., Gibson, F., Gray, T., Grimes, G., Hancock, John Mand Hardy, N. W., Hermjakob, H., Julian Jr, R. K., Kane, M., Kettner, C., Kinsinger, C., Kolker, E., Kuiper, M., Novère, N. L., Leebens-Mack, J., Lewis, S. E., Lord, P., Mallon, A.-M., Marthandan, N., Masuya, H., McNally, R., Mehrle, A., Morrison, N., Orchard, S., Quackenbush, J., Reecy, J. M., Robertson, D. G., Rocca-Serra, P., Rodriguez, H., Rosenfelder, H., Santoyo-Lopez, J., Scheuermann, R. H., Schober, D., Smith, B., Snape, J., Stoeckert Jr, C. J., Tipton, K., Sterk, P., Untergasser, A., Vandesompele, J., and Wiemann, S. (2008, https://doi.org/10.1038/nbt.1411).  Promoting coherent minimum reporting guidelines for biological and biomedical investigations: the mibbi project.  Nature Biotechnology, 889.
  * [Thessen et al., 2012] Thessen, A. E., Cui, H., and Mozzherin, D. (2012).  Applications of natural language processing in biodiversity science.  Advances in Bioinformatics, 2012(Article ID 391574):17 pages.
  * [Unger et al., 2014] Unger, C., Freitas, A., and Cimiano, P. (2014).  An Introduction to Question Answering over Linked Data, pages 100–140.  Springer International Publishing, Cham.
  * [W3C, 2012] W3C (2012).  OWL Working Group, OWL 2 Web Ontology Language.  https://www.w3.org/TR/owl2-overview/, accessed on 12.11.2019.
  * [Wilkinson et al., 2016] Wilkinson, M. D., Dumontier, M., Aalbersberg, I. J., Appleton, G., Axton, M., Baak, A., Blomberg, N., Boiten, J.-W., da Silva Santos, L. B., Bourne, P. E., Bouwman, J., Brookes, A. J., Clark, T., Crosas, M., Dillo, I., Dumon, O., Edmunds, S., Evelo, C. T., Finkers, R., Gonzalez-Beltran, A., Gray, A. J., Groth, P., Goble, C., Grethe, J. S., Heringa, J., ’t Hoen, P. A., Hooft, R., Kuhn, T., Kok, R., Kok, J., Lusher, S. J., Martone, M. E., Mons, A., Packer, A. L., Persson, B., Rocca-Serra, P., Roos, M., van Schaik, R., Sansone, S.-A., Schultes, E., Sengstag, T., Slater, T., Strawn, G., Swertz, M. A., Thompson, M., van der Lei, J., van Mulligen, E., Velterop, J., Waagmeester, A., Wittenburg, P., Wolstencroft, K., Zhao, J., and Mons, B. (2016).  The fair guiding principles for scientific data management and stewardship.  Scientific Data 3, (160018).
  * [Zenodo, 2019a] Zenodo (2019a).  https://zenodo.org/, accessed on 16.05.2019.
  * [Zenodo, 2019b] Zenodo (2019b).  Search.  https://zenodo.org/api/records, accessed on 30.11.2019.