Stabilization of multi-mode Schrodinger cat states via normal-mode dissipation engineering

Non-Gaussian quantum states have been deterministically prepared and autonomously stabilized in single- and two-mode circuit quantum electrodynamics architectures via engineered dissipation. However, it is currently unknown how to scale up this technique to multi-mode non-Gaussian systems. Here, we upgrade dissipation engineering to collective (normal) modes of nonlinear resonator arrays and show how to stabilize multi-mode Schrodinger cat states. These states are multi-photon and multi-mode quantum superpositions of coherent states in a single normal mode delocalized over an arbitrary number of cavities. We consider tailored dissipative coupling between resonators that are parametrically driven and feature an on-site nonlinearity, which is either a Kerr-type nonlinearity or an engineered two-photon loss. For both types of nonlinearity, we find the same exact closed-form solutions for the two-dimensional steady-state manifold spanned by superpositions of multi-mode Schrodinger cat states. We further show that, in the Zeno limit of strong dissipative coupling, the even parity multi-mode cat state can be deterministically prepared from the vacuum. Remarkably, engineered two-photon loss gives rise to a fast relaxation towards the steady state, protecting the state preparation against decoherence due to intrinsic single-photon losses and imperfections in tailored dissipative coupling, which sets in at longer times. The relaxation time is independent of system size making the state preparation scalable. Multi-mode cat states are naturally endowed with a noise bias that increases exponentially with system size and can thus be exploited for enhanced robust encoding of quantum information.


