Dissipative stabilization of entangled cat states using a driven Bose-Hubbard dimer
M. Mamaev1 , L. C. G. Govia2 , and A. A. Clerk2
1 Department of Physics, McGill University, Montréal, Québec, Canada.

arXiv:1711.06662v2 [quant-ph] 23 Mar 2018

2 Institute for Molecular Engineering, University of Chicago, 5640 S. Ellis Ave., Chicago, IL 60637

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 drivendissipative system, and use it show that the system permits the preparation and stabilizationof pure entangled non-Gaussian states, so-called
entangled cat states. Unlike previous proposals
for dissipative stabilization of such states, ourapproach requires only a linear coupling to asingle 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 [1], it is necessary to prepare andideally stabilize non-trivial entangled states of multiple bosonic modes. An ideal approach here is reservoir engineering [2], where controlled interactions with
tailored dissipation are used to accomplish both thesegoals. The simplest kind of entangled bosonic states
are the two-mode squeezed states, and there exist manyclosely-related reservoir engineering protocols which allow their stabilization (see e.g. [3–6]).
While useful in certain applications, two-mode
squeezed states are Gaussian, and hence describable bya 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 Schrödinger cat states [7]. These states are integral
to bosonic codes for quantum computing [8–10], andare equivalent to so-called “entangled coherent states”,
which have applications in many aspects of quantum
information [11].
Accepted in Quantum 2018-03-23, click title to verify

a

λ

b

λ

γ

γ

a

b

a

b

Reservoir

Figure 1: Left panel: Schematic of the system: two bosonicmodes â, b̂ are subjected to two-photon driving (amplitude λ),and have on-site Kerr-interactions which are repulsive(attractive) for mode â (mode b̂). The two modes couple to asingle engineered Markovian reservoir (e.g. a waveguide) whichmediates the dissipative equivalent of a tunneling interaction.
For appropriately matched parameters (see main text), thesystem has a unique pure steady state corresponding to an
entangled Schrödinger cat state. This state is an equalsuperposition 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.

Dissipative stabilization of non-Gaussian states is inprinciple 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 asingle mode: for example, a recently implemented protocol requires one to engineer a two-photon loss operator [9, 12]. There exist protocols for stabilizing entangled cat states [13, 14] 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.
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 wellstudied dissipative Bose-Hubbard dimer model (seee.g. [15–20]), where two driven bosonic modes with onsite Kerr interactions are coupled via coherent tunneling. We make some crucial modifications to this stan1

