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

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.


Introduction
To harness the full potential of continuous variable quantum computing [1], it is necessary to prepare and ideally 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 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. [3][4][5][6]).
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 [7]. These states are integral to bosonic codes for quantum computing [8][9][10], and are equivalent to so-called "entangled coherent states", which have applications in many aspects of quantum information [11]. 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.
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 [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 (see e.g. [15][16][17][18][19][20]), where two driven bosonic modes with onsite Kerr interactions are coupled via coherent tunneling. We make some crucial modifications to this stan-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), 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 [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, something that is not possible for the standard driven-dissipative Bose Hubbard dimer [17]. 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 catstate encoded qubit introduced in Ref. [9]. 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. [23][24][25], strong onsite Kerr interactions [26] and engineered dissipation [27][28][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][31][32]. This protocol was recently realized experimentally both in optomechanics [33][34][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 † , whered 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: with a jump operatorẑ given bŷ Γ represents the rate of the dissipative interaction (i.e. the cooling rate), and the parameters µ, ν are set to µ = cosh r, ν = sinh r, making the jump operatorẑ equal to the Bogoliubovmode annihilation operatorβ [r]. The unique steady state of this master equation is then the dark state of the jump operatorẑ, which is nothing else than the vacuum squeezed state with squeeze parameter r.

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 z. One simple way to do this is to make the parameter µ a photon-number dependent operator. We take µ = µ 1d †d , and consider a master equation of the form in Eq. (1) with jump operator The steady states of this master equation correspond to the dark states of the jump operatorẑ. 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 |α as we find that our jump operator has dark stateŝ where the amplitude α is given by 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 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ẑ As we show explicitly in Appendix A, for µ 0 = 0, the only possible dark states of this dissipator are evenparity states. We thus find a unique dark state (and hence pure steady state of Eq. (1)): 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 µ 0 , this state asymptotically approaches an even-parity cat state: Even for small but finite µ 0 /µ 1 , |C[µ 0 , µ 1 , ν] 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.

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ĉ which couples nonlinearly to the principle moded of interest. To that end, we consider a two-mode system whose non-dissipative dynamics is described by the Hamiltonian The quadratic terms in this Hamiltonian correspond to a 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ĉ mode (loss at rate γ), such that the full evolution of the two modes is given by In the limit where theĉ mode is heavily damped (i.e. γ Λ), standard adiabatic elimination yields a master equation of the form in Eq. (1) with jump operatorẑ cat and with Γ ∝ Λ 2 /γ. In this limit, our two-mode system has a unique steady state analogous to that in Eq. (9): As discussed earlier, the state |C[µ 0 , µ 1 , ν] d asymptotically approaches an even cat-state in the limit µ 0 /µ 1 → 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Ĥ 1 andĉ 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 γ, Λ, µ 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 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ĉ andd, 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â andb: The dissipative dynamics described in Eq. (12) stabilizes the state |ψ ss defined in Eq. (13). Transforming this state to the "localized"â,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 [37]. We thus find that in the localized basis (and for µ 0 µ 1 ), our steady state has the form of a two-mode entangled cat state: (15) To determine the feasibility of our scheme, we need to express the HamiltonianĤ 1 in terms of the localized mode operatorsâ,b. Before doing this, we first modifŷ H 1 to make it more symmetric, but without changing the steady state: We again consider the master equation It is easy to see that the state |ψ 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 Hamil-tonianĤ 2 in the localized mode basis. Defining for convenience we obtain 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 ±∆, 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 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 [21,22]. 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 [39][40][41][42]; our scheme leverages similar resources in two cavities (along with engineered dissipation) to now stabilize an entangled cat state. Note that the amplitude α of our cat state (c.f. Eq. (7)) in terms of the energy scales in the Hamiltonian is now α = i λ/U .

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ĉ,d mode basis. In this basis the steady state |ψ ss has theĉ mode in vacuum, while thed 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 thed mode is given by: Here |2n d 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 existence of a pure steady state; it is much simpler than more general 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 in cascaded quantum systems; in contrast to that work, our approach does not require cascaded interactions, or even any breaking of time-reversal symmetry. The analytic expression for |ψ ss 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 existence of a steady state, and |ψ ss tends to the desired even cat-state. Fig. 2 uses the exact expression of the steady state |C d to show this evolution as a function of the ratio ∆/U = µ 0 /µ 1 . The upper panel depicts the Wigner function for thed mode at three specific values of ∆/U . As described in the figure caption, this describes a cut through the two-mode Wigner function describing the state in the localizedâ,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 even for modest values of ∆/U , the fidelity with the desired even-parity cat state can be extremely high.

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  15)). The fidelity is defined by assess the impact of such losses, we add single-photon loss (at a rate κ in each cavity) to our master equation: 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|α| 2 . As a result, 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 (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 κ. We find that for ∆/U = 1, we get F ≈ 0.95 with the exact solution and F ≈ 0.86 with the ideal cat state for a loss rate κ = 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 ∼ 100-1000 kHz has been realized [26], while in stateof-the-art 3D cavities internal leakage rates are on the order of 10-100 Hz [45,46].
Note that for the results in Fig. (3), we have taken the coupling γ to the engineered dissipation to be γ = √ 2∆. 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 γ. 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.

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 ∆ = 0, one would instead stabilize a two dimensional manifold spanned by the entangled cat states (the states |2 and |1 defined in Eq. (27)).

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 steady state
We derive here the exact steady state of the system described by the HamiltonianĤ 2 in Eq. (16) and the master equation in Eq. (17). Motivated by the fact that only theĉ mode experiences any damping, we make an ansatz that there exists a pure-state steady state in which theĉ mode is in vacuum, i.e.
Here, α n are arbitrary normalized coefficients and |n d are Fock states. The above state is trivially a dark-state of the jump operatorĉ 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 |ψ withĤ 2 yields, where we have expressed µ 0 , µ 1 and ν in terms of ∆, U and λ using Eq. (18). The RHS of Eq. (25) has theĉ mode in the state |1 , while in |ψ it is in the state |0 . Thus, the only way that |ψ can be an eigenstate ofĤ 2 is if it has zero energy. We thus needĤ 2 to annihilate |ψ , and therefore require that the coefficient of each |n d state in Eq. (25) vanish. This leads to the recursion relation The above recursion relation does not mix the even and odd-parity subspaces. Note that for ∆ = 0, the n = 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 the solution 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 steady states. Thus, as discussed in the main text, a non-zero detuning ∆ is crucial in order to obtain a unique steady state.
While the above discussion has identified a pure steady state for our system when ∆ = 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 λ ≤ ∆, 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 ∆ is non-zero, the steady state described by Eq. (21) is indeed the unique steady state.

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 ∆ = 0, each cavity has exact zero-energy eigenstates |C ± (α) [39,40,42]. For the interesting limit of a large cat state amplitude α, 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: The states |2 , |4 have even total photon number parity, while the states |1 , |3 have an odd total photon number parity. It is useful to make a slight rotation of the even parity states: where for large |α| we have Note that in the large |α| limit, |2 → |2 , which is our expected entangled-cat steady state. The states |1 , |2 , |3 , |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 |1 and |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Ĥ 2 becomes a coherent tunneling term of am-plitude∆ between the odd parity states |1 and |3 , whereas the dissipator gives rise to the three depicted incoherent transitions (with rates Γ 3,4 , Γ 4,3 , Γ 2,3 ). One sees clearly from this schematic that, as expected, the dark state |2 will be our eventual steady state.
For large |α|, one finds that∆ and Γ 2,3 become exponentially small:∆ The exponential suppression of these rates corresponds to the extremely weak overlap of the coherent states |α and |−α . The above processes will then form the bottleneck in having populations relax to the state |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 α limit, we find that the rate-limiting slow relaxation between the dark˜ states |1 and |2 is described by where we have defined the effective incoherent tunneling rate between the state |1 and the {|3 , |4 } bright-state manifold as: The rate Γ rel describes an effective two step process where population first incoherently tunnels from |1 to the bright state manifold (rate Γ ∆ ), and then relaxes to the steady state |2 (rate Γ 2,3 ). Note that −1 , as the incoherent tunneling process described by Γ ∆ can also move population the wrong way, i.e. from |3 back to the state |1 .
It is instructive to consider the behaviour of Γ rel for small and large ∆/γ: As expected, the relaxation rate vanishes in both these extreme limits, as the relaxation process from |1 to |2 involves both the∆ ∝ ∆ tunneling process and 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 used in Fig. 3.

C Implementation issues C.1 Parametric drives
In the lab frame, the Hamiltonian for two parametrically driven cavities iŝ where ω a/b are the cavity resonance frequencies and ω 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 ω pa /2 for cavityâ and ω p b /2 for cavityb, and choose the parametric drive frequencies such that It is most convenient to consider a situation where the parametric drives are at the same frequency, i.e. ω pa = ω p b ≡ ω 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ω 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 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 ω 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â +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 timeindependent. In this case, a correlated dissipator can still be achieved using frequency conversion via engineered dissipation with active elements (see e.g. [27][28][29]).

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. [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 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λ (where m is an integer, andλ is the 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/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. [26]), 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 [47]). The effective Hamiltonian describing such a system has the form [48] 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.φ (n) is the junction phase (charge) operator. The parameter φ g describes a flux 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 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 [49], 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 [50] yields an effective self-Kerr interaction for the cavity given bŷ . (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 U qb is the anharmonicity of the transmon. For a transmon U qb ≈ −E C , i.e. U qb < 0. If δ < 0, or δ > 0 and δ + U qb /2 > 0, then U eff < 0 and the cavity nonlinearity is attractive. However, if δ > 0 and δ + U qb /2 < 0, which lies in the straddling regime (defined by δ + U qb < 0), then U 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 [51,52].

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 moded correspond to the entanglement between the local modesâ andb. As can be seen, the fringes are still clearly visible in the lossy case, indicating that entanglement between the local modes remains.

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   23)). We use system parameters of λ/U = 2, ∆/U = 0.5, γ/U = 1/ √ 2, which for (b) corresponds to a fidelity of F ≈ 0.9, as seen in Fig. 3 where δ∆, δλ, and δU are the dimensionless mismatch in the cavity detuning, parametric drive strength, and self-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 |0 c ⊗ 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.   F (c) 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 |0 c ⊗ C d of Eq. (21). The mismatch parameters δ∆, δλ, δU are dimensionless relative differences between the two cavities. For all curves the parametric drive strength λ/U = 2, the detuning ∆/U = 1, and the correlated dissipation rate γ/U = 2. As can be seen, F > 0.95 can be obtained for up to 5% relative mismatch.