I. INTRODUCTION
Schrödinger cat states-quantum superpositions of macroscopically distinct (or "classical") states-are a fundamental resource for quantum communication [1], quantum metrology [2][3][4] and quantum computation [5][6][7]. Recently, the development of bosonic quantum error-correcting codes has led to a surge of interest in exploiting Schrödinger cat states for encoding and manipulating quantum information in an hardware-efficient manner [6][7][8][9][10][11][12][13][14][15][16][17]. The two-dimensional subspace spanned by superpositions of Schrödinger cat states can be used to encode quantum information, implementing a socalled cat qubit. Interestingly, such a qubit is naturally endowed with a large noise bias, i.e. the qubit has a single dominant noise channel while all other types of noise are largely suppressed, which offers protection from errors in a quantum memory [8,11], and eases up the requirements of quantum error correction [13] and quantum annealing [18]. The implementation of bias-preserving quantum gates [7,[14][15][16] opens the door for universal fault-tolerant quantum computing based on cat qubits [6]. The potential for applications of Schrödinger cat states further increases when considering extensions to multiple modes. Two-mode Schrödinger cat states, also known as entangled coherent states [19], are a resource for continuous-variable quantum information [20,21] and metrology [22,23]. They have been experimentally prepared using feedback control in superconducting circuits [24]. A unique feature of these states is that they possess non-Gaussian entanglement, which allows circumventing many no-go theorems imposed by Gaussian quantum resources [25]. Extending the stabilization of non-Gaussian entangled states-and Schrödinger cat states in particular-to an arbitrary number of modes would render non-Gaussian quantum resources scalable [26].
The widespread application of Schrödinger cat states to photonic quantum technologies is also due to the unique control capabilities available in circuit quantum electrodynamics (circuit QED) architectures, which guarantee unprecedented control of nonlinear interactions and dissipation [27][28][29]. In this context, a particularly convenient approach is dissipation engineering [30,31], which exploits tailored interactions with a dissipative environment to attain deterministic and robust preparation of desired target states or operations. The dissipative preparation of single-mode Schrödinger cat states using engineered two-photon loss and two-photon (parametric) drive have been first envisioned in Ref. [32]. For a circuit QED setup, the autonomous stabilization of a twodimensional subspace spanned by single-mode cat states has been recently proposed in Ref. [9] and realized in Ref. [10]. It has also been proposed how to generate approximate twomode Schrödinger cat states via dissipation engineering [33]. However, engineering multi-mode Schrödinger cat states and the stabilization of steady-state manifolds thereof in drivendissipative parametrically coupled systems have remained an unexplored avenue.
In this work we propose how to dissipatively generate multi-mode cat states in a scalable fashion and show how to autonomously stabilize a decoherence-free subspace (DFS), which is spanned by multi-mode Schrödinger cat states. Our approach consists in upgrading dissipation engineering to collective modes (normal modes) of dissipatively coupled arrays of nonlinear resonators. We present our ideas in terms of two closely related models that are readily realizable with parametrically driven superconducting circuits. The first model is, as sketched in Fig. 1a, based on an array of parametrically driven Kerr resonators, also known as Kerr parametric resonators (KPRs). In the second model, the Kerr nonlinearity is replaced by engineered two-photon dissipation, depicted in Fig. 1b. In both cases, the only source of coupling is provided by engineered nonlocal dissipation connecting neighboring pairs of modes. This kind of non-local dissipation has recently been considered in the proposal for the generation of two-mode Schrödinger cat states [33] and in several other contexts, e.g. to realize non-reciprocal photon transport in a pair of modes [31] or in cavity arrays [34,35]. In contrast, the more conventional scenario of arrays of coherently (tunnel-) coupled KPRs have been studied in Refs. [36,37], in connection with the emergence of critical behavior. Fully connected networks of tunnel-coupled KPRs have also been proposed for solving combinatorial optimization problems [38,39].
As we shall see, the choice of the pairwise dissipative coupling is responsible for the unique features of our model. The idea is simple and yet extremely powerful, as it allows to select a single normal mode and completely suppress its decay, while all the other normal modes can be heavily damped. The interaction between the (driven) non-dissipative normal mode and the other dissipative modes leads to an asymptotically stable Schrödinger cat manifold for the non-dissipative mode, which we characterize analytically. Since these cat states are stabilized in a normal mode, they are indeed multi-mode superpositions delocalized over the entire array. We further show that the stabilization of multi-mode cat states is robust against intrinsic single-photon loss and imperfections in the dissipative coupling. In particular, the model featuring local two-photon loss [cf. Fig. 1b] guarantees a fast relaxation towards the asymptotic manifold, which allows for an efficient steady-state preparation of the even parity cat state from the vacuum even in the presence of intrinsic single-photon loss. Multi-mode cat states can be stabilized and prepared in the Zeno limit even in the presence of imperfections in non-local dissipation, albeit with decreased fidelity as an additional decoherence channel sets in. The relaxation time towards the Schrödinger cat manifold is independent of system size, making the state preparation scalable. Crucially, the noise bias of multi-mode cat states exponentially increases with system size, which can be exploited for realizing a quantum memory with enhanced protection with respect to single-mode implementations.
The remainder of this work is structured as follows. In Sec. II we introduce our first model, namely the Kerrresonator array of Fig. 1a. In Sec. III we derive an analytic expression for the steady states of the model, which take the form of a two-dimensional DFS spanned by even and odd multi-mode Schrödinger cat states. In Sec. IV we discuss how to target such states, and we show that, for sufficiently strong non-local dissipation, the even multi-mode cat state can be faithfully prepared by initializing the array in the vacuum. In Sec. V we explain this behavior by developing an effective single-mode theory, valid in the Zeno limit of strong nonlocal dissipation. We describe the long-time dynamics of the nondissipative mode and extract an analytical expression of the dissipative gap, dictating the slowest convergence time to the FIG. 1. Two cavity arrays considered for the stabilization of multimode Schrödinger cat states; straight (wiggly) lines represent Hamiltonian (dissipative) contributions. (a) One-dimensional array of N resonatorsâj with on-site Kerr nonlinearity U , non-local dissipation at rate γ, parametric two-photon drive with amplitudes Gj = Ge −iθj , intrinsic photon loss at rate κ, and periodic boundary con-ditionsâN+1 =â1. The inset shows that non-local dissipation is implemented via coupling to an engineered reservoir. Note that each dissipative process (wiggly lines) is due to coupling to an independent reservoir. (b) Alternative model where the Kerr nonlinearity is replaced by local two-photon loss at rate η.
steady state. The expression shows that the requirements of preparing a cat state from the vacuum with near-unit fidelity and guaranteeing a fast convergence cannot be simultaneously satisfied. In Sec. VI we then show how this drawback can be remedied by replacing the Kerr non-linearity with local twophoton loss as shown in Fig. 1b. In Sec. VII we account for the effect of unwanted decoherence mechanisms. We include in both models intrinsic single-photon loss, which causes decoherence within the steady-state manifold. In Sec. VIII, we discuss effects of decoherence due to imperfections in nonlocal dissipation. In Sec. IX, we discuss the noise bias of the multi-mode cat manifold and how it can be exploited for a robust encoding of quantum information. In Sec. X, we estimate experimentally relevant parameters of our model based on state-of-the-art circuit QED technology. Sec. XI contains concluding remarks and possible directions for future investigation.

II. MODEL
We consider a one-dimensional array of N identical Kerr resonators, as depicted in Fig. 1a, each subject to a parametric two-photon drive. In a frame rotating at the resonator frequency and for the pump frequency matching twice the resonator frequency, the system is described by the Hamiltonian where U is the strength of the Kerr nonlinearity and G j = G e −iθj is the amplitude of the two-photon drive at site j with a phase difference θ between neighboring resonators; without loss of generality we assume G > 0.
The distinctive feature of our model lies in the coupling between resonators, which is not described by a Hamiltonian term, but has instead a dissipative nature, inasmuch as it correlates photon emission events from neighboring resonators. This is formally described by a non-local dissipation term in the Lindblad master equatioṅ where γ is the rate of non-local dissipation, L is the Liouvillian superoperator and D[Ô]ρ =ÔρÔ † − 1 2 Ô †Ôρ +ρÔ †Ô . Loosely speaking, this unconventional dissipator can be interpreted as a dissipative analogue of a hopping term. In fact, it corresponds to indirect hopping via an auxiliary, fast decaying, mode, which acts as an engineered reservoir, see the inset of Fig. 1a. It can be implemented in different ways, e.g. via coupling the neighboring resonatorsâ j andâ j+1 to a transmission line with respective coupling amplitudes having a relative phase difference φ [31,33].
We assume periodic boundary conditionsâ N +1 =â 1 . In this case, non-local dissipation can be expressed in the plane- where the rate of dissipation γ k = 2γ [1 − cos (k − φ)] depends on the quasi-momentum k which takes the values k = 2π N , 2 2π N , ..., 2π. From Eq. (4) we see that the effect of engineered nonlocal dissipation is best understood in momentum space, where it takes the form of a non-uniform damping for the normal modes. In particular, the normal mode with k = φ is a dark mode of the non-local dissipator, which does not experience any dissipation as γ φ = 0. Therefore the relative phase φ allows to select a single non-dissipative normal mode. This behavior is to be contrasted with intrinsic, i.e. non-engineered, single-photon loss, which instead leads to homogeneous decay of all normal modes. For now we assume that (4) is the only dissipation channel acting on the system; effects of intrinsic photon loss will be accounted for in Sec. VII.
Writing the Hamiltonian in the plane-wave basis we get where, δ k1+k2,k3+k4 is Kronecker delta and the arguments k 1 + k 2 and k 3 + k 4 are defined modulo 2π. We can see that the Kerr nonlinearityĤ U corresponds to a four-wave mixing process, which scatters a photon pair k 1 and k 2 to a photon pair k 3 and k 4 conserving the total quasi-momentum. The two-photon driveĤ G corresponds to the creation of a photon pair with a total quasi-momentum θ.

III. MULTI-MODE DECOHERENCE FREE SUBSPACE
Our first goal is to characterize the steady states of the master equation (3). In particular, we look for steady states in the form of pure states |Ψ , which are referred to as dark states [40]. Combining Eqs. (3) and (4) it readily follows that dark states must satisfy the following two conditionŝ The condition (7) is satisfied if all modes k = φ are in the vacuum, i.e. |Ψ = k =φ |0 k ⊗ |ψ φ , from which it is clear that the role of engineered dissipation is to select a single normal mode that is (possibly) populated in the infinite-time limit. The dark state condition (8) can then be written as where ζ = i N G U and we hereafter set θ = 2φ [41]. For = Gζ 2 , overall destructive interference is guaranteed if the degenerate photon pairs annihilated in mode φ satisfy the following condition (10) in which case both sides of Eq. (9) vanish.
Eq. (10) represents a simplified dark state condition. From its factorized form we can deduce that the two coherent statesb φ | ± ζ φ = ±ζ| ± ζ φ and their superpositions are dark states. All possible dark states span a two-dimensional subspace shown in Fig. 2a, whose basis states are superpositions of coherent states, known as Schrödinger cat states |C ± φ = N ± (|ζ φ ± | − ζ φ ), where N ± = 2 1 ± e −2|ζ| 2 −1/2 . Note that the cat states |C ± φ are exactly orthogonal in contrast to the coherent states |±ζ φ , which have a finite overlap. For convenience, from now on we drop the subscript φ from the state of the non-dissipative where the coefficients c µ depend on the initial state. Even though steady states which are neither dark states nor their incoherent mixtures can in principle exist [40], we numerically verified that the master equation (3) has no steady states lying outside of the DFS (11) (see Appendix A for more details). Notice that DFSs with a similar structure have recently been realized in circuit QED [7,10,14,17,42]. However, compared to these cases, here the crucial difference is that the DFS is spanned by multi-mode Schrödinger cat states. In fact, the ground state of a single non-dissipative KPR supports a twodimensional manifold spanned by a single-mode cat states of amplitude ζ = i G/U [12,43]. The steady state of a parametrically driven resonator with engineered two-photon dissipation (hereafter simply referred to as two-photon driven dissipation) [9] also supports a similar manifold, in this case with amplitude ζ = −iG/η, where η is the rate of two-photon dissipation (see Sec. VI for more details). We now show that in our case the states in the DFS (11) are multi-mode.
The cat states are encoded in an arbitrary normal mode k = φ. As the normal modes are delocalized over the entire array, the cat states that we discuss here are highly nonlocal. This can be seen explicitly by moving into the basis of local modes where ζ j = ζ √ N e −ijφ [44]. Controlling the phase φ (and that of the drive θ = 2φ) it is possible to prepare different instances of these cat states. In particular, if φ is an even multiple of π, Eq. (12) resembles a GHZ-like state states, with real and opposite amplitudes; if φ is an odd multiple of π, the amplitudes being superimposed are staggered.
These pure states are the multi-mode version of standard Schrödinger cat states [45]. They have the remarkable property of being both multi-photon and multi-mode quantum superpositions, and possess genuine multipartite entanglement [21,26,46]. Moreover, such entanglement is of the most useful kind, namely non-Gaussian, which makes (12) a resource for quantum metrology [2,22,23], quantum information based on continuous variables [20,21,25] and continuous-variable quantum computation [47]. We also stress the potential of multi-mode cat states is still largely unexplored. In fact, theoretical studies of multi-mode cat states has been mostly limited to N = 2 (also known as entangled coherent states [21]) due to the current lack of preparation protocols that are at the same time scalable and robust to errors.
To contrast the multi-mode cat states with standard singlemode cat states, we plot in Fig. 2b, 2c and 2d the N -mode Wigner function W (x 1 , ..., x N , p 1 , ..., p N ) of the cat states |C ± stabilized in the normal mode φ = 0 [48]. The unique feature of multi-mode cat states is the quantum interference in multi-mode phase space, which manifests itself in interference fringes of the Wigner function in the x j -x l plane (for any j = l) depicted in 2c and 2d for |C + and |C − , respectively. This feature is absent for single-mode cat states.
An interesting observation comes from comparing our protocol with known techniques for the generation of multi-mode cat states in linear optics [2]. In a linear optics network, multimode cat states can be prepared by mixing a single-mode cat state |C ± with N − 1 vacua. The optics network is represented by a unitaryT performing the linear transformation a j =Tb 2πj/NT † = k T jkbk of the input modesb k , where T jk = 1 √ N e −ijk . In our case, as depicted in Fig. 2e, this linear transformation is automatically implemented by engineering the state of a normal mode. The obvious differences are that our method allows to stabilize a cat state in stationary (cavity) modes rather than itinerant modes, and that it does not rely on an external resource (i.e., the input single-mode cat state). Besides that, our strategy possesses yet another fundamental advantage. For a single-mode cat state of amplitude ζ at the input of an interferometer, the corresponding multi-mode cat state at the output has rescaled local amplitudes ζ/ √ N ; this feature reflects the passive character of the optical network, which does not pump energy into the system. On the other hand, in the driven-dissipative process that we exploit to stabilize the multi-mode cat state, the total power ∝ N G injected into the array increases with the number of resonators N . This in turn allows for the amplitude of the state |ζ| ∝ √ N to increase with N . As a result, the multi-mode cat states can be distributed over an arbitrary number of local modes N with the local amplitudes |ζ j | = |ζ|/ √ N = G/U independent of N . This scaling of amplitudes allows for increasing noise bias by using multi-mode cat states instead of single-mode cat states, which we discuss in detail in Sec. IX.

IV. STEADY-STATE PREPARATION OF MULTI-MODE CAT STATES
In the last section we showed that the steady states of Eq. (3) belong to a DFS spanned by multi-mode cat states. Therefore, depending on the initial stateρ in , the state of the array asymptotically converges towards a pure dark state α|C + +β|C − or their incoherent mixtures. Given the unique properties displayed by these cat states [cf. the discussion following Eq. (12)], we now focus on the stabilization of |C ± , i.e., we wish to find initial statesρ in for which the array asymptotically approaches these states. At the same time, we should keep in mind that the choice of the initial stateρ in must be experimentally feasible.
In analogy with the preparation of a single-mode even parity cat state, which can be stabilized starting from the vacuum [9,12], we initialize the array in the N -mode vacuum stateρ in = |0 0|. In Fig. 3 we plot the fidelity c ++ (c −− ) of the steady state with the multi-mode cat state |C + (|C − ) as a function of the non-local decay rate for different two-photon pump strengths (see Appendix A for more details about the numerical simulations). We can see that c ++ (blue lines in Fig. 3a) is in general smaller than unity and c −− (green lines in Fig. 3b) is non-zero. Therefore, at variance with the singlemode case, starting from the vacuum does not guarantee to reach the target steady state |C + . Nonetheless, Fig. 3 provides a clear indication that, for sufficiently strong non-local dissipation, fidelities with |C + arbitrarily close to unity can be achieved, which renders the preparation of multi-mode cat states from vacuum feasible in the strong dissipation regime.
The cat state |C + can be transformed into a superposition cos α|C + + i sin α|C − of cat states |C + and |C − by a rotation within the stabilized DFS induced by a weak singlephoton drive as demonstrated for single-mode cat states [14] (see Appendix C for more details).
The target steady state |C + is not deterministically approached from the initial vacuum for all γ since photon parity is not conserved during the time evolution [49]. The multimode cat states |C ± spanning the DFS have a well-defined parity as they are ±1 eigenstates of the generalized parity op-eratorP = exp iπ N j=1â † jâ j . However, the photon parity is not a conserved quantity of the master equation (3) -a conserved quantity is an operatorĴ such that L †Ĵ = 0, while For the initial vacuum, the system starts off in the evenparity subspace but during the time evolution it leaks to the odd-parity subspace leading to non-unit (non-zero) values of the fidelity c ++ (c −− ) shown in Fig. 3. For γ ∼ G, U this leakage leads to a decreased purity p of the steady state (orange lines in Fig. 3a), while for sufficiently strong dissipation the leakage is suppressed and the generalized parity is approximatively conserved (to an arbitrary degree of accuracy). This further confirms that the regime of strong nonlocal dissipation is relevant for steady-state preparation of the even cat state |C + . In this regime, the dynamics at all times can be captured by an effective single-mode theory as we will show in the next section. On the other hand, due to nonconservation of parity and the many-body interactionĤ U , for values of the dissipation rate comparable to the Kerr nonlinearity we expect the transient dynamics to exhibit genuine many-body features.

V. THE ZENO LIMIT
We now have numerical evidence that, in the limit of dominant non-local dissipation, the even multi-mode cat state |C + can be deterministically prepared from the N -mode vacuum. In order to better explain this behavior, in this section we develop an analytical insight into the regime of large nonlocal decay rate γ. For γ G, U , the many-body interac-tionĤ U between normal modes is suppressed: all decaying modes k = φ are strongly damped and only the single normal mode φ is populated. Transitions from the vacuum state to excited states of the decaying modes k = φ are inhibited, i.e. they are in the vacuum at all times. This is analogous to the quantum Zeno effect [51, 52], i.e. the inhibition of quantum transitions due to a frequent projective measurement of a quantum system. In this analogy, the measurement is enacted by the engineered environment via the non-local dissipation , which continuously projects all decaying modes onto the vacuum state. In the following, we refer to the limit of a large non-local decay rate γ G, U as the Zeno limit.
In Appendix B we show that in this limit the reduced density matrixρ φ of the non-decaying mode φ is described by the effective master equatioṅ where L φ is the effective single-mode Liouvillian. This effective master equation is derived by treating the Hamilto-nianĤ in (3) as a perturbation to the dominant dissipation L γ and projecting all decaying modes k = φ onto the vacuum state. Perturbation theory to first order inĤ describes a unitary Zeno dynamics governed by the effective Hamiltonian namely that of a singlemode KPR [12]. On the other hand, processes which are of second order inĤ lead to an effective two-photon driven dissipation term described by the dissipator D b2 φ − ζ 2 and occurring at a rate Γ = 4 U 2 N 2 k =φ 1 γ k ; this effective dissipation is facilitated by virtual excitations of the decaying modes k = φ. We thus come to an important conclusion: the emergence of effective two-photon dissipation is the mechanism responsible for the relaxation towards the steady state of the non-decaying mode.
We note that, since it involves second-order processes, dissipation sets in only at long times t ∼ 1 Γ ∼ γN 2 U 2 . At short times t ∼ N U , the mode φ undergoes a unitary Zeno dynamics governed by the effective HamiltonianĤ φ . A general (model-independent) treatment of Lindblad master equations with strong dissipation acting on a subset of the degrees of freedom, is discussed in detail in Ref. [53], including the derivation of the unitary Zeno dynamics and effective weak dissipation of non-decaying modes.
In Fig. 3b, we plot the difference ρ ss −ρ d ⊗ρ φ,ss (purple lines) between the steady stateρ ss of the full master equation (3) and the steady stateρ d ⊗ρ φ,ss of the effective master equation (13) as a function of the decay rate γ for different values of G U , whereρ d = k =φ |0 k 0| is the (N − 1)-mode vacuum state and Â 2 = Tr Â †Â is the Hilbert-Schmidt norm. We can see that the difference vanishes as O U 2 γ 2 with the increasing decay rate γ. This confirms that, for large enough γ, the effective master equation (13) accurately describes the dynamics of the Kerr-resonator array at long times.
In the Zeno limit, the photon parity P φ = (−1)b † φb φ is a conserved quantity of the effective Liouvillian, i.e., L † φ P φ = 0, as it commutes with both the HamiltonianĤ φ and the jump operatorb 2 φ − ζ 2 [50]. Therefore, in this limit there is no leakage to the odd-parity subspace. This provides a neat explanation for the behavior observed in Fig. 3a for large γ and proves that the dissipative dynamics of the master equation (3) allows for the steady-state preparation of the even cat state |C + from the vacuum state |0 with a fidelity arbitrarily close to unity.
The rate of convergence towards the steady state is given by the dissipative gap ∆ d , which is determined from the spectrum of the Liouvillian as the smallest non-vanishing real part of the eigenvalues [50]. Since the rate of two-photon dissipation is inversely proportional to the non-local decay rate, i.e., Γ ∝ γ −1 , the dissipative gap ∆ d ∝ γ −1 retains the same dependence (see Appendix F for more details). For |ζ| > ∼ 1, the dissipative gap can be approximated as Eq. (14) is an important result as it expresses a trade-off between the suppression of the leakage to the odd-parity subspace and the time required to reach the steady-state manifold. The price of an accurate (i.e., near-unit-fidelity) preparation of the even multi-mode cat state |C + (from the vacuum) is a slow convergence. Although this does not per se preclude the stabilization of |C + , it becomes problematic in the presence of competing decoherence mechanisms, such as intrinsic photon loss. In Sec. VII, we will discuss the impact of this source of decoherence in the Kerr-resonator array and how it limits the steady-state preparation. Moreover, increasing the number of modes N in the superposition leads to further increase of the convergence time, which poses a theoretical limit to the scalability of our approach.

VI. ALTERNATIVE MODEL INCLUDING TWO-PHOTON LOSS
To avoid the trade-off between the suppression of the leakage to the odd-parity subspace and the rate of convergence towards the steady state, we consider an alternative model which allows the stabilization of multi-mode cat states |C ± while at the same time featuring a large dissipative gap. This model is based on an array of dissipatively coupled linear resonators (U = 0), subject to local two-photon loss at rate η and twophoton pump with the amplitude G, as sketched in Fig. 1b. The model only requires a small adjustment with respect to our previous model, namely replacing the on-site Kerr nonlinearity with engineered two-photon loss. This alternative is known to stabilize a DFS whose basis states are even and odd singlemode cat states [9,32]. The local dissipator has been realized experimentally in Refs. [10,14,17].
The dissipatively coupled array of resonators described by the master equation (15) exhibits the same DFS as the Kerrresonator array [see Eq. (11)] with a modified amplitude of the cat states ζ →ζ = e −i 3π 4 ζ (see Appendix D for a detailed derivation of the DFS). In the Zeno limit of strong non-local dissipation γ G, η, the time evolution of mode φ is described by an effective master equatioṅ which is derived by employing the first-order perturbation the- leading to the convergence towards the DFS. This is in contrast to the Kerr-resonator array, where the convergence towards the DFS is due to secondorder processes involving virtual excitations of modes k = φ.
For this second model we can completely suppress the leakage to the odd-parity subspace in the Zeno limit (the photon parity P φ is a conserved quantity) and at the same time maintain a large dissipative gap which increases with the strength G of the two-photon drive (see Appendix F for more details). Note that, in contrast to the Kerr-resonator array (14), the dissipative gap ∆ d is independent of γ in the Zeno limit. Importantly, the gap is also independent of system size N , making the preparation of cat states in large arrays feasible.

VII. EFFECTS OF INTRINSIC SINGLE-PHOTON LOSS
In all current circuit QED implementations, unwanted intrinsic photon loss is the dominant source of decoherence, both for the case of the Kerr-resonator array [7,54] and for the model with two-photon driven dissipation [10,14,17]. Photon loss at each resonator is described by the Liouvil- We take intrinsic loss into account by including it in the master equation (3), thus obtaininġ ρ = (L + L κ )ρ. In contrast to the non-local dissipation L γ , which leads to non-uniform dissipation in momentum space, intrinsic photon loss induces uniform dissipation of all normal modes, i.e., L κρ = κ k D[b k ]ρ. For modes k = φ, intrinsic photon loss only increases the rate of single-photon dissipation γ k → γ k + κ which results in their faster decay.
Crucially, intrinsic photon loss gives rise to the single- We remind the reader that the latter is due to local two-photon loss in the model of Sec. VI and arises due to the coupling to modes k = φ in the Kerr-resonator array, see Sec. V. The single-photon dissipation opens a decoherence channelb φ |C ± = ζ|C ∓ within the DFS (11) at rate |ζ| 2 κ [9]. This results in quantum jumps between even and odd cat states, leading to decoherence but not to leakage out of the DFS; even in the presence of singlephoton loss, the steady state is confined within the cat manifold. In addition, intrinsic photon loss also causes a decrease in the amplitude |ζ| of the cat states. However, this decrease is negligible provided that the intrinsic loss rate κ G is small (which is the case for state-of-the-art superconducting cavities [54]) as the two-photon drive quickly repumps depleted photons.
Due to the decoherence, the array (for both models) converges at times t ∼ 1/|ζ| 2 κ towards the fully mixed state ρ ss ≈ 1 2 (|C + C + | + |C − C − |). Therefore, for the deterministic preparation of |C + , it is important how the rate of decoherence |ζ| 2 κ compares to the rate of convergence towards this state, which corresponds to the dissipative gap ∆ d .
We start by investigating the effects of intrinsic photon loss on the Kerr-resonator array in the Zeno limit, i.e., we consideṙ We plot the fidelity c ++ of the instantaneous stateρ d ⊗ρ φ (t) with the even cat state |C + (blue lines) and the fidelity c −− with the odd cat state |C − (green lines) in Fig. 4a as a function of time for different values of the intrinsic loss rate κ and starting from the vacuum state. We can see that, small intrinsic loss κ/U = 10 −4 is required to reach a large fidelity c ++ ≈ 0.91 (blue solid line) before both fidelities c ++ and c −− converge towards the stationary values c ++ = c −− ≈ 0.5 and the fully mixed state is approached. On the other hand, for a larger loss rate κ/U = 10 −2 , the fidelity c ++ (blue dot-dashed line) does not exceed the value c ++ = 0.5 at any time as the rate of decoherence |ζ| 2 κ is larger than the rate ∆ d of convergence towards the even cat state. Since the dissipative gap ∆ d (14) is small in the Zeno limit γ G, U , the steady-state preparation in the Kerr-resonator array is largely limited by decoherence.
We now consider the impact of photon loss on our second model in the Zeno limit described by the master equatioṅ Fig. 4b we can see that a large fidelity c ++ builds up at times orders of magnitude shorter than the characteristic decoherence time 1 κ (|ζ| 2 ∼ 1 here) as the array quickly converges towards the even cat state. The quick convergence is guaranteed by a large dissipative gap ∆ d compared to the decoherence rate |ζ| 2 κ with the ratio For all values of the intrinsic loss rate κ, the maximal fidelity is orders of magnitude closer to unity than that for the Kerr-resonator array and is reached at orders of magnitude shorter times.
We conclude that local two-photon loss leads to a fast convergence towards the DFS, which allows for an efficient steady-state preparation in the presence of intrinsic photon loss. This is in contrast to the steady-state preparation in the Kerr-resonator array, which is severely limited by decoherence due to the slow convergence towards the DFS.

VIII. ROBUSTNESS AGAINST IMPERFECTIONS IN NON-LOCAL DISSIPATION
Leveraging strong and tunable non-local dissipation plays a crucial role in our preparation protocol, allowing to select a single non-dissipative normal mode and at the same time strongly damping the remaining modes. We now study the robustness of our protocol against imperfections in the nonlocal dissipation, which can appear due to the imperfect tuning of the couplings to the engineered reservoirs.
As the first example of imperfection, we consider a finite accuracy in tuning the phase of the non-local dissipators. Imagine we select a plane-wave mode φ ≡ k (whose dissipation we aim to suppress) but due to finite precision, the externally imposed phase φ = φ + δφ of the non-local jump operators differs from the target quasi-momentum φ by a small offset δφ, where we assume δφ > 0 for simplicity. Due to the quasi-momentum mismatch φ = φ , mode φ is no longer a dark mode of the non-local dissipator L γ and, as a consequence, it exhibits non-vanishing single-photon loss at the rate γ φ = 2γ[1 − cos(δφ)]. As we discussed in Sec. VII, single-photon loss of the selected mode leads to decoherence within the DFS, which limits the fidelity of prepared cat states. In contrast to the intrinsic photon loss discussed in Sec. VII, which is independent of γ, the single-photon loss rate γ φ is proportional to γ and thus it can reach significant values in the Zeno limit of large γ. The dependence of the unwanted singlephoton loss on γ leads to a trade-off between suppressing the excitations of remaining normal modes in the Zeno limit and keeping the unwanted single-photon loss rate small. As both of these competing processes decrease the fidelity of prepared cat states, γ has to be optimized to attain efficient cat-state preparation.
We also notice that the second smallest single-photon loss rate γ φ +2π/N = 2γ[1 − cos(2π/N − δφ)] decreases due to the offset δφ. Therefore, in order to achieve strong damping of the remaining normal modes, i.e. γ φ +2π/N η, G, we need not only a large rate of non-local dissipation γ but also a small offset δφ 2π/N , i.e., the non-local phase should be tunable with a sufficiently high resolution. The second condition becomes increasingly challenging for large system sizes N .
So far we have considered a common offset δφ for all the cavities, so that the overall translational invariance of the model is preserved. We now move to address the effects of disorder. To this aim, we consider a general form of imperfections in the engineered non-local dissipation L γ = j γ j D[â j − e φjâ j+1 ] leading to rates γ j and phases φ j of jump operators that vary at different bonds between resonators j and j + 1. These imperfections break the translational symmetry of the array. As a result, normal modeŝ b q = N j=1 T * qjâ j of the non-local dissipator L γ are no longer plane waves, where T jq is a transformation matrix and q = 1, 2, ...N [55]. As a consequence, pure steady states of the master equations (3) and (15) in the form of dark states do, in general, not exist. However, it is still possible to show that if all but one normal modes q = 1 experience large single-photon lossγ q =1 , cat states can be stabilized and prepared in the selected mode q = 1 in the Zeno limit. Here we ordered normal modes q by the loss rateγ q from smallest to largest, so that q = 1 is now the closest mode to achieve the ideal non-dissipative condition. For concreteness, we focus here on the alternative model (15) with local two-photon loss and with a vanishing two-photon pump phase θ = 0 [56]. In the Zeno limit γ q =1 η, G,γ 1 , the time evolution of mode q = 1 is described by an effective master equatioṅ T * 2 j1 andζ = −iNḠ/η (see Appendix G for more details). Despite the broken translational invariance, we still recover two-photon driven dissipation at rateη, which stabilizes a DFS spanned by cat states with amplitudesζ. We then get to a remarkable conclusion: in the Zeno limit, disorder affects only quantitatively the nature of two-photon driven dissipation, leading to a modified amplitude of cat statesζ = −iNḠ/η and of the dissipative gap ∆ d ≈ 2|Ḡ|. This provides an analytical insight into the robustness of the Zeno limit to disorder. However, from Eq. (18) we also see a qualitative change with respect to the disorderless case (16), namely the appearance of additive single photon loss term. Due to imperfections, non-local dissipation can give rise to single-photon loss of the selected mode q = 1. In particular, for a non-trivial overall phase Φ = N j=1 φ j = 2nπ, where n is any integer, the dark mode condition, T (j+1)1 = e φj T j1 for all j, is not compatible with periodic boundary conditions and, as a consequence, the selected mode exhibits a non-vanishing single-photon loss rateγ 1 = 0.
We now consider disorder in non-local dissipation rates γ j = γ(1 + σδγ j ) and jump-operator phases φ j = π 3 σδφ j , where δγ j and δφ j are normally distributed random variables and σ is the strength of disorder. The cat-state amplitudeζ and the dissipative gap ∆ d are only marginally modified in the presence of moderate disorder σ < ∼ 0.25 (see Appendix G for more details). On the other hand, the single-photon loss rate can reach considerable values even for moderate disorder, leading to significant decoherence within the DFS and, as a result, to a decreased fidelity of prepared cat states. The fidelity of prepared state is determined by the ratio of the decoherence rate |ζ| 2γ 1 and the dissipative gap ∆ d , see Sec. VII. We plot the distribution of the ratio in a disordered array as a function of the disorder strength in Fig. 5a. We can see that with the increasing strength of disorder it is more likely to obtain a decoherence rate |ζ| 2γ 1 ∼ ∆ d comparable to the dissipative gap ∆ d , which severely limits the fidelity of prepared cat states. On the other hand, the decoherence rate |ζ| 2γ 1 remains, on average (blue dashed line), order of magnitude smaller than the dissipative gap ∆ d even for moderate disorder σ ∼ 0.25, allowing the preparation of cat states with a large fidelity.
We now study the effects of disorder in arrays with increasing system size N . In order to keep the smallest loss rate of the remaining modesγ 2 ≈ γ approximately constant with increasing systems size, we increase the non-local dissipa- We plot the distribution of the ratio between the decoherence rate |ζ| 2γ 1 and the dissipative gap ∆ d for a constant disorder strength σ = 0.075 as a function of system size N in Fig. 5b. We can see that, with increasing system size, the ratio increases on average (blue dashed line) and it is more likely to obtain a decoherence rate |ζ| 2γ 1 ∼ ∆ d comparable to the dissipative gap ∆ d . The decoherence rate |ζ| 2γ 1 ∝ N increases due to the increasing catstate amplitude |ζ| ∝ √ N while the single-photon loss rateγ 1 and the dissipative gap ∆ d stay approximately constant with increasing system size (see Appendix G).
We note that for some disorder realizations the overall phase Φ is vanishingly small and, as a consequence, the single photon loss rateγ 1 is negligible. For these disorder realiza- . Single-photon amplification is a second-order perturbation process in the two-photon pumpĤ G and it is facilitated by virtual excitations of strongly decaying normal modes q = 1. However, single-photon amplification does not severely limit the fidelity of prepared cat states since it is suppressed in the Zeno limit as its rate is inversely proportional to γ (see Appendix H for more details).
In conclusion, due to imperfections in non-local dissipation, the selected normal mode is not a dark mode of the non-local dissipator. Remarkably, cat states can still be stabilized and prepared in the selected mode in the Zeno limit, albeit with decreased fidelity as decoherence due to unwanted single-photon loss sets in. Since the single-photon loss rate is proportional to the non-local dissipation rate γ, we face a trade-off between suppressing the excitations of remaining normal modes and keeping the decoherence rate small. The effects of decoherence are more pronounced for large system sizes due to the increasing cat-state amplitude. However, we can conclude that for moderate disorder strengths and moderate system sizes, the decoherence rate remains, on average, order of magnitude smaller that the dissipative gap, thus allowing for the preparation of cat states with a large fidelity even in the presence of imperfections.

IX. NOISE BIAS OF MULTI-MODE CAT STATES
A manifold spanned by Schrödinger cat states is advantageous for encoding and protecting quantum information because it experiences biased noise [9], which means that states in the manifold display an asymmetric response to different kinds of noise channels. This can considerably simplify the quantum-error-correcting protocols required to protect quantum information from decoherence. In this Section we first recall the most important features of biased noise and then show that multi-mode cat states offer an enhanced noise bias over single-mode cats considered so far.
For concreteness, let us assume we encode the logical quantum states into the following superposition of singlemode cat states, |0/1 L ≡ (|C + ± |C − )/ √ 2 ≈ | ± ζ (the last approximation is excellent for |ζ| > ∼ 2), as for instance realized with a single KPR cat state of amplitude ζ = i G/U [12]. The overlap between the coherent states | +ζ| − ζ | = exp(−2|ζ| 2 ) is exponentially suppressed with the number of photons |ζ| 2 . This simple fact endows the encoding with a natural protection against noise. Indeed, by increasing the amplitude, the coherent states move further apart and it becomes extremely unlikely for any noise process to cause a jump between |0 L and |1 L , i.e., to generate a socalled bit-flip error. Note that increasing the amplitude comes at the price of an increased phase-flip error rate, i.e., noiseinduced flips between |± L ≡ |C ± eigenstates. However, it can be shown that the phase-flip error rate increases only linearly with respect to the number of photons in the mode. As a result, there is a net bias exp(−2|ζ| 2 )/|ζ| 2 between the error rates experienced by a single resonator [9]. Such noise bias has been recently experimentally observed and tuned in Ref. [17].
Let us now consider a similar encoding, but replacing the single mode cat states with multi-mode cat states |C ± φ . This choice allows for a further exponential improvement in the noise bias, achieving exp(−2N |ζ| 2 )/(N |ζ| 2 ); the expression is written in terms of single-mode amplitude, for instance for the case of Kerr nonlinearity ζ N -mode ≡ i N G U = √ N ζ. The bias exponentially increases with respect to both the number of coherent photons in each resonator (set by the ratio between the two-photon pump and either the Kerr non-linearity or the two-photon loss rate) and the number of resonators in the array. While the first feature leads to 'standard noise bias' and is already exploited for qubit encoding in a single-mode cat manifold, the second feature is unique to our model. The physical reason for this is that the total number of photons in the multi-mode cat state stabilized in the resonator array increases with the number of resonators as discussed in Sec. III. Therefore, the multi-mode DFS (11) can encode a qubit with bit-flip errors exponentially suppressed with system size N , gaining an exponential improvement in the noise bias compared to single-mode cat states. This can be exploited to realize a protected quantum memory which only suffers from phase-flip errors, which can be in turn corrected by simple classical error correction techniques (e.g. a simple repetition code) [15].

X. EXPERIMENTAL PARAMETERS
For both models of cat-state preparation, the regime γ G, U, η κ is required to suppress excitations of remaining normal modes (γ G, U, η, as discussed in Sec. V) while at the same time reducing the effects of decoherence due to intrinsic single-photon loss (U, η κ, as discussed in Sec. VII). Strong Kerr nonlinearities can be effectively induced both in three-dimensional microwave cavities [54] and on-chip resonators [42], via coupling to a Josephson junction, resulting in photon-photon interactions far exceeding the photon loss rate. The induced Kerr nonlinearity U ∼ 100 − 1000 kHz can be larger than intrinsic loss rates κ ∼ 10 − 100 Hz of state-ofthe-art three-dimensional cavities [57,58] by orders of magnitude, with feasible ratios κ/U ∼ 10 −5 − 10 −3 . The stabilization of Schrödinger cat states in a single Kerr resonator using a parametric two-photon drive has been demonstrated in Ref. [7] reporting an amplitude |ζ| ∼ 1, which corresponds to G/U ∼ 1. In an array, the amplitude |ζ| ∝ √ N of the cat states increases with the number of Kerr resonators. Local two-photon loss can be realized in a similar three-dimensional cQED architecture [9] at rate η ∼ 100 kHz [10] with a feasible ratio κ/η = 10 −2 [14]. It has been employed for the stabilization of a single-mode Schrödinger cat state with an amplitude |ζ| ∼ 1 [10].
In addition to coupling neighboring resonators to a transmission line as discussed in [31], nonlocal dissipation can be also implemented by the coupling to a strongly decaying microwave mode with a feasible decay rate κ A ∼ 100 MHz [10]. Microwave modes can be coupled using parametrically driven elements consisting of several Josephson junctions reaching coupling strengths g ∼ 100 MHz as demonstrated in planar superconducting circuits [59,60]. Using a weaker parametric drive such that g κ A , one can achieve effective nonlocal dissipation at rate γ ∼ g 2 /κ A by adiabatically eliminating the strongly decaying mode. It is also feasible to achieve the Zeno limit of strong non-local dissipation γ G, U, η as the parametric coupling g and the decay rate κ A can be sufficiently large.
Finally, we stress that, while multi-mode cat states can be stabilized and prepared in resonator arrays with a large system size N , all characteristic features of the multi-mode cat-state preparation discussed in this manuscript-and especially the exponentially enhanced noise bias-can be observed already for a small system size N = 3.

XI. CONCLUSIONS
We showed that a two-dimensional manifold spanned by superpositions of multi-mode Schrödinger cat states can be stabilized in an array of resonators coupled via non-local dissipation. The required non-linearity, which is either of the Kerr type or an engineered two-photon loss, acts locally on each resonator while the dissipative coupling is linear, which makes our proposal particularly convenient for experimental implementations. The two models we put forward are readily realizable with state-of-the-art circuit QED architectures.
In the Zeno limit of strong non-local dissipation, we showed that the even-parity multi-mode cat state can be prepared from the initial vacuum state with a fidelity approaching unity. In the Kerr-resonator array, the steady-state preparation is limited by decoherence due to intrinsic photon loss, which sets in at long times, as the relaxation towards the steady state is due to a weak effective two-photon dissipation resulting in a slow rate of convergence. On the other hand, local two-photon loss gives rise to a quick convergence towards the multi-mode cat states allowing for their efficient preparation. Multi-mode cat states can be stabilized and prepared in the Zeno limit even in the presence of imperfections in non-local dissipation, albeit with decreased fidelity as an additional decoherence channel sets in. Importantly, the rate of convergence towards the cat states is independent of system size making the steady-state preparation scalable. Scaling towards large system sizes requires a good suppression of imperfections in nonlocal dissipation as they lead to decoherence which is more pronounced in large systems. Our protocol exploits N local two-photon pumps allowing for the amplitude |ζ| ∝ √ N of the multi-mode cat states to increase with the number of resonators in the array N . This in turn leads to an exponentially increasing noise bias of the cat manifold, which can be exploited for the implementation of a quantum memory with enhanced protection.
Being able to prepare multi-mode cat states in a deterministic, robust and scalable fashion is of direct relevance for applications in quantum metrology, quantum computation and quantum information. In quantum computation, especially in the context of bosonic codes based on cat qubits, multimode cat states correspond to GHZ states, which are a required resource for universal fault-tolerant quantum computing, e.g. by enabling Toffoli state preparation [6]. In quantum metrology, multi-mode cat states can provide a source of non-Gaussian multimode light needed for attaining Heisenberg scaling [22,23]; they are also a source of non-Gaussian multipartite entanglement to generate non-Gaussian cluster states [61]. Moreover, our proposal may also foster multimode extensions of existing protocols currently limited by the lack of non-Gaussian resources beyond the single-or twomode case. Multi-mode cat states are also relevant to fundamental aspects of quantum theory, e.g. to test predictions of decoherence theory and explore the quantum-to-classical transition [62].
Our treatment of dissipatively coupled cavity arrays opens new avenues in quantum reservoir engineering [30] as it allows for tailored dissipation in momentum space. Interesting future directions are the stabilization of other non-Gaussian multi-mode entangled states using different on-site models, reservoir engineering with multiple normal modes with different quasi-momenta by employing different non-local reservoirs, and tailoring dissipation in momentum space of higher dimensions. Finally, for the specific models considered in this work, even though the steady states can be described analytically for any system parameters and the transient dynamics can be well described in the Zeno limit by an effective master equation for a single normal mode, much less has been understood about the driven-dissipative many-body dynamics for the non-local dissipation rate comparable to the Kerr nonlinearity, which will be subject of future investigation.

NOTE ADDED
During the final stage of this project, a related paper appeared, investigating mixed steady states of dissipatively coupled Kerr resonator arrays, which correspond to the ground states of a frustrated antiferromagnet [63]. After the completion of this project, another work appeared, studying the generation and detection of two-mode cat states in dissipativelycoupled parametric oscillators [64], however only treating the particular case N = 2.
Appendix A: Numerical simulations of the master equation (3) In this appendix we provide details about numerical simulations of the master equation (3) and its conserved quantities.
In general, steady states which are neither dark states nor their incoherent mixtures can exist [40]. However, we verified that the master equation (3) has no steady states lying outside of the DFS (11) by numerically solving for Lρ ss = 0 for N = 3, 4 and various values of the non-local decay rate γ/U and two-photon drive strength G/U . We truncate the infinitely dimensional Hilbert space considering only M k lowest energy levels of each mode k.
Each coefficient c µ = Tr ρ inĴµ of the DFS (11) is associated with a conserved quantityĴ µ and encodes the information aboutρ in that is preserved during the time evolution. As a result, conserved quantities can be used to find the particular steady state for a given initial state. Conserved quantitiesĴ µ are formally defined as solutions of L †Ĵ µ = 0, where L † is the adjoint Liouvillian and and {·, ·} is the anticommutator [50]. We solve numerically for the conserved quantities. Finding conserved quantities, which are bi-orthogonal to the basis operatorsξ ±± = |C ± C ± | of the DFS, i.e. Tr[Ĵ † µξβ ] = D µ δ µβ , we can determine coefficients c µ = Tr[Ĵ † µρin ] of the particular steady state for a given initial stateρ in . Normalization D µ = 1, for all µ, guarantees thatρ is a density matrix. We plot in Fig. 3 the fidelities with the even cat state and the odd cat state as well as the purity of the particular steady state for the initial vacuum state, where we used the following truncation of the Hilbert space: M k=φ = 12 and M k =φ = 3 for G/U = 0.5, M k=φ = 16 and M k =φ = 3 for G/U = 0.75, and M k=φ = 18 and M k =φ = 3 for G/U = 1.

Appendix B: Effective master equation in the Zeno limit
In this appendix, we describe the effective Zeno dynamics of mode k = φ in the regime of strong non-local dissipation, which acts as a continuous measurement projecting all modes k = φ onto the vacuum state. We derive an effective master equation for mode φ employing the Dyson series of the Liouvillian dynamics, which was described in detail in Ref. [53].
We start by rescaling time τ = γ t in the master equation (3) obtaining In the limit of strong dissipation γ k U, G, for all k = φ, we can treat K as a perturbative term. A general solution of equation (B1) can be written asρ where U(τ ) is the propagator satisfying By iterating this equation, we obtain the Dyson series for the propagator where the ellipsis denotes terms of third and higher order in K.
Mode φ does not experience any dissipation (Tr dLd )ρ φ = ρ φ , whereρ φ = Tr dρ is the reduced density matrix of mode φ and Tr d denotes the trace over all decaying modes k = φ. The dissipatorL d targets a unique stateρ d in the reduced Hilbert space of all decaying modes k = φ, i.e.
where Tr φ is the trace over mode φ. The unique steady state is the vacuum stateρ d = k =φ |0 k 0|. The projection onto the kernel ofL d is P d = lim τ →∞ exp(L d τ ) obeying following relationsL and for an arbitrary density matrix or operatorX. For strong dissipation γ k U, G, the dynamics of the system is constrained to the dissipation-free subspaceρ(τ ) = ρ d ⊗ρ φ (τ ) at all times τ 1 [53]. The dynamics of the reduced density matrixρ φ (τ ) is described by the propagator P d U(τ ) P d . Writing the Dyson series for the propagator up to the second order in K, we obtain where we used equation (B6). We explicitly evaluate the second-order term in the Dyson series in the Supplementary Material to obtain We now consider that we initially start with a state in the dissipation-free subspace, i.e.ρ(0) =ρ d ⊗ρ φ (0). The time evolution within the dissipation-free subspace is described bŷ To obtain a master equation in a differential form, we use the equivalence ∂ρ φ /∂τ = lim τ →0 {ρ φ (τ ) −ρ φ (0)} /τ . Using the Dyson series (B9) and rescaling back time t = τ /γ, we obtain the effective master equation (13) in Lindblad form. The effective HamiltonianĤ φ = Tr d ρ d ⊗Î φ Ĥ is the projection of the full HamiltonianĤ onto the dissipation-free subspace, whereÎ φ is the identity operator in the reduced Hilbert space of mode φ.

Appendix C: Rotation within the stabilized DFS
We consider a weak single-photon driveĤ λ = λ(e ijφâ j + e −ijφâ † j ) applied to an arbitrary single resonator j. In the Zeno limit γ U, G, λ, the dynamics of mode φ is described by the effective master equatioṅ To derive the effective master equation, we followed the same steps as in Appendix B for the Kerr resonator array without the single-photon drive. The single-photon driveĤ λ of resonator j gives rise to the single-photon driveĤ φ, λ = λ √ N (b φ +b † φ ) of mode φ, which is a first order process inĤ λ . Note that the single-photon driveĤ λ in combination with the Kerr non-linearityĤ U and the two-photon driveĤ G does not lead to any second-order effect due to the strong suppression of exactions of decaying modes and the conservation of total quasimomentum.
It has been shown in Ref. [9] that the single-photon drivê H φ, λ induces the rotation within the DFS stabilized by two-photon driven dissipation , provided that the drive strength λ/ √ N is small compared to the rate Γ of two-photon driven dissipation. The even cat state |C + , which is dissipatively prepared from the initial vacuum state, can be transformed into a superposition R(α)|C + = cos α|C + + i sin α|C − of cat states |C + and |C − by the rotationR(α).
Appendix D: Decoherence-free subspace of the alternative model (15) In this appendix, we discuss the steady states of the alternative model with local two-photon los described by the master equation (15).
We start by writing the dissipator in a form of the non-Hermitian HamiltonianĤ eff and the jump term Jρ which can be expressed in the plane wave basis aŝ where the arguments k 1 + k 2 and k 3 + k 4 are defined modulo 2π. Similarly to the Kerr-resonator array, we look for pure steady statesρ SS = |Ψ Ψ| with all modes k = φ being in the vacuum state, in which case the non-local dissipators k γ k D[b k ]ρ SS = 0 in Eq. (15) vanish. Usingb kρSS =ρ SSb † k = 0 for all k = φ, we evaluate all remaining terms in the master equation (15) obtaininĝ All three remaining terms vanish if the pure steady state is one of the coherent statesb φ |ψ =b φ | ±ζ = ±ζ|ζ or an arbitrary superposition of these coherent states. We conclude that an arbitrary superposition of coherent states | ±ζ is a steady state of the master equation satisfying Lρ SS = 0. Similarly as for the array of Kerr resonators, all possible superpositions of coherent states | ±ζ and their incoherent admixtures form the DFS (11). We numerically confirm that there are no steady states lying outside the DFS (11).
Appendix E: Effective master equation in the Zeno limit for the alternative model (15) In this appendix, we derive the effective master equation for mode φ for the alternative model (15) with local two-photon loss in the Zeno limit of strong non-local dissipation. We follow the same steps in the derivation as for the Kerr-resonator array in Appendix B.
We start by rescaling time τ = γ t in the master equation (15) obtaining where In the limit of strong non-local dissipation γ k η, G, for all k = φ, we can treat K as a perturbative term. A general solution of the master equation isρ(τ ) = U(τ )ρ(0), where the propagator U(τ ) can be expanded in the Dyson series (B4).
The dissipatorL d is identical to the dominant dissipator for the Kerr-resonator array in the Zeno limit. It targets a unique stateρ d in the reduced Hilbert space of modes k = φ. For strong non-local dissipation γ k η, G, the dynamics of the system is constrained to the subspaceρ(τ ) =ρ d ⊗ρ φ (τ ) at all times τ 1 [53]. The dynamics of the reduced density matrixρ φ (τ ) is described by the propagator We used the Dyson series (B4) as well as the relation (B6) in the first equality and, in the second equality, we explicitly evaluated the first-order Note that the first-order term in the Dyson series leads to decay towards the DFS (11) at times t ∼ N η [32]. As a result, it is sufficient to consider the Dyson series only up to the first order in K as higher-order terms set in at longer times t ∼ γ η 2 . This is in contrast to the Zeno dynamics of the Kerr-resonator array, where the first-order term in the Dyson series leads to the unitary dynamics and decay towards the DFS is due to the second-order term. We now consider that we initially start with a state in the subspaceρ d ⊗ρ φ . The time evolution within this subspace is described byρ φ (τ ) = Tr d [P d U(τ ) P d (ρ d ⊗ρ φ (0))]. We obtain the effective master equation (16) in Lindblad form where we used the Dyson series (E2) and rescaled back time t = τ /γ.

Appendix F: Dissipative gap in the Zeno limit
In this appendix we study the dissipative gap ∆ d in the Zeno limit, which is determined from the spectrum of the effective Liouvillian L φ as the smallest non-vanishing real part of the eigenvalues [50]. We investigate the dissipative gap first for the Kerr-resonator array and then for the alternative model with local two-photon loss.
For the Kerr-resonator array, we numerically determine the spectrum of the Liouvillian, and from that we extract the dissipative gap. We plot the dissipative gap in Fig. 6a as a function of the non-local decay rate γ for several values of the twophoton pump strength G. From the numerical data, we infer that the dissipative gap is ∆ d = Γ 2 N U 1 , where 1 is the energy of the first excited state of the effective HamiltonianĤ φ . For large |ζ| > ∼ 1, the energy of the first excited state approaches 1 ≈ 4 U N |ζ| 2 and, as a result, the dissipative gap can be ap-proximated by Eq. (14). In Fig. 6a, we can see that the exact values (solid lines) approach the values given by Eq. (14) (dashed lines) as the strength G of the two-photon drive and, as a consequence, |ζ| increase. We now study the dissipative gap ∆ d for the alternative model, which is extracted from the numerically calculated spectrum of the Liouvillian L φ . Note that, in contrast to the Kerr-resonator array, the Liouvillian and, as a consequence, the dissipative gap are independent of γ in the Zeno limit for large γ. We plot the dissipative gap ∆ d in Fig. 6b as a function of the pump power G. It shows that the dissipative gap grows as ∆ d ≈ 2η|ζ| 2 /N with the amplitude ζ.
Appendix G: Effective master equation (18) in the presence of imperfections In this appendix, we discuss the effective master equation for mode q = 1 for the alternative model (15) with local twophoton loss in the Zeno limit of strong non-local dissipation in the presence of imperfections.
To derive the effective master equation in the Zeno Starting initially in the subspacê ρ d ⊗ρ 1 , the dynamics of the reduced density matrixρ 1 is described by the effective master equation (18) in Lindblad form within first-order perturbation theory. In the derivation of the effective master equation, we followed the same steps as in Appendix E for the alternative model without imperfections.
The amplitude of cat statesζ and the dissipative gap ∆ d are modified due to disorder in non-local decay rates γ j and jump-operator phase φ j . We plot the distribution of the catstate amplitude and the dissipative gap in Figs. 7a and 7b, respectively, in the presence of random disorder as a function of the disorder strength σ. We can see that for moderate disorder σ < ∼ 0.25, the cat-state amplitude and the dissipative gap retain values that are comparable to those of the disorder-free array for σ = 0. This is in contrast to the single-photon loss rateγ 1 , which reaches considerable values due to even moderate disorder σ ∼ 0.25 and leads to a significant decoherence rate |ζ| 2γ 1 in comparison to the dissipative gap as shown in Fig. 5 in the main text.
We plot the distribution of the single-photon loss rateγ 1 and the dissipative gap ∆ d in Figs. 7c and 7d, respectively, for the constant disorder strength σ = 0.075 as a function of system size N . We can see that the dissipative gap moderately decreases on average (blue dashed lines) with increasing system size but it retains values that are comparable to the value ∆ d ≈ 2G of the disorder-free array. The single-photon loss rate moderately increases on average (blue dashed lines) with increasing system size. This is in contrast to the cat-state amplitude |ζ| 2 ∝ N which increases proportionally to system size N leading to a rapidly increasing decoherence rate |ζ| 2γ 1 as discussed in the main text. In this appendix, we discuss single-phonon amplification which is a second order process inĤ G in the Zeno limit. Crucially, single-photon amplification is the dominant source of decoherence within the DFS for some realizations of disorder in non-local dissipation. As a result, single-photon amplification has to be taken into account to accurately describe the dynamics of the selected mode in the Zeno limit.
In particular, we consider the alternative model with vanishing two-photon pump phase θ = 0 [56] and with disordered non-local dissipation rates γ j and jump-operator phases φ j . Within second-order perturbation theory, the dynamics of the selected mode q = 1 in the Zeno limit is described by the effective master equatioṅ Single-photon amplification at rate is a second order process in the two-photon pumpĤ G and is absent within the first order perturbation theory, see Eq. (18).
To derive the master equation (H1), we evaluate the second order term in the Dyson series (B8) where we consider only the perturbation in a form of the twophoton pump K = − i γ [Ĥ G , · ]. We neglect the remaining perturbations as they give rise to only marginal second-order corrections.
Single-photon amplification is facilitated by virtual excitations of strongly decaying modes q = 1. For an ideal array with vanishing jump-operator phases φ j = 0, the simultaneous excitation of the selected mode q = 1 and another mode q = 1 via the two-photon pump is not allowed due to the conservation of quasi-momentum. As a result, the single-photon amplification rate µ vanishes. On the other hand, imperfections in the form of non-vanishing jump-operator phases φ j = 0 lead to the breaking of the quasi-momentum conservation. As a consequence, these imperfections enable the simultaneous excitation of the selected mode q = 1 and another mode q = 1 giving rise to non-vanishing single-photon amplification.
We will now show that for imperfections leading to vanishingly small overall phase Φ = N j=1 φ j , single-photon amplification is a dominant source of decoherence as single-photon loss rateγ 1 is negligible. We consider a disordered array with the vanishing overall phase Φ = 0 and, as a consequence, γ 1 = 0. In Fig. 8a we plot the fidelity c ++ of the instantaneous stateρ(t) with the even cat state |C + (blue lines) and the fidelity c −− with the odd cat state |C − (green lines) as a function of time for different values of the non-local dissipation rare γ and starting initially from the vacuum state. In Fig. 8b, we plot the purity p of the instantaneous state. We can see that at long times the array approaches the fully mixed statê ρ 1 ≈ 1 2 (|C + C + | + |C − C − |) with c ++ ≈ c −− ≈ p ≈ 0.5 as decoherence due to single-photon amplification sets in. We compare full master equation (15) simulations (dark blue and dark green lines) to the simulations of the effective master equation (H1) in the Zeno limit (light blue and light green lines) and we plot in Fig. 8b difference ρ −ρ d ⊗ρ 1 (purple lines) between the instantaneous stateρ according to the full master equation and the instantaneous stateρ 1 according to the master equation (H1) for the effective Zeno dynamics. We can see that in the Zeno limit for large non-local dissipation rate γ, the effective master equation (H1) describes well the time-evolution of the array as the difference ρ −ρ d ⊗ρ 1 is small throughout the entire time evolution.
Single-photon amplification excites the array out of the DFS. As single-photon amplification µ is a second-order processes, it is, in the Zeno limit, week compared to the twophoton driven dissipationη, which stabilizes the DFS. As a result, upon excitation out of the DFS the array quickly converges back towards the DFS. However, the single-photon amplification process transforms odd-parity eigenstates to evenparity eigenstates and vice versa as it injects a single photon into the array. As a result, single-photon amplification leads to leaking between parity subspaces and convergence towards the fully mixed state.
Even though, single-photon amplification leads to decoherence within DFS, it does not severely limit the fidelity of prepared cat states in the Zeno limit. Since single-photon amplification is a second-order process facilitated by virtual excitations of strongly decaying modes, its rate is inversely proportional to γ, see Eq. (H2). As a result, single-photon amplification is suppressed in the Zeno limit. This is in contrast to single-photon loss whose rate is proportional to γ leading to strong decoherence in the Zeno limit. [48] The joint Wigner function of the N -mode stateρ is W ( x, p) = π −N R N d y exp(−2 i y · p) x + y|ρ| x − y , where x, p and y are real-valued N -dimensional vectors, and xj and pj are quadratures of resonator j. For cat states stabilized in mode φ = 0, the amplitudes ζj and ζj+1 at neighboring resonators have a relative phase difference φ [see Eq. (12)]. This results in a relative rotation between the local coherent states |ζj j and |ζj+1 j+1 in phase space by an angle φ.
[49] Note that the dynamics of an array with only two resonators, i.e. N = 2, is fundamentally different from that of larger arrays with N ≥ 3 as the photon parityP φ = exp iπb † φb φ of mode φ is conserved, where the quasi-momentum can take only values φ = 0, π. As a result, the even-parity cat state |C + is deterministically approached from the initial vacuum state for any γ. However, in this manuscript we focus on arrays N ≥ 3 for which the photon paritiesP andP φ are not conserved.