dard 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), whichmediates the dissipative equivalent of a tunnel interaction. This dissipative coupling could be generated inmany different ways, including by simply coupling bothcavities passively to a low-loss waveguide [21, 22]. 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, somethingthat is not possible for the standard driven-dissipative
Bose Hubbard dimer [17]. Using our exact solution, we
identify a simple parameter regime where the steadystate is a pure entangled cat state. Our scheme is alsoeasily modified so that it stabilizes a two-dimensional
manifold spanned by two even-parity cat states; thisallows a potential two-cavity generalization of the catstate encoded qubit introduced in Ref. [9]. We stressthat all the key ingredients required for our scheme
are available in state of the art circuit QED setups(namely parametric drives, see e.g. [23–25], strong onsite Kerr interactions [26] and engineered dissipation[27–29]). 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 wellunderstood reservoir engineering protocol for stabilizing
a single-mode squeezed state [30–32]. This protocol wasrecently realized experimentally both in optomechanics[33–35] as well as in a trapped ion system [36]. The
protocol is based on the observation that a squeezed
state is the vacuum of a bosonic Bogoliubov operatorβ̂[r] = cosh(r)dˆ + sinh(r)dˆ† , where dˆ is the annihilation operator for the mode of interest. Hence, by simply cooling the mode β̂[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 ρ̂ of the mode:


1 †
d
†
ρ̂ = Γ ẑ ρ̂ẑ − {ẑ ẑ, ρ̂} ≡ ΓL[ẑ]ρ̂,(1)
dt2
Accepted in Quantum 2018-03-23, click title to verify

with a jump operator ẑ given byẑ ≡ µdˆ + ν dˆ† .

(2)

Γ represents the rate of the dissipative interaction(i.e. the cooling rate), and the parameters µ, ν are settoµ = cosh r, ν = sinh r,(3)making the jump operator ẑ equal to the Bogoliubovmode annihilation operator β̂[r]. The unique steadystate of this master equation is then the dark state of
the jump operator ẑ, which is nothing else than thevacuum squeezed state with squeeze parameter r.

2.2 Dissipative stabilization of a single-mode catstate
To modify this scheme to stabilize a non-Gaussian state,
we now introduce nonlinearity into the jump operatorẑ. One simple way to do this is to make the parameterµ a photon-number dependent operator. We take µ =ˆ and consider a master equation of the form inµ1 dˆ† d,
Eq. (1) with jump operatorẑ = µ1 dˆ† dˆ dˆ + ν dˆ† .
(4)
The steady states of this master equation correspondto the dark states of the jump operator ẑ. Despite thenonlinearity, these dark states are easily found: theycorrespond to a two-dimensional subspace spanned byeven and odd Schrödinger cat states. Defining even/odd
cat-states in terms of coherent states |αi as|C± (α)i = p|αi ± |−αi2(1 ± e−2|α|2 ),(5)

we find that our jump operator has dark states√ẑ |C± ( 2α)i = 0,

(6)where the amplitude α is given byp
α = i ν/(2µ1 ).

(7)

For some applications, the ability to stabilize a manifold of dark-states can be extremely useful, see e.g. [9].
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 andodd cat-states above, and stabilize, e.g., only the evenparity 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ẑ = µ0 + µ1 dˆ† dˆ dˆ + ν dˆ† ≡ ẑcat .
(8)2

As we show explicitly in Appendix A, for µ0 6= 0, theonly possible dark states of this dissipator are evenparity states. We thus find a unique dark state (andhence pure steady state of Eq. (1)):
ẑcat C̃[µ0 , µ1 , ν] = 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 thatin the limit of small µ0 , this state asymptotically approaches an even-parity cat state:
D
√ E
(10)
lim
C̃[µ0 , µ1 , ν] C+ [ 2α] → 1.
µ0 /µ1 →0

Even for small but finite µ0 /µ1 , |C̃[µ0 , µ1 , ν]i 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 asingle-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 ĉ which couples nonlinearly to the principle mode dˆ of interest. To that end,
we consider a two-mode system whose non-dissipative
dynamics is described by the Hamiltonian†Ĥ1 = Λ ĉ† ẑcat + ẑcatĉhiˆ dˆ + ν dˆ† + h.c.
= Λĉ† µ0 dˆ + µ1 (dˆ† d)(11)
The quadratic terms in this Hamiltonian correspond toa simple beam-splitter interaction (∝ µ0 ) and two-mode
squeezing interaction (∝ ν). The remaining nonlinear
term (∝ µ1 ) describes photon-number dependent tunneling.
We now add damping of the ĉ mode (loss at rate γ),such that the full evolution of the two modes is givenbydρ̂ = −i[Ĥ1 , ρ̂] + γL[ĉ]ρ̂.
(12)dt
In the limit where the ĉ mode is heavily damped(i.e. γ  Λ), standard adiabatic elimination yields a
master equation of the form in Eq. (1) with jump operator ẑcat and with Γ ∝ Λ2 /γ. In this limit, our two-mode
system has a unique steady state analogous to that in
Eq. (9):|ψss i = |0ic |C̃[µ0 , µ1 , ν]id .
Accepted in Quantum 2018-03-23, click title to verify(13)

As discussed earlier, the state |C̃[µ0 , µ1 , ν]id asymptotically approaches an even cat-state in the limit µ0 /µ1 →0.
We next make a crucial observation: even if one isfar 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 Ĥ1 and ĉ bothannihilate 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 γ, Λ, µ0 , µ1 and ν.

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 BoseHubbard dimer model. Surprisingly, we will find thatthe resulting setup is much simpler than that describedin Eq. (11).
ˆ withWe again start with our two modes ĉ and d,
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 â and b̂:1ˆ b̂ = √1 (ĉ − d).
ˆâ = √ (ĉ + d),
2
2

(14)

The dissipative dynamics described in Eq. (12) stabilizes the state |ψss i defined in Eq. (13). Transforming
this state to the “localized” â, b̂ mode basis is equivalent
to passing the state through a beam splitter. It is wellknown that two-mode entangled coherent states can begenerated by sending single-mode cat states into thearms of a 50/50 beam-splitter [37]. We thus find thatin the localized basis (and for µ0  µ1 ), our steadystate has the form of a two-mode entangled cat state:1|ψss i ≈ √ [|C+ (α)ia |C+ (α)ib − |C− (α)ia |C− (α)ib ] .
2(15)
To determine the feasibility of our scheme, we needto express the Hamiltonian Ĥ1 in terms of the localized
mode operators â, b̂. Before doing this, we first modify
Ĥ1 to make it more symmetric, but without changing
the steady state:Ĥ2 = Ĥ1 + Λµ1 ĉ† ĉ† ĉ dˆ + h.c.
hiˆ dˆ + ν dˆ† + h.c. (16)= Λĉ† µ0 dˆ + µ1 (ĉ† ĉ + dˆ† d)

3

4 Exact solution for the steady state

We again consider the master equationdρ̂ = −i[Ĥ2 , ρ̂] + γL[ĉ]ρ̂.
dt

(17)

It is easy to see that the state |ψss i (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 Ĥ2 in the localized mode basis. Defining forconvenience∆ = Λµ0 , λ = Λν/2, U = Λµ1 ,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 ĉ, dˆ mode basis. In this basis
the steady state |ψss i has the ĉ mode in vacuum, while
the 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 dˆ modeis given by:C̃[µ0 , µ1 , ν] d = Nc̃2n =dρ̂ = −i[Ĥ2 , ρ̂] + γL[â + b̂]ρ̂,
dt

(19)

with
Ĥ2 = ∆â† â + λ(â† â† + h.c.) + U â† â† ââ− ∆b̂† b̂ − λ(b̂† b̂† + h.c.) − U b̂† b̂† b̂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. Eachcavity has a Kerr nonlinearity (strength U, −U ) andis subject to a parametric drive (drive detunings ±∆,drive amplitudes ±λ). 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. [38]) in many different contexts. For thesimple case where each cavity has the same resonancefrequency, it can be generated by coupling both cavities to a transmission line with an appropriately chosendistance between the cavities [21, 22]. Further details
about implementation issues (including how to realize arepulsive 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 emergenceof cat states and their adiabatic preparation in a single
parametrically-driven Kerr cavity is the subject of several recent works [39–42]; our scheme leverages similar
resources in two cavities (along with engineered dissipation) to now stabilize an entangled cat state. Notethat the amplitude α of our cat state (c.f. Eq. (7)) intermspof the energy scales in the Hamiltonian is now
α = i λ/U .
Accepted in Quantum 2018-03-23, click title to verifyc̃2n |2nid ,

(21)

n=0

(18)

we obtain

∞
X

−ν
2µ1

n s



µ0
1
Γ
+
2
2µ1
(2n − 1)!!
.
(2n)!! Γ n + 1 + µ0
2
2µ1

(22)

Here |2nid denotes a Fock state with 2n photons, Γ(z) is
the Gamma function, and N is an overall normalization
constant. Our ability to easily find an analytic form
relies on the symmetry of the problem and the existenceof a pure steady state; it is much simpler than moregeneral methods based on the positive-P function [17,41, 43]. Our solution method is reminiscent of that used
in Ref. [44], which studied dissipative entanglement incascaded quantum systems; in contrast to that work,
our approach does not require cascaded interactions, oreven any breaking of time-reversal symmetry.
The analytic expression for |ψss i allows us to confirm
the expected parameter dependence of the state. Forµ0  µ1 , nonlinearity plays a minor role, and the state
tends to a Gaussian squeezed state. In contrast, forµ0  µ1 , the nonlinearity is crucial to the existenceof a steady state, and |ψss i tends to the desired evencat-state. Fig. 2 uses the exact expression of the steadystate |C̃id to show this evolution as a function of the ratio ∆/U = µ0 /µ1 . The upper panel depicts the Wigner
function for the dˆ mode at three specific values of ∆/U .
As described in the figure caption, this describes a cut
through the two-mode Wigner function describing thestate in the localized â, b̂ basis. One clearly sees the
cross-over from an even cat-state to a squeezed state as∆ is increased. This figure also demonstrates that evenfor modest values of ∆/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. To4

pa -pb

0.3
0.2
0.1
0
-0.1
-0.2
-0.3

x a -x b(d) Fidelity, F

1.

order of 10-100 Hz [45, 46].
Note that for the results in Fig. (3), we have taken√ thecoupling γ to the engineered dissipation to be γ = 2∆.
As discussed in more detail in Appendix B, this choiceminimizes the timescale required to prepare the steadystate, and thus helps minimize the effects of internal
loss. However, the fidelity is not extremely sensitive tothis precise tuning of the value of γ.
Fidelity, F
1

0.8
0.6
0.4

λ/U = 20.2

λ/U = 3

1

2

0.9
0.8
3

4

5

6

7

8

Δ/U

Figure 2: (a-c) Wigner functions of the steady state(c.f. Eq. (13)) for the non-local mode dˆ at snapshots of ∆/U .
This is equivalent to a cut through the two-mode Wigner
function for the local modes,
defined by the√√ combined
quadratures (x̂a − x̂b )/ 2 and (p̂a − p̂b )/ 2, where x̂k and p̂k
are the usual quadrature operators for mode k. As ∆/U
increases, the evolution from an even-parity cat state to asqueezed state is clearly visible. (d) Fidelity of the exact steadystate (c.f. Eq. (21)) with both the squeezed state (dashed lines)
obtained for ∆/U → ∞, and the even cat-state (solid lines)
obtained for ∆/U → 0 (c.f. Eq. (15)). The fidelity is defined by
F (|ψ1 i , |ψ2 i) = |hψ1 | ψ2 i|2 .

assess the impact of such losses, we add single-photon
loss (at a rate κ in each cavity) to our master equation:dρ̂ = −i[Ĥ2 , ρ̂] + γL[â + b̂]ρ + κL[â]ρ̂ + κL[b̂]ρ̂. (23)dtWith the inclusion of single particle loss, the steadystate will no longer be pure, and we cannot find ananalytic closed form solution. We instead numericallysimulate the master equation. At a heuristic level, thesensitivity 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, this2timescale can be extremely long, scaling as e4|α| . As aresult, low internal loss levels (relative to U , λ, and ∆)are required for good performance.
Fig. 3 plots the fidelity of the numerically-obtained
steady state with both the ideal, loss-free steady state|ψss i (c.f. Eq. (13)), as well as with the ideal two-mode
entangled cat-state of Eq. (15), as a function of theinternal loss rate κ. We find that for ∆/U = 1, we get
F ≈ 0.95 with the exact solution and F ≈ 0.86 with theideal cat state for a loss rate κ = 10−3 U . This valueof 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 ofU ∼ 100-1000 kHz has been realized [26], while in stateof-the-art 3D cavities internal leakage rates are on the
Accepted in Quantum 2018-03-23, click title to verify

0.7
0.6

Δ/U = 0.1

0.5

Δ/U = 0.5
Δ/U = 1

0.4
10-5

10-4

10-3

10-2

κ/U

Figure 3: Effect of intrinsic cavity losses κ (identical for bothcavities, c.f. Eq. (23)) on the steady state, for various values ofthe parametric drive detuning parameter ∆. Solid lines are thefidelity 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 ∆/U → 0).
For all curves the parametric drive√ strength is set to λ/U = 2(thus setting the amplitude α = i 2), and the coupling
rate γ√to the correlated dissipation is taken to be γ = 2∆ tomaximize the speed of internal system dynamics (see Appendix
p√
√  2
B). The fidelity is defined by F (ρ1 , ρ2 ) = Trρ1 ρ2 ρ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 themagnitude 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 matchingof the Hamiltonian parameters of the two cavities, thedesired behaviour is still found for small mismatches.
For the same parameter choices as in Fig. 3, a driveor 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 thatuses a modified version of a driven-dissipative Bose
Hubbard dimer to stabilize a two-cavity entangled cat
5state. The scheme uses detuned parametric drives applied to each cavity, and couples them only through alinear dissipative tunneling interaction that is mediatedby an external reservoir. While we have focused on thecase where there is a unique steady state, by setting
the parametric drive detunings ∆ = 0, one would instead stabilize a two dimensional manifold spanned by
the entangled cat states (the states |2i and |1i definedin 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
A Analytic expression for the steadystate
We derive here the exact steady state of the systemdescribed by the Hamiltonian Ĥ2 in Eq. (16) and themaster equation in Eq. (17). Motivated by the fact that
only the ĉ mode experiences any damping, we makean ansatz that there exists a pure-state steady state inwhich the ĉ mode is in vacuum, i.e.
|ψi = |0ic ⊗∞X

αn |nid .

(24)

n=0

Here, αn are arbitrary normalized coefficients and |nid
are Fock states.
The above state is trivially a dark-state of the jump
operator ĉ appearing in the dissipator of the masterequation in Eq. (17). To be a steady state, it must thennecessarily also be an eigenstate of the Hamiltonian.
Acting on |ψi with Ĥ2 yields,
Ĥ2 |ψi = |1ic ⊗

∞
X
n=0

αn [∆ + U (n − 1)]√

n |n − 1id√+ 2λ n + 1 |n + 1id .

(25)

where we have expressed µ0 , µ1 and ν in terms of ∆, U
and λ using Eq. (18).
The RHS of Eq. (25) has the ĉ mode in the state|1i, while in |ψi it is in the state |0i. Thus, the only
way that |ψi can be an eigenstate of Ĥ2 is if it has
zero energy. We thus need Ĥ2 to annihilate |ψi, andtherefore require that the coefficient of each |nid state
Accepted in Quantum 2018-03-23, click title to verify

in Eq. (25) vanish. This leads to the recursion relation(0(∆ + U n) n + 1αn+1 =√−2λ nαn−1√

n = 0,n ≥ 1.
(26)

The above recursion relation does not mix the evenand odd-parity subspaces. Note that for ∆ 6= 0, then = 0 case of Eq. (22) forces α1 = 0, and as a result, forces all odd αj to be zero. The solution thus
has only even-parity Fock states. Solving Eq. (26) for
the remaining even parity coefficients directly gives thesolution presented in Eq. (22), with α2n ≡ c̃n .
Note also that if ∆ = 0, then the n = 0 case of
Eq. (22) is trivially satisfied regardless of the value of
α1 , and one can find both even and odd parity steadystates. Thus, as discussed in the main text, a non-zerodetuning ∆ is crucial in order to obtain a unique steadystate.
While the above discussion has identified a puresteady state for our system when ∆ 6= 0, the question
remains whether this is a unique steady state. In thecase where U = 0, the system is linear and one canexactly solve the system. The linear system is stablefor λ ≤ ∆, and in this regime the unique steady stateis exactly that given by Eq. (22). For non-zero U , wehave no rigorous analytic argument precluding the existence of multiple steady states (i.e. in addition to thestate we describe analytically here). We have however
investigated the model numerically over a wide range ofparameters, and find that as long as ∆ is non-zero, thesteady state described by Eq. (21) is indeed the uniquesteady state.

B Dynamical timescales
We consider in this appendix the parameter dependenceof 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, thistimescale sets the sensitivity to unwanted perturbations
such as internal loss. Note first that when the drive detuning ∆ = 0, each cavity has exact zero-energy eigenstates |C± (α)i [39, 40, 42]. For the interesting limit ofa large cat state amplitude α, the slowest timescalesin our system correspond to dynamics within the four
dimensional subspace spanned by these cat states. Entangled cat states form an orthonormal basis for this6

subspace:EVEN

1|2̃i = √ (|C+ (α)i |C+ (α)i − |C− (α)i |C− (α)i) ,2
1|4̃i = √ (|C+ (α)i |C+ (α)i + |C− (α)i |C− (α)i) ,2
1|1i = √ (|C+ (α)i |C− (α)i − |C− (α)i |C+ (α)i) ,2
1|3i = √ (|C+ (α)i |C− (α)i + |C− (α)i |C+ (α)i) .
2

ODD
3,4

|4i

|3i

˜

4,3
2,3

|2i

|1i

(27)

The states |2̃i , |4̃i have even total photon number parity, while the states |1i , |3i have an odd total photon
number parity.
It is useful to make a slight rotation of the even parity
states:|2i ≡ cos θ |2̃i + sin θ |4̃i ,(28)|4i ≡ − sin θ |2̃i + cos θ |4̃i .

(29)

Figure 4: Schematic showing the effective slow dynamics in thefour dimensional manifold spanned by even and odd cat-statesin each cavity. There are two dark states of the dissipator in thissubspace, denoted |2i and |1i, with |2i the eventual steady state.
The detuning ∆ of the parametric drives induces tunnelingbetween the odd parity states |1i and |3i, while dissipationinduces the indicated transitions between states of different
parity. Exponentially small rates are indicates in dark red.

states |1i and |2i is described bywhere for large |α| we have2

sin θ ' 2e−2|α| .

(30)

Note that in the large |α| limit, |2i → |2̃i, which is our
expected entangled-cat steady state.
The states |1i , |2i , |3i , |4i will form a useful basis forour cat state subspace. One finds that there are two distinct dark states in this subspace: they correspond tothe basis states |1i and |2i. If we now restrict our masterequation to this cat-state subspace, we find the effective
dynamics shown schematically in Fig. 4. The detuning
term in Ĥ2 becomes a coherent tunneling term of am˜ between the odd parity states |1i and |3i,
plitude ∆whereas the dissipator gives rise to the three depictedincoherent transitions (with rates Γ3,4 , Γ4,3 , Γ2,3 ). One
sees clearly from this schematic that, as expected, thedark state |2i will be our eventual steady state.
˜ and Γ2,3 become exFor large |α|, one finds that ∆ponentially small:˜ ' −4|α|2 e−2|α|2 ∆,∆2

Γ2,3 ' 16|α|2 e−4|α| γ.

(31)2Γ∆ + Γ2,3 −q

2(2Γ∆ ) + (Γ2,3 )

2


,

(33)

where we have defined the effective incoherent tunnelingrate between the state |1i and the {|3i, |4i} bright-state
manifold as:

Γ∆ =

2
˜2
2 ∆
4∆
.
' 16|α|2 e−4|α|Γ2,3 + Γ4,3γ

(34)

The rate Γrel describes an effective two step process
where population first incoherently tunnels from |1i
to the bright state manifold (rate Γ∆ ), and then relaxes to the steady state |2i (rate Γ2,3 ). Note that
−1 −1Γrel 6= Γ−1, as the incoherent tunneling pro∆ + Γ2,3
cess described by Γ∆ can also move population thewrong way, i.e. from |3i back to the state |1i.
It is instructive to consider the behaviour of Γrel forsmall and large ∆/γ:(32)

The exponential suppression of these rates corresponds
to the extremely weak overlap of the coherent states|αi and |−αi. The above processes will then form thebottleneck in having populations relax to the state |2i,
and will determine the slow timescales characterizing
state preparation.
The superoperator describing this reduced masterequation can be readily diagonalized, yielding the relevant dynamical rates. In the large α limit, we findthat the rate-limiting slow relaxation between the dark
Accepted in Quantum 2018-03-23, click title to verify

1Γrel =
2(
Γrel '

2

8|α|2 e−4|α| γ22
16|α|2 e−4|α| ∆γ

if ∆  γ,if ∆  γ.

(35)

As expected, the relaxation rate vanishes in both these
extreme limits, as the relaxation process from |1i to˜ ∝ ∆ tunneling process and the|2i involves both the ∆dissipative processes Γ2,3 ∝ γ. One finds that for fixed∆ and α, Γrel is maximized
when the damping rate γ is√chosen to be γ = 2∆. This motivates the choice usedin Fig. 3.
7

C Implementation issues
C.1 Parametric drives
In the lab frame, the Hamiltonian for two parametrically driven cavities is
ĤLab = ωa â† â + λ(â† â† e−iωpa t + h.c.)+ ωb b̂† b̂ − λ(b̂† b̂† e−iωpb t + h.c.),(36)

where ωa/b are the cavity resonance frequencies andωpa/b are the frequencies of the parametric drives oneach cavity. To obtain the quadratic part of the Hamiltonian in Eq. (20), we go to a frame rotating at ωpa /2
for cavity â and ωpb /2 for cavity b̂, and choose the parametric drive frequencies such that ∆ = ωa − ωpa /2 =
ωpb /2 − ωb .
It is most convenient to consider a situation where theparametric drives are at the same frequency, i.e. ωpa =ωpb ≡ ωp , as this simplifies phase locking of the two
drives, which is required to maintain the relative signdifference in the two drive amplitudes. In this case, thecavities must have different frequencies, and the parametric drive frequency is simply taken to be the average
cavity frequency: 2ωp = ωa + ωb . This results in the detuning parameter ∆ being given by ∆ = (ωa − ωb )/2.
Choosing identical pump frequencies also has a strong
advantage when considering the correlated dissipatorneeded for our scheme, as now, this dissipator can be
implemented in a completely passive manner. One canfor example simply couple the two cavities to a waveguide, as discussed in the next subsection. For this choiceof pump frequencies, we can go to an interaction pictureat ωp /2 for the waveguide coupled to the cavities, withthe result that the cavity-waveguide interaction is time
independent for both cavities. As discussed in the next
subsection, one can then derive the needed â + b̂ dissipator by arranging the distance between the cavities
appropriately.
If instead the cavities are not pumped at thesame 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 timeindependent. In this case, a correlated dissipator canstill be achieved using frequency conversion via engineered dissipation with active elements (see e.g. [27–29]).

C.2 Non-local dissipator
As discussed in the main text, the non-local dissipatorin Eq. (19) (corresponding to correlated single photonloss) can be generated by coupling the two cavities toa one-dimensional transmission line or waveguide. A
Accepted in Quantum 2018-03-23, click title to verify

full derivation of this is provided in Ref. [21] and in
Appendix B of Ref. [22]. 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 ofthese two kinds of couplings is a function of the distanced between the attachment points of the cavities to the
waveguide. For d = mλ̃ (where m is an integer, and λ̃ isthe wavelength of a waveguide excitation at frequency
ω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/vg , where vg is the
waveguide group velocity. The slow dynamics of oursystem (see Appendix B) means that this condition caneasily 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. [26]), it is also possible to generate the repulsive interaction required inour scheme. One option is to use a qubit based on an
inductively shunted Josephson junction (such as a flux
qubit or fluxonium [47]). The effective Hamiltonian describing such a system has the form [48]1Ĥqb = EC n̂2 + EL φ̂2 − EJ cos φ̂ + φg ,2(37)

where EC is the capacitive charging energy, EL is theinductive energy associated with the shunt inductance,and EJ is the Josephson energy. φ̂ (n̂) is the junction
phase (charge) operator. The parameter φg describes aflux bias of the system, and is proportional to the externally applied magnetic flux. If one applies half a flux
quantum, φg = π, then one effectively changes the signof the Josephson potential. Expanding the cosine potential in this case yields a quartic term with a positivecoefficient, corresponding to a repulsive Kerr interaction. If one now weakly and non-resonantly couples acavity to this qubit circuit (i.e. the dispersive regime),
the cavity will inherit this repulsive nonlinearity.
An alternate approach would be to use a standardtransmon qubit [49], but work in the so-called straddling regime, where the cavity frequency is between thefrequencies for the ground to first excited state and firstexcited 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 [50] yields an effective self-Kerr8interaction for the cavity given by
g 4 Uqb2δ 3 (δ +

Uqb
2 )

.

(38)

Here δ = ω01 −ωc is the detuning, and g the coupling between the transmon’s ground to first exited state transition and the cavity, and Uqb is the anharmonicity of thetransmon. For a transmon Uqb ≈ −EC , i.e. Uqb < 0. If
δ < 0, or δ > 0 and δ + Uqb /2 > 0, then Ueff < 0 andthe cavity nonlinearity is attractive. However, if δ > 0
and δ + Uqb /2 < 0, which lies in the straddling regime(defined by δ + Uqb < 0), then Ueff > 0 and the cavity
will inherit a repulsive nonlinearity from the transmon.
Other, more recent circuit designs also allow forstrong repulsive Kerr nonlinearities, and even the ability
to tune the strength of the nonlinearity in situ [51, 52].

D Wigner Function Negativity
While the fidelity is a good quantitative measure of thesuccess of our scheme in the presence of internal loss, itwould also be interesting to understand how the quantum correlations in the steady state diminish as a resultof internal loss. To that end, we compare the Wigner
functions of the steady state with no internal loss, andthat with finite internal loss rate, where the lossy steadystate (panel (b)) has 90% fidelity with the ideal steadystate.
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 dˆ correspond to the entanglement between the local modes â and b̂. As can be seen, the fringes are stillclearly visible in the lossy case, indicating that entanglement between the local modes remains.

E Parameter mismatch
The pure steady state of our system requires matchingof parameters between two Kerr cavities. In this appendix, we quantify the loss of fidelity from having thesecond cavity imperfectly matched with the first. To
do so, we numerically calculate the steady state of the
HamiltonianhiĤmis = ∆ â† â − (1 + δ∆)b̂† b̂
hi+ λ â† â† − (1 + δλ)b̂† b̂† + h.c.
(39)hi
+ U â† â† ââ − (1 + δU )b̂† b̂† b̂b̂ .
Accepted in Quantum 2018-03-23, click title to verify

pa -pbĤKerr = Ueff â† â† ââ, Ueff =

0.3
0.2
0.1
0
-0.1
-0.2
-0.3

x a -x b
Figure 5: Wigner function of the steady state of the non-local
mode dˆ for (a) zero internal loss (c.f. Eq. (13)) and (b) finite
internal loss κ/U = 2 × 10−4 (identical for both cavities,
c.f. Eq. (23)). We use √system parameters of λ/U = 2,∆/U = 0.5, γ/U = 1/ 2, which for (b) corresponds to a fidelityof F ≈ 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).

where δ∆, δλ, and δU are the dimensionless mismatchin the cavity detuning, parametric drive strength, andself-Kerr nonlinearity respectively. The original system
is recovered for δ∆ = δλ = δU = 0. Figure 6 shows the
fidelity of the numerically-computed steady state with
the ideal steady state |0ic ⊗ C̃ d of Eq. (21). We findthat F > 0.95 is maintained for mismatches as large as
5%. Having imperfectly matched cavities thus poses nosignificant concern for the scheme.
F
1.

(a)

-0.1

-0.05

F
1.

(b)

F
1.

(c)

0.98

0.98

0.98

0.96

0.96

0.96

0.94

0.94

0.

0.05

0.1

δΔ

-0.1

-0.05

0.

0.94
0.05

0.1

δλ

-0.1

-0.05

0.

0.05

0.1

δU

Figure 6: Fidelity F for our system with mismatched (a) cavity
detuning ∆, (b) parametric drive λ and (c) self-Kerr nonlinearity
U , as described by the modified Hamiltonian in Eq. (39). We
compare to the exact steady state |0ic ⊗ C̃of Eq. (21). Thedmismatch parameters δ∆, δλ, δU are dimensionless relative
differences between the two cavities. For all curves theparametric drive strength λ/U = 2, the detuning ∆/U = 1, andthe correlated dissipation rate γ/U = 2. As can be seen,
F > 0.95 can be obtained for up to 5% relative mismatch.

References[1] S. L. Braunstein and P. van Loock, Rev. Mod.
Phys. 77, 513 (2005).
[2] J. F. Poyatos, J. I. Cirac, and P. Zoller, Phys. Rev.
Lett. 77, 4728 (1996).
[3] A. S. Parkins, E. Solano, and J. I. Cirac, Phys.
Rev. Lett. 96, 053602 (2006).
[4] H. Krauter, C. A. Muschik, K. Jensen,
W. Wasilewski, J. M. Petersen, J. I. Cirac,
9and E. S. Polzik, Phys. Rev. Lett. 107, 080503(2011).
[5] C. A. Muschik, E. S. Polzik, and J. I. Cirac, Phys.
Rev. A 83, 052312 (2011).
[6] M. J. Woolley and A. A. Clerk, Phys. Rev. A 89,
063805 (2014).
[7] S. Haroche and J.-M. Raimond, Exploring the
Quantum: Atoms, Cavities and Photons (Oxford
University Press, Oxford, 2006).
[8] 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).
[9] M. Mirrahimi, Z. Leghtas, V. V. Albert,
S. Touzard, R. J. Schoelkopf, L. Jiang, and M. H.
Devoret, New J. Phys. 16, 045014 (2014).
[10] 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).
[11] B. C. Sanders, J. Phys. A: Math. Theor. 45, 244002(2012).
[12] 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).
[13] A. Sarlette, Z. Leghtas, M. Brune, J. M. Raimond,and P. Rouchon, Phys. Rev. A 86, 012114 (2012).
[14] C. Arenz, C. Cormick, D. Vitali, and G. Morigi,
J. Phys. B: At. Mol. Opt. Phys. 46, 224001 (2013).
[15] T. Pudlik, H. Hennig, D. Witthaut, and D. K.
Campbell, Phys. Rev. A 88, 063606 (2013).
[16] T. Grujic, S. R. Clark, D. Jaksch, and D. G. Angelakis, Phys. Rev. A 87, 053846 (2013).
[17] B. Cao, K. W. Mahmud, and M. Hafezi, Phys.
Rev. A 94, 063805 (2016).
[18] W. Casteels and C. Ciuti, Phys. Rev. A 95, 013812(2017).
[19] W. Casteels and M. Wouters, Phys. Rev. A 95,
043833 (2017).
[20] C. Noh and D. G. Angelakis, Rep. Prog. Phys 80,
016401 (2017).
[21] D. E. Chang, L. Jiang, A. V. Gorshkov, and H. J.
Kimble, New J. Phys. 14, 063003 (2012).
[22] A. Metelmann and A. A. Clerk, Phys. Rev. X 5,021025 (2015).
[23] 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).
[24] P. Krantz, A. Bengtsson, M. Simoen, S. Gustavsson, V. Shumeiko, W. D. Oliver, C. M. Wilson,
Accepted in Quantum 2018-03-23, click title to verify

P. Delsing, and J. Bylander, Nat. Comm. 7, 11417
EP (2016).
[25] 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).
[26] 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).
[27] K. M. Sliwa, M. Hatridge, A. Narla, S. Shankar,
L. Frunzio, R. J. Schoelkopf, and M. H. Devoret,
Phys. Rev. X 5, 041020 (2015).
[28] F. Lecocq, L. Ranzani, G. A. Peterson, K. Cicak,
R. W. Simmonds, J. D. Teufel, and J. Aumentado,
Phys. Rev. Applied 7, 024028 (2017).
[29] F. Lecocq, J. B. Clark, R. W. Simmonds, J. Aumentado, and J. D. Teufel, Phys. Rev. Lett. 116,
043601 (2016).
[30] J. I. Cirac, A. S. Parkins, R. Blatt, and P. Zoller,
Phys. Rev. Lett. 70, 556 (1993).
[31] P. Rabl, A. Shnirman, and P. Zoller, Phys. Rev.
B 70, 205304 (2004).
[32] A. Kronwald, F. Marquardt, and A. A. Clerk,
Phys. Rev. A 88, 063833 (2013).
[33] E. E. Wollman, C. U. Lei, A. J. Weinstein, J. Suh,
A. Kronwald, F. Marquardt, A. A. Clerk, andK. C. Schwab, Science 349, 952 (2015).
[34] F. Lecocq, J. B. Clark, R. W. Simmonds, J. Aumentado, and J. D. Teufel, Phys. Rev. X 5, 041037(2015).
[35] J.-M. Pirkkalainen, E. Damskägg, M. Brandt,
F. Massel, and M. A. Sillanpää, Phys. Rev. Lett.
115, 243601 (2015).
[36] 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).
[37] J. K. Asbóth, P. Adam, M. Koniorczyk, andJ. Janszky, Eur. Phys. J. B 30, 403 (2004).
[38] K. Fang, J. Luo, A. Metelmann, M. H. Matheny,
F. Marquardt, A. A. Clerk, and O. Painter, Nature
Physics 13, 465 EP (2017).
[39] H. Goto, Phys. Rev. A 93, 050301 (2016).
[40] H. Goto, Sci. Rep. 6, 21686 (2016).
[41] F. Minganti, N. Bartolo, J. Lolli, W. Casteels, andC. Ciuti, Sci. Rep. 6, 26987 (2016).
[42] S. Puri, S. Boutin, and A. Blais, npj Quant. Inf.
3, 18 (2017).
[43] P. D. Drummond and D. F. Walls, J. Phys. A:
Math. Gen. 13, 725 (1980).
[44] K. Stannigel, P. Rabl, and P. Zoller, New J. Phys.
14, 063014 (2012).
10

[45] 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).
[46] 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).
[47] V. E. Manucharyan, J. Koch, L. I. Glazman, andM. H. Devoret, Science 326, 113 (2009).
[48] 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.
[49] 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).
[50] G. Zhu, D. G. Ferguson, V. E. Manucharyan, andJ. Koch, Phys. Rev. B 87, 024510 (2013).
[51] N. E. Frattini, U. Vool, S. Shankar, A. Narla, K. M.
Sliwa, and M. H. Devoret, Appl. Phys. Lett. 110,222603 (2017).
[52] W. Zhang, W. Huang, M. E. Gershenson, and
M. T. Bell, Phys. Rev. Applied 8, 051001 (2017).

Accepted in Quantum 2018-03-23, click title to verify

11