$\dagger$$\dagger$footnotetext: These authors contributed equally.

# Candidate for a passively protected quantum memory in two dimensions

Simon Lieu† Joint Quantum Institute, NIST/University of Maryland, College
Park, Maryland 20742, USA  Joint Center for Quantum Information and Computer
Science, NIST/University of Maryland, College Park, Maryland 20742 USA  Yu-Jie
Liu† Department of Physics, Technical University of Munich, 85748 Garching,
Germany  Munich Center for Quantum Science and Technology (MCQST),
Schellingstr. 4, 80799 München, Germany  Alexey V. Gorshkov  Joint Quantum
Institute, NIST/University of Maryland, College Park, Maryland 20742, USA
Joint Center for Quantum Information and Computer Science, NIST/University of
Maryland, College Park, Maryland 20742 USA

###### Abstract

An interesting problem in the field of quantum error correction involves
finding a physical system that hosts a “passively protected quantum memory,”
defined as an encoded qubit coupled to an environment that naturally wants to
correct errors. To date, a quantum memory stable against finite-temperature
effects is only known in four spatial dimensions or higher. Here, we take a
different approach to realize a stable quantum memory by relying on a driven-
dissipative environment. We propose a new model, the photonic-Ising model,
which appears to passively correct against both bit-flip and phase-flip errors
in two dimensions: A square lattice composed of photonic “cat qubits” coupled
via dissipative terms which tend to fix errors locally. Inspired by the
presence of two distinct $\mathbb{Z}_{2}$-symmetry-broken phases, our scheme
relies on Ising-like dissipators to protect against bit flips and on a driven-
dissipative photonic environment to protect against phase flips. We also
discuss possible ways to realize the photonic-Ising model.

Quantum error correction remains one of the biggest challenges towards
building a practical quantum computer [1, 2]. One of the leading candidates
for realizing fault tolerance is the family of quantum stabilizer codes [3],
including the surface code [4, 5, 6] and the GKP code [7]. These error-
correcting schemes are based on fast error recovery controlled by the feedback
from repetitive syndrome measurements.

A prominent alternative is the finite-temperature quantum memory: Certain
thermal environments naturally evolve arbitrary initial states into a qubit
subspace of interest at low temperature, thus eliminating the need for active
measurements and correcting operations. Many recent studies have investigated
thermal self-correcting properties [8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18,
19, 6]. To date, the only known models that host a passive quantum memory via
this mechanism are topological codes in four dimensions (4D) and higher, e.g.
the 4D toric code [19, 6].

A separate line of research aims to uncover a passively protected quantum
memory via engineered “driven-dissipative” systems [20, 21, 22, 23, 24, 25,
26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39]. Such passive
protection includes but is not limited to the finite-temperature case, since a
thermal-equilibrium steady state is not required. The memory is dynamically
protected against certain noise channels by (local) Markovian dissipation.
This has led to a number of new ideas for passive error correction, such as
the autonomously corrected cat qubit [40, 41] and the dissipative Toom’s rule
[42, 27, 6]. Unfortunately, none of these models can protect a quantum memory
for an exponentially-long time as a function of the system size (in less than
four dimensions).

In this work, we study a model with engineered dissipation which appears to
protect against both bit flips and phase flips and lives in two spatial
dimensions. Instead of relying on topological order, we suggest that the model
should belong to a phase that spontaneously breaks two different
$\mathbb{Z}_{2}$ symmetries. Each $\mathbb{Z}_{2}$-symmetry-broken phase
protects a “classical bit,” which together form a robust qubit. The proposed
model provides an example of a robust quantum memory which, at low
temperature, can be exponentially long-lived in system size parameters and has
challenging, yet realistic physical requirements.

Quantum memory.—Consider a Hilbert space $\mathcal{H}$, and define two
encoded, logical states $|\bar{0}\rangle,|\bar{1}\rangle\in\mathcal{H}$ that
span the codespace $\mathcal{C}$. We assume the system is always initilized in
the codespace: $\rho_{i}=|\psi\rangle\langle\psi|$ where
$\ket{\psi}\in\mathcal{C}$.

A local continuous-time Markovian generator $\mathcal{L}$ in Lindblad form is
defined by

$\displaystyle\frac{d\rho}{dt}=\mathcal{L}(\rho)=-i[H,\rho]+\sum_{j}\left(L_{j}\rho
L_{j}^{{\dagger}}-\frac{1}{2}\\{L^{{\dagger}}_{j}L_{j},\rho\\}\right),$ (1)

where $H$ is the Hamiltonian of the system and $L_{j}$ are local dissipators
which arise due to the system-environment coupling [43]. We consider a
dynamical process that can be decomposed into two parts, an “error” generator
and a “recovery” generator: $\mathcal{L}=\mathcal{L}_{e}+\mathcal{L}_{r}$. The
error generator describes the main channels of physical noise which move the
initial state out of the codespace. The recovery generator stabilizes the
codespace: $\mathcal{L}_{r}(\rho_{i})=0$, i.e. any state in the codespace is a
steady state of the recovery. We allow for this noisy process to occur for a
time $t$, which generically sends $\rho_{i}$ to a mixed state
$\rho_{m}(t)=e^{\mathcal{L}t}(\rho_{i})$.

Finally, we employ a “single-shot” decoding quantum channel $\mathcal{E}_{r}$
which sends every state in the Hilbert space back to the codespace [44]. The
final state is

$\rho_{f}(t)=\mathcal{E}_{r}e^{\mathcal{L}t}(\rho_{i}).$ (2)

We wish to find systems where the difference between the initial and final
states is exponentially small in the system size:

$\displaystyle 1-\Tr[\rho_{i}\rho_{f}(t)]=O(e^{-\gamma M})\text{ as
}M\to\infty,$ (3)

where $\gamma>0$ is a time-independent constant and $M$ is some system size
parameter. A system described by $\mathcal{L}$ hosts a passively protected
quantum memory for any finite time $t$ if Eq. (3) holds as the thermodynamic
limit is approached.

The bit-flip and phase-flip errors of a two-level system are generated via the
Pauli operators $X,Z$ respectively. A good quantum memory should thus protect
against both sources of noise. Recent work [39] has described the connection
between $\mathbb{Z}_{2}$ symmetry breaking and error correction: A symmetry-
broken phase protects quantum information against $X$ or $Z$ errors, but not
both. This leads to a protected classical bit, which can be viewed as a
quantum bit experiencing biased noise [45].

In this work, we attempt to glue two different classical bits together to form
a robust qubit. Our strategy involves studying a system that passively
corrects against bit flips due to Ising-like dissipators which tend to align
qubits locally. Furthermore, phase flips will passively correct due to driven-
dissipative stabilization of the photonic cat code. We begin by describing
spontaneous symmetry breaking in the cat code and in the Ising model
separately. We then describe a model which inherits both protecting features.

Photonic cat code.—Let us briefly review $\mathbb{Z}_{2}$ spontaneous symmetry
breaking in the photonic cat code [28, 46]. [For a detailed analysis, we refer
to Ref. [39].] Consider a driven-dissipative photonic cavity in the presence
of two-photon drive and two-photon loss. The rotating-frame Hamiltonian and
dissipator read $H=\lambda\left(a^{2}+(a^{\dagger})^{2}\right),\
L_{2}=\sqrt{\kappa_{2}}a^{2}.$ Here $a$ is the annihilation operator for a
cavity photon, $\lambda$ is the drive strength, and $\kappa_{2}$ is the two-
photon loss rate. While the model has $\mathbb{Z}_{2}$ symmetry
$[H,Q]=[L_{2},Q]=0$ generated by parity $Q=e^{i\pi a^{\dagger}a}$, the steady
state can violate this symmetry:

$\displaystyle\rho_{ss}=|\psi\rangle\langle\psi|,\qquad|\psi\rangle=c_{0}|\alpha_{e}\rangle+c_{1}|\alpha_{o}\rangle,$
(4)

for $|c_{0}|^{2}+|c_{1}|^{2}=1$, where
$|\alpha_{e}\rangle\sim|\alpha\rangle+|-\alpha\rangle,|\alpha_{o}\rangle\sim|\alpha\rangle-|-\alpha\rangle$,
and $|\alpha\rangle$ is a coherent state with amplitude
$\alpha=e^{-i\pi/4}\sqrt{N}$ and $N\equiv\lambda/\kappa_{2}$ photons. The even
and odd cat states $|\alpha_{e/o}\rangle$ represent logical 0 and 1,
respectively.

The cat code is protected against phase-flip errors generated by photon
dephasing $L_{d}=\sqrt{\kappa_{d}}a^{\dagger}a$. Indeed, the phase-flip
logical error rate scales as $e^{-\gamma N}$ where $\gamma$ is a constant
[28]. The symmetry-broken states
$|\pm\alpha\rangle\approx(|\alpha_{e}\rangle\pm|\alpha_{o}\rangle)/\sqrt{2}$
have an exponentially-long lifetime in the limit of large $N$, ensuring that
logical phase flips are unlikely.

The dominant decoherence mechanism for the cat qubit stems from the bit flip,
generated via single-photon loss $L_{1}=\sqrt{\kappa_{1}}a$:
$a|\alpha_{e/o}\rangle\sim|\alpha_{o/e}\rangle$, which reduces the qubit
steady state structure to a classical bit: $\rho_{ss}\approx
c|+\alpha\rangle\langle+\alpha|+(1-c)|-\alpha\rangle\langle-\alpha|,c\in[0,1]$
[39]. More generally, perturbations that commute with photon parity (e.g.
$[L_{d},Q]=0$) are expected to be passively corrected, while terms which
explicitly break the symmetry (e.g. $\\{L_{1},Q\\}=0$) are not.

2D Ising model.—We now turn our attention to a system that has the opposite
problem: $\mathbb{Z}_{2}$ symmetry breaking will protect against bit flips but
not phase flips. We consider the 2D Ising model on an $M\times M$ square
lattice with periodic boundary conditions. The Hamiltonian reads

$\displaystyle H_{is}=-\sum_{x,y=1}^{M}(Z_{x,y}Z_{x+1,y}+Z_{x,y}Z_{x,y+1})\,,$
(5)

where $Z_{x,y}$ is the $Z$ Pauli operator on site $(x,y)$. The ferromagnetic
states are the ground states of this model and span the codespace:
$|\bar{0}\rangle\equiv|\downarrow\downarrow\downarrow\ldots\rangle,|\bar{1}\rangle\equiv|\uparrow\uparrow\uparrow\ldots\rangle$,
with $Z\ket{\downarrow}=\ket{\downarrow}$ and
$Z\ket{\uparrow}=-\ket{\uparrow}$.

We define local dissipators that describe the thermalization of the Ising
Hamiltonian. (For simplicity, we set the Hamiltonian in the master equation to
zero.) Consider dissipators that are a product of a spin flip ($X$) with a
projector onto a particular domain-wall configuration. These jumps will cause
a spin to flip sign according to a local “majority rule,” i.e. only if more
than two of the neighboring spins are misaligned. Specifically:

$\displaystyle L_{x,y}^{(4)}$
$\displaystyle=\sqrt{\kappa}X_{x,y}P_{x,y;\rightarrow}^{-}P_{x,y;\uparrow}^{-}P_{x-1,y;\rightarrow}^{-}P_{x,y-1;\uparrow}^{-}\,,$
$\displaystyle L_{x,y}^{(3)}$
$\displaystyle=\sqrt{\tilde{\kappa}}X_{x,y}P_{x,y;\rightarrow}^{+}P_{x,y;\uparrow}^{-}P_{x-1,y;\rightarrow}^{-}P_{x,y-1;\uparrow}^{-}\,,$
(6)

where $\tilde{\kappa}=\sqrt{\Delta\kappa+\Delta^{2}}-\Delta$ and
$P_{x,y;\rightarrow}^{\pm}=(1\pm Z_{x,y}Z_{x+1,y})/2$ ,
$P_{x,y;\uparrow}^{\pm}=(1\pm Z_{x,y}Z_{x,y+1})/2$ are projectors onto
particular local configurations of spins. The superscripts indicate the number
of domain walls which the projector is checking for, and we neglect to write
jumps related by rotational invariance (i.e. there are 4 different $L^{(3)}$
operators per site) [47]. We also consider an error process in the form of a
uniform bit flip rate on each lattice site:
$L_{x,y}^{\prime}=\sqrt{\Delta}X_{x,y}.$.

Figure 1:  The total rate of transitioning from a configuration with 4 domain
walls to a configuration with 0 domain walls satisfies detailed balance:
$\kappa_{4\rightarrow 0}/\kappa_{0\rightarrow 4}=e^{8\beta}$.

We have chosen our dissipators above such that the steady state of the model
is the thermal state of the 2D classical Ising model:

$\displaystyle\rho_{ss}=\frac{e^{-\beta H_{is}}}{\Tr[e^{-\beta
H_{is}}]},\qquad\beta=\frac{1}{8}\ln\left[\frac{\kappa+\Delta}{\Delta}\right],$
(7)

with the effective temperature set by the relative ratio of the correction
rate to the bit-flip rate. Within the quantum jump picture [48, 49], the rates
of transitioning between different classical configurations respect detailed
balance. (See e.g. Fig. 1.)

While the thermal state (7) is always a steady state of the model, it is not
unique. All dissipators commute with the parity operator
$Q=\prod_{i=1}^{M^{2}}X_{i}$: $[L_{j},Q]=0$. This means that the dynamics
preserves the parity of the state (called a “strong $\mathbb{Z}_{2}$ symmetry”
[50]). In the thermodynamic limit of the low-temperature (symmetry-broken)
phase, a qubit can be stored in the steady state [39].

We can confirm this picture via numerical simulations. Suppose we initialize
our system in a ferromagnetic state:
$|\psi\rangle=|\bar{0}\rangle=(|E_{0}^{+}\rangle+|E_{0}^{-}\rangle)/\sqrt{2}$
where $|E_{0}^{\pm}\rangle$ are ground states in the different parity sectors
[51]. We then quench the system with the noisy Lindbladian for a time $T$ much
larger than the inverse of the dissipative gap, so that the system settles
into its steady state. Finally, we apply a single-shot decoder which brings
the state back to the codespace by measuring all domain walls in the system
then flipping all bits in the smaller domain. Our results are summarized in
Fig. 2. In the low-temperature phase, the overlap starts to approach the ideal
value of 1 exponentially fast in $M$. Qualitatively different behavior occurs
in the high-temperature phase [$\beta>\beta_{c}=\ln(1+\sqrt{2})/2\approx
0.44$; red dots], where the success rate stays at 50% for a wide range of $M$.

Figure 2:  (a) The overlap between the initial and final states for the
protocol given in the main text, for a Lindbladian in the high-temperature
phase (red dots), and in the low-temperature phase (black and blue dots). As
linear system size $M$ grows, the overlap approaches one only in the low-
temperature (symmetry-broken) phase corresponding to $\beta>\beta_{c}\approx
0.44$. (b) Same black data points on a log plot; the overlap tends to one
exponentially fast in $M$. In both (a) and (b), the quench time is
$T=800/\kappa$, i.e. long enough to reach the steady state. The simulation
employs the quantum jump approach by averaging over $10^{5}$ trajectories.

Unfortunately, the stored qubit is unstable to noise that violates the strong
symmetry. In particular, the presence of dephasing $L_{i}\sim Z_{i}$ (phase
flips), reduces the strong $\mathbb{Z}_{2}$ symmetry to a “weak
$\mathbb{Z}_{2}$ symmetry” (defined at the level of the superoperator:
$[\mathcal{L},\mathcal{Q}]=0$, where $\mathcal{Q}(\rho)=Q\rho Q^{{\dagger}}$),
such that only a classical bit can be stored in the steady state. In this
case, the steady state at low-temperature has the structure $\rho_{ss}\approx
c|\bar{0}\rangle\langle\bar{0}|+(1-c)|\bar{1}\rangle\langle\bar{1}|$, for
$c\in[0,1]$. In analogy with the cat qubit in the presence of single-photon
loss, $Z$ dephasing destroys the coherence between Ising ferromagnetic states.

2D photonic-Ising model.—We see that the cat code passively corrects against
phase flips but not bit flips, and that the 2D Ising model passively corrects
against bit flips but not phase flips. Is it possible to combine the
protecting features of both models to construct a system that passively
corrects against both sources of noise?

Consider an $M\times M$ square lattice of photonic cavities. Each cavity
undergoes a two-photon drive process and a two-photon loss process:

$\displaystyle H_{x,y}=\lambda(a_{x,y}^{2}+(a^{\dagger}_{x,y})^{2}),\qquad
L_{2,x,y}=\sqrt{\kappa_{2}}a_{x,y}^{2}\,,$ (8)

where $a_{x,y}$ is the annihilation operator on site $(x,y)$. Next, let us
consider a parity-parity interaction between neighboring cavities:
$H_{S}=-\sum_{\langle ij\rangle}Q_{i}Q_{j}$, where $Q_{j}$ is the photon
parity operator at site $j$. Similar to the Ising model, at low temperatures,
such interaction will tend to align the parities of neighboring cavities via
the following local dissipators (for a microscopic derivation of the
dissipators, see the SM [52]):

$\displaystyle L_{x,y}^{(4)}$
$\displaystyle=\sqrt{\kappa_{nn}}a_{x,y}P_{x,y;\rightarrow}^{-}P_{x,y;\uparrow}^{-}P_{x-1,y;\rightarrow}^{-}P_{x,y-1;\uparrow}^{-}\,,$
$\displaystyle L_{x,y}^{(3)}$
$\displaystyle=\sqrt{\tilde{\kappa}_{nn}}a_{x,y}P_{x,y;\rightarrow}^{+}P_{x,y;\uparrow}^{-}P_{x-1,y;\rightarrow}^{-}P_{x,y-1;\uparrow}^{-}\,,$
(9)

where $a_{x,y}$ is the annihilation operator for the cavity at site $x,y$,
$\tilde{\kappa}_{nn}=\sqrt{\kappa_{1}\kappa_{nn}+\kappa_{1}^{2}}-\kappa_{1}$,
$\kappa_{1}$ is the single-photon loss rate (corresponding to the dissipator:
$L_{1,x,y}=\sqrt{\kappa_{1}}a_{x,y})$, $P^{\pm}_{x,y;\rightarrow}=(1\pm
Q_{x,y}Q_{x+1,y})/2,P^{\pm}_{x,y;\uparrow}=(1\pm Q_{x,y}Q_{x,y+1})/2,$ and
$Q_{x,y}=e^{i\pi a_{x,y}^{\dagger}a_{x,y}}$. The following states are the
steady states of the model in the absence of errors ($\kappa_{1}=0$) and span
the codespace:

$\displaystyle|\psi\rangle=c_{0}|\alpha_{e}\rangle|\alpha_{e}\rangle|\alpha_{e}\rangle\ldots+c_{1}|\alpha_{o}\rangle|\alpha_{o}\rangle|\alpha_{o}\rangle\ldots\,,$
(10)

for $|c_{0}|^{2}+|c_{1}|^{2}=1$.

For thermal systems, the existence of a passively-correcting quantum memory is
related to the presence of an extensive energy barrier which local errors must
overcome in order to create a logical bit-flip or phase-flip operation [53].
In the model described above, a logical bit-flip operation can be created via
local single-photon loss $L_{1,x,y}=\sqrt{\kappa_{1}}a_{x,y}$ only by passing
through a configuration with an extensive number of domain walls, which is
exponentially unlikely in the limit of large lattice size
$M\rightarrow\infty$. Similarly, a phase-flip error can only be generated by
taking the state $|\alpha_{e}\rangle\pm|\alpha_{o}\rangle$ to
$|\alpha_{e}\rangle\mp|\alpha_{o}\rangle$ for any of the cavities. However,
such a process is also unlikely to occur via dephasing perturbations
$L_{d,x,y}=\sqrt{\kappa_{d}}a^{\dagger}_{x,y}a_{x,y}$ which perturb states
locally in phase space, since the states
$|\pm\alpha\rangle\approx|\alpha_{e}\rangle\pm|\alpha_{o}\rangle$ are well
separated in phase space and an unstable fixed point sits between them [54].
The logical phase-flip errors are again exponentially unlikely as
$N\rightarrow\infty$.

The single-photon loss and the dephasing lead to terms proportional to
$a^{{\dagger}}a$ and $(a^{{\dagger}}a)^{2}$ in the Lindbladian, which result
in leakage out of the effective two-level codespace for each cavity into other
states of the cavity. This leakage poses a challenge for numerical simulation
since (unlike the Ising model) we need to keep track of more than two degrees
of freedom per lattice site. Nevertheless, we shall provide evidence for a
stable quantum memory by employing a variety of approximations.

First, let us consider an approximation that allows us to map the dynamics of
the photonic-Ising model directly to the classical-Ising model studied above.
Specifically, we introduce an idealized model by replacing the single-photon
loss dissipator $L_{1}=\sqrt{\kappa_{1}}a$ with $E_{1}=\sqrt{\kappa_{1}}b$,
where $b=aV$ and $V$ is the projector onto the codespace:
$V=|\alpha_{e}\rangle\langle\alpha_{e}|+|\alpha_{o}\rangle\langle\alpha_{o}|$.
We also assume an absence of dephasing errors, i.e. $\kappa_{d}=0$. This
allows us to treat each site as an effective two-level system
$|0\rangle=|\alpha_{e}\rangle,|1\rangle=|\alpha_{o}\rangle$, avoiding any
leakage out of the codespace. We similarly replace $a\rightarrow b$ in the
nearest-neighbor coupling dissipators (9) (except in the definition of $Q$).
The operator $b$ can be regarded as an “idealized bit flip” since, for $N\gg
1$, it takes the form
$b\approx\alpha(|\alpha_{e}\rangle\langle\alpha_{o}|+|\alpha_{o}\rangle\langle\alpha_{e}|)$.
The idealized model maps exactly to the Ising model studied above, with an
effective bit-flip error rate of $N\kappa_{1}$, an effective Ising-correction
rate of $N\kappa_{nn}$, and an inverse temperature
$\beta=\ln\left[(\kappa_{nn}+\kappa_{1})/\kappa_{1}\right]/8.$ We therefore
find that this model passively corrects against bit flips in the limit
$M\rightarrow\infty$ of the low-temperature phase. In the limit of large
driving strength and small single-photon loss, we expect the photonic-Ising
model to be well approximated by the idealized model since the state rarely
leaves the codespace. We provide quantitative evidence for this in the
Supplemental Material (SM) [52].

Dephasing, single-photon loss, and bit-flip recovery jumps ($L^{(3)}_{x,y}$
and $L^{(4)}_{x,y}$) cause leakage out of the codespace which is neglected
within the idealized model. It is natural to ask whether this leakage is
detrimental to the passively protected memory when the idealized model is no
longer a good approximation. We provide evidence that this is not the case by
studying a toy model which resembles the 2D model. Consider a single cavity
coupled to a spin-$1/2$ particle (described by Pauli operators $X,Y,Z$),
leading to two logical states $\ket{\downarrow}\ket{\alpha_{e}}$ and
$\ket{\downarrow}\ket{\alpha_{o}}$. The Hamiltonian and jump operators read
$h=\lambda(a^{2}+(a^{\dagger})^{2}),\ l_{2}=\sqrt{\kappa_{2}}a^{2},\
l_{1}=\sqrt{\kappa_{1}}Xa,\ l_{d}=\sqrt{\kappa_{d}}a^{{\dagger}}a,\
l_{nn}=\sqrt{\kappa_{nn}}\frac{1}{2}X(1-Z)a$. The model assumes that single-
photon loss is accompanied by a spin flip, while two-photon drive and
dephasing are not. The flip-recovery jump $l_{nn}$ is triggered by a flipped
spin state $\ket{\uparrow}$, similar to the bit-flip recovery jump caused by a
parity misalignment in 2D. Importantly, leakage caused by the noise processes
$l_{1},l_{d}$, and the flip-recovery jump is captured by this model. In the SM
[52], we analyze this model numerically and analytically. We find that the
initial state can always be perfectly restored via a decoder (up to
corrections exponentially small in $N$).

Figure 3: The mean-field phase diagram for $\kappa_{d}=\kappa_{1}$. The top
right corner shades the region where both $\langle Q\rangle$ and $\langle
a^{2}\rangle$ are non-zero. Both phase and bit-flip errors are protected. When
$\langle a^{2}\rangle\neq 0$ but $\langle Q\rangle=0$, we expect protection
only for phase errors. When $\langle a^{2}\rangle=0$, we expect the memory to
become fragile under either noise.

Finally, the stability of the memory can also be understood as the coexistence
of two order parameters: $\langle Q\rangle=\langle e^{i\pi
a^{{\dagger}}a}\rangle\neq 0$ indicates the ferromagnetic phase and therefore
suppression of bit-flip errors, while $\langle a^{2}\rangle\neq 0$ indicates
that the cat states are stabilized, implying suppression of phase-flip errors.
We use a product-state mean-field ansatz
$\rho=\bigotimes^{M}_{x,y=1}\rho_{x,y}$, where each $\rho_{x,y}$ is a density
matrix for a two-level system in the basis of $\ket{\pm\alpha_{MF}}$ for some
mean-field coherent parameter $\alpha_{MF}$. A non-trivial dissipative phase
of the system is identified by non-zero fixed points of $\langle Q\rangle$ and
$\langle a^{2}\rangle$. The mean-field solutions suggest that, for small
$\kappa_{1},\kappa_{d}$, the memory is protected against both phase and bit-
flip errors. When $\kappa_{1}$ or $\kappa_{d}$ exceeds a threshold, the order
parameters undergo two second-order phase transitions and the quantum memory
is no longer stable (see the SM [52]). The mean-field phase diagram is
sketched in Fig. 3.

Implementing the photonic-Ising dissipators. The key ingredients for our
proposal are the microscopic dissipators defined in Eq. (9). A direct approach
to achieve such terms involves engineering an Ising-like interaction between
cavity modes: $H_{S}\propto-\sum_{\langle ij\rangle}Q_{i}Q_{j}$. The natural
system-bath interaction of the form $\sum_{i}(a_{i}+a^{{\dagger}}_{i})\otimes
B_{i}$ (where $B_{i}$ acts on bath degrees of freedom) would then give rise to
the model described above (within the standard Born-Markov approximation)
[52]. The parity-parity interaction $H_{S}$ can be engineered from coupling
between high-impedance cavity modes and Josephson junctions [55, 56], as we
review in the SM [52].

Inspired by the microscopic dissipators Eq. (9), an alternative approach to
protect the memory involves digitally implementing a stochastic local error
decoder. In the SM [52], we provide an explicit description of how to achieve
such a local decoder autonomously without the need of measurements; however,
it requires that the local error decoding should be carried out at a rate that
scales linearly in the photon number $N$ of each cavity. The implementation
can be achieved simultaneously with the dephasing protection Eq. (8), making
it fault-tolerant [54]. Note that this is different from the active repetition
cat code in 1D [54] as we avoid the processing of non-local syndrome
information.

Discussion and outlook.—We proposed a photonic-Ising model that hosts robust
quantum memory under both single-photon loss and dephasing noise. We can
estimate the logical error rates in the photonic-Ising model as follows. While
the bit-flip error rate becomes extensive ($\sim O(N)$) in the limit of large
cavity photon number, the Ising-type interaction gives rise to an
exponentially-suppressed error rate $O(e^{-\gamma M})$ with $\gamma>0$ [57,
58, 59], resulting in a logical bit-flip error rate of $O(Ne^{-\gamma M})$.
Similarly, a single cavity yields a phase-flip error rate of
$O(e^{-\gamma^{\prime}N})$ with $\gamma^{\prime}>0$, while this is made
extensive by the spatially-extended lattice configuration, resulting in a
logical phase flip error rate of $O(M^{2}e^{-\gamma^{\prime}N})$. Harmonic
oscillators with small non-linearities and outstanding coherence
properties—and thus with large achievable $N$—can be found in a variety of
photonic and phononic systems (e.g. [38, 60]).

The realization of the parity-parity coupling based on Josephson junctions and
a high-impedance cavity mode is experimentally challenging. Future efforts
should consider other experimental schemes that can lead to the same effective
model. The photonic-Ising model can be generalized to adapt the Toom’s rule
[42], or to higher dimensions [61] for a more robust perturbative stability.
The full perturbative stability of the model remains an interesting open
question.

###### Acknowledgements.

Acknowledgment.—We sincerely thank Oles Shtanko, Victor Albert, and Daniel
Slichter for useful discussions. We thank Fernando Brandão for pointing out
Ref. [61] regarding the perturbative stability of a three-dimensional model.
S.L. was supported by the NIST NRC Research Postdoctoral Associateship. Y.-J.L
was supported by the Max Planck Gesellschaft (MPG) through the International
Max Planck Research School for Quantum Science and Technology (IMPRS-QST).
A.V.G. acknowledges funding by NSF QLCI (award No. OMA-2120757), DoE QSA, DoE
ASCR Accelerated Research in Quantum Computing program (award No. DE-
SC0020312), ARO MURI, AFOSR, DARPA SAVaNT ADVENT, DoE ASCR Quantum Testbed
Pathfinder program (award No. DE-SC0019040), U.S. Department of Energy Award
No. DE-SC0019449, NSF PFCQC program, and AFOSR MURI.

## References

  * Nielsen and Chuang [2002] M. A. Nielsen and I. Chuang,  _Quantum computation and quantum information_ (American Association of Physics Teachers, 2002).
  * Lidar and Brun [2013] D. A. Lidar and T. A. Brun,  _Quantum error correction_ (Cambridge university press, 2013).
  * Gottesman [1997] D. Gottesman,  _Stabilizer codes and quantum error correction_ , Ph.D. thesis, California Institute of Technology (1997).
  * Bravyi and Kitaev [1998] S. B. Bravyi and A. Y. Kitaev, arXiv preprint quant-ph/9811052 (1998).
  * Fowler  _et al._ [2012] A. G. Fowler, M. Mariantoni, J. M. Martinis, and A. N. Cleland, Phys. Rev. A 86, 032324 (2012).
  * Dennis  _et al._ [2002] E. Dennis, A. Kitaev, A. Landahl, and J. Preskill, J. Math. Phys. 43, 4452 (2002).
  * Gottesman  _et al._ [2001] D. Gottesman, A. Kitaev, and J. Preskill, Phys. Rev. A 64, 012310 (2001).
  * Bacon [2006] D. Bacon, Phys. Rev. A 73, 012340 (2006).
  * Yoshida [2011] B. Yoshida, Ann. Phys. 326, 2566 (2011).
  * Roberts and Bartlett [2020] S. Roberts and S. D. Bartlett, Phys. Rev. X 10, 031041 (2020).
  * Terhal [2015] B. M. Terhal, Rev. Mod. Phys. 87, 307 (2015).
  * Brown  _et al._ [2016] B. J. Brown, D. Loss, J. K. Pachos, C. N. Self, and J. R. Wootton, Rev. Mod. Phys. 88, 045005 (2016).
  * Bombin  _et al._ [2013] H. Bombin, R. W. Chhajlany, M. Horodecki, and M. A. Martin-Delgado, New J. Phys. 15, 055023 (2013).
  * Haah [2011] J. Haah, Phys. Rev. A 83, 042330 (2011).
  * Bravyi and Haah [2013] S. Bravyi and J. Haah, Phys. Rev. Lett. 111, 200501 (2013).
  * Chesi  _et al._ [2010] S. Chesi, B. Röthlisberger, and D. Loss, Phys. Rev. A 82, 022305 (2010).
  * Breuckmann  _et al._ [2016] N. P. Breuckmann, K. Duivenvoorden, D. Michels, and B. M. Terhal, arXiv preprint arXiv:1609.00510 (2016).
  * Alicki  _et al._ [2009] R. Alicki, M. Fannes, and M. Horodecki, J. Phys. A: Math. Theor. 42, 065303 (2009).
  * Alicki  _et al._ [2010] R. Alicki, M. Horodecki, P. Horodecki, and R. Horodecki, Open Sys. Inf. Dyn. 17, 1 (2010).
  * Paz and Zurek [1998] J. P. Paz and W. H. Zurek, Proc. R. Soc. A: Math. Phys. Eng. Sci. 454, 355 (1998).
  * Barnes and Warren [2000] J. P. Barnes and W. S. Warren, Phys. Rev. Lett. 85, 856 (2000).
  * Ahn  _et al._ [2002] C. Ahn, A. C. Doherty, and A. J. Landahl, Phys. Rev. A 65, 042301 (2002).
  * Sarovar and Milburn [2005] M. Sarovar and G. J. Milburn, Phys. Rev. A 72, 012306 (2005).
  * Oreshkov and Brun [2007] O. Oreshkov and T. A. Brun, Phys. Rev. A 76, 022318 (2007).
  * Kerckhoff  _et al._ [2010] J. Kerckhoff, H. I. Nurdin, D. S. Pavlichin, and H. Mabuchi, Phys. Rev. Lett. 105, 040502 (2010).
  * Lihm  _et al._ [2018] J.-M. Lihm, K. Noh, and U. R. Fischer, Phys. Rev. A 98, 012317 (2018).
  * Pastawski  _et al._ [2011] F. Pastawski, L. Clemente, and J. I. Cirac, Phys. Rev. A 83, 012304 (2011).
  * Mirrahimi  _et al._ [2014] M. Mirrahimi, Z. Leghtas, V. V. Albert, S. Touzard, R. J. Schoelkopf, L. Jiang, and M. H. Devoret, New J. Phys. 16, 045014 (2014).
  * Kapit [2016] E. Kapit, Phys. Rev. Lett. 116, 150501 (2016).
  * Reiter  _et al._ [2017] F. Reiter, A. S. Sørensen, P. Zoller, and C. A. Muschik, Nat. Commun. 8, 1822 (2017).
  * Fujii  _et al._ [2014] K. Fujii, M. Negoro, N. Imoto, and M. Kitagawa, Phys. Rev. X 4, 041039 (2014).
  * Ofek  _et al._ [2016] N. Ofek, A. Petrenko, R. Heeres, P. Reinhold, Z. Leghtas, B. Vlastakis, Y. Liu, L. Frunzio, S. M. Girvin, L. Jiang, M. Mirrahimi, M. H. Devoret, and R. J. Schoelkopf, Nature 536, 441 (2016).
  * Gyenis  _et al._ [2021] A. Gyenis, P. S. Mundada, A. Di Paolo, T. M. Hazard, X. You, D. I. Schuster, J. Koch, A. Blais, and A. A. Houck, PRX Quantum 2, 010339 (2021).
  * Campagne-Ibarcq  _et al._ [2020] P. Campagne-Ibarcq, A. Eickbusch, S. Touzard, E. Zalys-Geller, N. E. Frattini, V. V. Sivak, P. Reinhold, S. Puri, S. Shankar, R. J. Schoelkopf, L. Frunzio, M. Mirrahimi, and M. H. Devoret, Nature 584, 368 (2020).
  * Leghtas  _et al._ [2015] Z. Leghtas, S. Touzard, I. M. Pop, A. Kou, B. Vlastakis, A. Petrenko, K. M. Sliwa, A. Narla, S. Shankar, M. J. Hatridge, _et al._ , Science 347, 853 (2015).
  * Heeres  _et al._ [2017] R. W. Heeres, P. Reinhold, N. Ofek, L. Frunzio, L. Jiang, M. H. Devoret, and R. J. Schoelkopf, Nat. Commun. 8, 94 (2017).
  * Lescanne  _et al._ [2020] R. Lescanne, M. Villiers, T. Peronnin, A. Sarlette, M. Delbecq, B. Huard, T. Kontos, M. Mirrahimi, and Z. Leghtas, Nat. Phys. 16, 509 (2020).
  * Berdou  _et al._ [2022] C. Berdou, A. Murani, U. Reglade, W. Smith, M. Villiers, J. Palomo, M. Rosticher, A. Denis, P. Morfin, M. Delbecq, _et al._ , arXiv preprint arXiv:2204.09128 (2022).
  * Lieu  _et al._ [2020] S. Lieu, R. Belyansky, J. T. Young, R. Lundgren, V. V. Albert, and A. V. Gorshkov, Phys. Rev. Lett. 125, 240405 (2020).
  * Gertler  _et al._ [2021] J. M. Gertler, B. Baker, J. Li, S. Shirol, J. Koch, and C. Wang, Nature 590, 243 (2021).
  * Leghtas  _et al._ [2013] Z. Leghtas, G. Kirchmair, B. Vlastakis, R. J. Schoelkopf, M. H. Devoret, and M. Mirrahimi, Phys. Rev. Lett. 111, 120501 (2013).
  * Toom [1974] A. L. Toom, Problems Inform. Transmission 10, 239 (1974).
  * Lindblad [1976] G. Lindblad, Commun. Math. Phys. 48, 119 (1976).
  * fn [1] If the code space spans the kernel of $\mathcal{L}_{r}$, then a natural choice for the decoder is $\mathcal{E}_{r}=\lim_{t\to\infty}e^{\mathcal{L}_{r}t}$. Such a decoder is typically used for the cat code. However, if $\mathcal{L}_{r}$ has other steady states which are not in the code space, such a decoder is not ideal. We do not use this decoder for the Ising model, since its $\mathcal{L}_{r}$ has steady states outside of the code space.
  * Chamberland  _et al._ [2022] C. Chamberland, K. Noh, P. Arrangoiz-Arriola, E. T. Campbell, C. T. Hann, J. Iverson, H. Putterman, T. C. Bohdanowicz, S. T. Flammia, A. Keller, G. Refael, J. Preskill, L. Jiang, A. H. Safavi-Naeini, O. Painter, and F. G. Brandão, PRX Quantum 3, 010329 (2022).
  * Gilles  _et al._ [1994] L. Gilles, B. M. Garraway, and P. L. Knight, Phys. Rev. A 49, 2785 (1994).
  * fn [3] We could also include jumps that flip a spin if there are two domain walls surrounding it, e.g. $L^{(2)}=\sqrt{\gamma}X_{x,y}P^{+}_{x,y;\rightarrow}P^{+}_{x,y;\uparrow}P^{-}_{x-1,y;\rightarrow}P^{-}_{x,y-1;\uparrow}$, since such a process does not change the energy and hence respects detailed balance (for any rate $\gamma$). Consequently, such jumps do not change the exact thermal steady state solutions. We choose not to include such jumps since the mean-field solution is more accurate without them.
  * Plenio and Knight [1998] M. B. Plenio and P. L. Knight, Rev. Mod. Phys. 70, 101 (1998).
  * Daley [2014] A. J. Daley, Adv. Phys. 63, 77 (2014).
  * Buča and Prosen [2012] B. Buča and T. Prosen, New J. Phys. 14, 073007 (2012).
  * fn [2] If the symmetry-broken initial state $|\bar{0}\rangle$ is recoverable, then so is any state in the code space. This is because an arbitrary state in the code space can be expressed as: $c_{0}|\bar{0}\rangle+c_{1}|\bar{1}\rangle=(c_{0}-c_{1})|\bar{0}\rangle+c_{1}(|\bar{0}\rangle+|\bar{1}\rangle)$. Since the dynamics is linear, and the symmetric state $|\bar{0}\rangle+|\bar{1}\rangle$ can always be properly recovered in the absence of $Z$ dephasing, this implies that any state in the code space is recoverable if $|\bar{0}\rangle$ is recoverable.
  * [52] See the Supplemental Material.
  * Bravyi and Terhal [2009] S. Bravyi and B. Terhal, New J. Phys. 11, 043029 (2009).
  * Guillaud and Mirrahimi [2019] J. Guillaud and M. Mirrahimi, Phys. Rev. X 9, 041053 (2019).
  * Cohen [2017] J. Cohen,  _Autonomous quantum error correction with superconducting qubits_ , Ph.D. thesis, Ecole Normale Superieure (2017).
  * Cohen  _et al._ [2017] J. Cohen, W. C. Smith, M. H. Devoret, and M. Mirrahimi, Phys. Rev. Lett. 119, 060503 (2017).
  * Thomas [1989] L. E. Thomas, Commun. Math. Phys. 126, 1 (1989).
  * Schonmann [1987] R. H. Schonmann, Commun. Math. Phys. 112, 409 (1987).
  * Randall [2006] D. Randall, in  _Proceedings of the Seventeenth Annual ACM-SIAM Symposium on Discrete Algorithm_ , SODA ’06 (Society for Industrial and Applied Mathematics, USA, 2006) p. 870–879.
  * Home  _et al._ [2011] J. P. Home, D. Hanneke, J. D. Jost, D. Leibfried, and D. J. Wineland, New J. Phys. 13, 073026 (2011).
  * Poulin  _et al._ [2019] D. Poulin, R. G. Melko, and M. B. Hastings, Phys. Rev. B 99, 094103 (2019).
  * Breuer and Petruccione [2002] H. Breuer and F. Petruccione,  _The Theory of Open Quantum Systems_ (Oxford University Press, 2002).
  * Weimer  _et al._ [2010] H. Weimer, M. Müller, I. Lesanovsky, P. Zoller, and H. P. Büchler, Nat. Phys. 6, 382 (2010).
  * Barreiro  _et al._ [2011] J. T. Barreiro, M. Müller, P. Schindler, D. Nigg, T. Monz, M. Chwalla, M. Hennrich, C. F. Roos, P. Zoller, and R. Blatt, Nature 470, 486 (2011).

Supplemental Material for “Candidate for a passively protected quantum memory
in two dimensions”

The Supplemental Material is organized as follows: In Sec. 1, we provide
numerical evidence that the idealized bit flip approximation introduced in the
main text is reasonable in the limit of large drive, small single-photon loss,
and no dephasing. In Sec. 2, we study a “toy model”, which was introduced in
the main text, which mimics the dynamics of the 2D photonic-Ising model, and
which is tractable both numerically and analytically. This model suggests that
leakage out of the codespace arising from single-photon loss and dephasing is
not detrimental to passive correction. In Sec. 3, we provide details on the
mean-field theory order parameters described in the main text. In Sec. 4 we
show that the model studied in the main text can be achieved for an Ising-like
Hamiltonian interaction between cavity modes in the presence of the natural
system-bath coupling. In Sec. 5 we describe a way to achieve an Ising-like
interaction between cavity modes in a superconducting circuit scheme. In Sec.
6, we provide an alternative way of achieving the model in the main text;
instead of engineering a Hamiltonian interaction, we describe a way to
engineer only the desired dissipators. This scheme relies on the presence of
fast unitary gates and ancilla resets.

##  I 1\. Idealized bit flip approximation

In this section, we elaborate on the idealized bit flip approximation used in
the main text. In experiments, the bit flip error for a single photonic cat
qubit is generated via single-photon loss $L_{1}=\sqrt{\kappa_{1}}a$. However,
in order to map our many-body-cat-qubit system to the 2D Ising model, we must
replace this noise generator with an “idealized bit flip”, represented via the
jump operator: $E_{1}=\sqrt{\kappa_{1}}aV$ where $V$ is a projector onto the
codespace. We provide evidence that $E_{1}$ is a reasonable approximation for
$L_{1}$ in the limit of small single-photon loss and large two-photon drive
(compared to the two-photon loss rate), which is the relevant regime for
modern experiments involving photonic cat qubits [45]. We also assume the
absence of photon dephasing. To this end, we shall present two models for a
single cavity and show that their steady states and dissipative gaps converge
in this limit.

Model 1 has the standard single-photon loss term which is expected to appear
in experiment. Model 2 has the “idealized bit flip” which is needed to make
numerical progress.

Model 1: Let us consider a single photonic cavity in the presence of two-
photon drive $H=\lambda[a^{2}+(a^{\dagger})^{2}]$, two photon loss
$L_{2}=\sqrt{\kappa_{2}}a^{2}$, and single-photon loss
$L_{1}=\sqrt{\kappa_{1}}a$. It is convenient to utilize the gauge freedom of
the Lindbladian to eliminate the Hamiltonian by incorporating it in a
dissipative term. The following two dissipators share the same master equation
as the model just described:

$\displaystyle L_{c}$
$\displaystyle=\sqrt{\kappa_{2}}(a^{2}-\alpha^{2}),\qquad\alpha=\sqrt{\frac{\lambda}{\kappa_{2}}}e^{-i\pi/4}$
(S1) $\displaystyle L_{1}$ $\displaystyle=\sqrt{\kappa_{1}}a.$ (S2)

The dissipator $L_{c}$ will cause states in the Hilbert space to evolve
towards the coherent states $|\pm\alpha\rangle$, which are dark states of
$L_{c}$. We thus find that $L_{c}$ generates the “recovery” part of the
Lindbladian, while $L_{1}$ generates bit flip errors and causes leakage out of
the codespace.

From the perspective of quantum trajectories, single-photon loss causes the
amplitude of a coherent state to decay due to the non-Hermitian Hamiltonian
term proportional to $\kappa_{1}a^{\dagger}a$ which (by itself) causes the
coherent state parameter to decay via $\alpha e^{-\kappa_{1}t}$. The two-
photon drive process ensures that the steady state amplitude remains non-zero,
but nevertheless the photon population decreases due to the single-photon
loss. Within mean-field theory, the average number $\bar{n}$ of photons in the
cavity satisfies

$\bar{n}=\frac{2\lambda-\kappa_{1}}{2\kappa_{2}}.$ (S3)

This suggests that, in the limit of $\lambda/\kappa_{2}\gg 1$, the steady
state of the system should start to converge to a coherent state
$|\pm\mu\rangle$ with a shifted amplitude:

$a|\pm\mu\rangle=\pm\mu|\pm\mu\rangle,\qquad\mu=\sqrt{\frac{2\lambda-\kappa_{1}}{2\kappa_{2}}}e^{-i\pi/4}.$
(S4)

Numerics suggest that the true steady state of the system will be a mixture of
several pure states [39]. However, the steady state will have large overlap
with the states $|\pm\mu\rangle$. In the limit $\kappa_{1}/\kappa_{2}\ll 1$,
the steady state will start to converge to a mixture of the states
$|\pm\mu\rangle$.

We can confirm this via numerical simulations. In Fig. S1 we plot the overlap
of the steady state with $|\mu\rangle$ as a function of the drive strength
$\lambda/\kappa_{2}$, for different choices of $\kappa_{1}/\kappa_{2}$. We
find that the steady state of the system approaches $|\mu\rangle$ in the limit
$\lambda/\kappa_{2}\gg 1,\kappa_{1}/\kappa_{2}\ll 1$. These parameters are in
a regime that is relevant for modern experiments [45]. We also plot the
dissipative gap, which scales linearly with the drive strength.

Figure S1: Model 1: (a) Expectation value of $|\mu\rangle\langle\mu|$ in the
steady state of the model described in Eqs. (S1), (S2) with
$\lambda/\kappa_{2}=N$ for different choices of $\kappa_{1}/\kappa_{2}$. In
the limit $\lambda/\kappa_{2}\gg 1,\kappa_{1}/\kappa_{2}\ll 1$, the system
converges to the coherent state $|\mu\rangle$. We use exact Lindblad evolution
starting from the initial state $|\alpha\rangle$ and evolving for a time
$t=200/\kappa_{2}$ to reach the steady state. (b) The dissipative gap
$\Lambda_{g}$ scales linearly as a function of the drive strength, for
$\kappa_{1}/\kappa_{2}=10^{-3}$. Figure S2: Model 2: (a) Expectation value of
$|\mu\rangle\langle\mu|$ in the steady state of the model described in Eqs.
(S6), (S7) with parameters: $\lambda/\kappa_{2}=N$ for different choices of
$\kappa_{1}/\kappa_{2}$. In the limit, $\lambda/\kappa_{2}\gg
1,\kappa_{1}/\kappa_{2}\ll 1$ the system converges to the coherent state
$|\mu\rangle$. We use exact Lindblad evolution starting from the initial state
$|\alpha\rangle$ and evolving for a time $t=200/\kappa_{2}$ to reach the
steady state. (b) The dissipative gap $\Lambda_{g}$ scales linearly as a
function of the drive strength, for $\kappa_{1}/\kappa_{2}=10^{-3}$.

Beyond a shift in the coherent state amplitude, single-photon loss also has
the effect of reducing the qubit-steady-state structure to a classical-bit-
steady-state structure. Only classical mixtures of coherent states are stable,
while off-diagonal coherences have a finite lifetime:

$\rho_{ss}\approx c|\mu\rangle\langle\mu|+(1-c)|-\mu\rangle\langle-\mu|.$ (S5)

for $c\in[0,1],\lambda/\kappa_{2}\gg 1,\kappa_{1}/\kappa_{2}\ll 1$. The steady
state is thus two dimensional, enough only to store a classical bit.

Model 2: Let us now consider a different model which will have the same steady
state and dissipative gap in the limit $\lambda/\kappa_{2}\gg
1,\kappa_{1}/\kappa_{2}\ll 1$, but will involve the “idealized bit flip”
rather than single-photon loss. Consider the dissipators

$\displaystyle L_{c}$
$\displaystyle=\sqrt{\kappa_{2}}(a^{2}-\alpha^{2}),\qquad\alpha=\sqrt{\frac{\lambda}{\kappa_{2}}}e^{-i\pi/4}$
(S6) $\displaystyle E_{1}$
$\displaystyle=\sqrt{\kappa_{1}}b=\sqrt{\kappa_{1}}aV,\qquad
V=|\alpha_{e}\rangle\langle\alpha_{e}|+|\alpha_{o}\rangle\langle\alpha_{o}|$
(S7)

where
$|\alpha_{e}\rangle\sim|\alpha\rangle+|-\alpha\rangle,|\alpha_{o}\rangle\sim|\alpha\rangle-|-\alpha\rangle$.
In this model, the dissipator $E_{1}$ does not cause any leakage of photons
out of $|\alpha\rangle$. This is because the non-Hermitian Hamiltonian term
proportional to $E_{1}^{\dagger}E_{1}$ keeps superpositions of
$|\pm\alpha\rangle$ in this subspace (due to the projector $V$). Nevertheless,
the term $E_{1}$ ensures that quantum superpostions of $|\pm\alpha\rangle$ are
unstable, while classical mixtures are stable. The steady state starts to
converge to the following state in the limit of large drive
$\lambda/\kappa_{2}\gg 1$:

$\rho_{ss}\approx
c|\alpha\rangle\langle\alpha|+(1-c)|-\alpha\rangle\langle-\alpha|,$ (S8)

for $c\in[0,1]$.

The overlap between $|\alpha\rangle$ and $|\mu\rangle$ satisfies

$|\langle\alpha|\mu\rangle|^{2}=\exp[-\frac{\kappa_{1}^{2}}{16\kappa_{2}\lambda}]\approx
1-\frac{\kappa_{1}^{2}}{16\kappa_{2}\lambda}+\ldots$ (S9)

This implies that the deviation from unity scales as $\kappa_{1}^{2}$ when
$\kappa_{2}\lambda\gg\kappa_{1}^{2}$. We confirm this in Fig. S2: The
deviation between the steady state of Model 2 and $|\mu\rangle$ scales
quadratically with $\kappa_{1}$ in the limit of large drive. We also plot the
dissipative gap, which again scales linearly with the drive strength.

We have shown that Models 1 and 2 converge to each other in terms of their
steady state and their dissipative gap in the limit $\lambda/\kappa_{2}\gg
1,\kappa_{1}/\kappa_{2}\ll 1$. This suggests that Model 2 is a reasonable
approximation for Model 1 in this regime. Intuitively, this happens because
the system quickly evolves toward the codespace, such that the projector term
$V$ acts trivially on the state. In the main text, we demonstrated that Model
2 passively corrects against bit flip errors via the Ising-like dissipators
described above. We expect Model 1 to behave in qualitatively the same manner
after the replacement of $b\rightarrow a$.

We note that, although we used the limit $\lambda/\kappa_{2}\gg
1,\kappa_{1}/\kappa_{2}\ll 1$ to establish the exact mapping to the Ising
model, we do not expect that this limit is needed to preserve quantum
information in general. Rather, the system only needs to stay within the
ordered phase (see Fig. 3 in the main text and SM Sec. 3). A relatively small
$\kappa_{1}$ ensures that the steady state of the dynamics is a mixed state.
Nevertheless, we expect that this mixed state will be a “noiseless subsystem”,
which implies that it can be decoded with a channel superoperator at the end
of the dynamics.

##  II 2\. Toy model

The Ising-inspired bit-flip recovery jump operators [Eqs. (9) in the main
text] by themselves will not give rise to protection against single-photon
loss in the absence of a drive, since single-photon loss will cause the system
to evolve to a vacuum state. In this section, we argue that, when the bit-flip
recovery is coupled with the driving, the resulting environment is able to
protect against both dephasing and single-photon loss errors.

Ideally, we would like to numerically simulate the 2D array of $M^{2}$ cat
qubits introduced in the main text. However, such a simulation is
computationally expensive. We restrict ourselves to the toy model introduced
in the main text: a single cat qubit coupled to a two-level system, the latter
described by Pauli operators $X,Y,Z$. The logical states of this toy system
are defined as $\ket{\downarrow}\ket{\alpha_{o}}$ and
$\ket{\downarrow}\ket{\alpha_{e}}$, where $\ket{\alpha_{e}},\ket{\alpha_{o}}$
are the logical states for a single cat qubit. The noise and recovery jump
operators are modified to

$\displaystyle l_{c}$
$\displaystyle=\sqrt{\kappa_{2}}(a^{2}-\alpha^{2}),\qquad\alpha=\sqrt{\frac{\lambda}{\kappa_{2}}}e^{-i\pi/4}$
(S10) $\displaystyle l_{1}$ $\displaystyle=\sqrt{\kappa_{1}}Xa,\quad
l_{d}=\sqrt{\kappa_{d}}a^{{\dagger}}a,$ (S11) $\displaystyle l_{nn}$
$\displaystyle=\sqrt{\kappa_{nn}}\frac{1}{2}X(1-Z)a,$ (S12)

where $l_{c}$ generates a Lindbladian that is equivalent to the combined
action of $h$ and $l_{2}$ in the main text. In this toy model, the spin-$1/2$
particle is essentially a “classical bit” that takes the discrete value of up
or down. Any single-photon loss event is always accompanied by a flip of the
spin. A bit-flip recovery for the cat qubit can then be achieved by checking
the orientation of the spin: an annhilation operator $a$ is applied to the
cavity if the spin points upwards, otherwise nothing happens. This mimics the
full 2D case where a bit-flip recovery jump is triggered by a parity
misalignment between nearest-neighbor cat qubits. The difference between the
2D model and the toy model is that the latter always knows when an odd number
of single photon-loss events has occurred. What remains to be tested is
whether the errors can be corrected by introducing the bit-flip recovery jump.

Suppose we initialize the dynamics with a generic state in the codespace. We
consider the following two scenarios: We choose the model with (i)
$\kappa_{2}=1,\kappa_{d}=0.1,\kappa_{1}=0.1,\kappa_{nn}=0$ and (ii)
$\kappa_{2}=1,\kappa_{d}=0.1,\kappa_{1}=0.1,\kappa_{nn}=0.3$. The system size
parameter is $N=\lambda/\kappa_{2}$ with $N\to\infty$ representing the
thermodynamic limit. The initial state is first evolved with this Lindbladian
for duration $T=15$, then followed by the corresponding noiseless Lindbladian
evolution ($\kappa_{d}=\kappa_{1}=0$) for another $T=15$. In the end, we
compute the fidelity between the final state and the initial state. The
results for the two scenarios are shown in Fig. S3 for different $N$.

The results clearly show distinct behaviors. For case (i), where
$\kappa_{nn}=0$, the single-photon loss causes uncorrectable errors in the
stored memory, leading to a saturated fidelity of $1/2$ (due to an equal
mixture of the flipped and unflipped states) as $N$ increases. For case (ii),
where $\kappa_{nn}\neq 0$, increasing $N$ leads to a fidelity exponentially
close to the ideal value of 1.

As a sanity check, let us consider the same numerical simulation but modify
case (ii) by setting $\kappa_{nn}$ = 0 during the noiseless dynamics (while
still keeping $\kappa_{nn}=0.3$ during the noisy dynamics). The results of the
simulation are shown in Fig. S4. In case (i), which is identical to the one
studied in Fig. S3, the fidelity relaxes to $1/2$ regardless of the system
size as before. The modified case (ii) shows a saturated fidelity between
$1/2$ and $1$, suggesting a partial preservation of the initial quantum
memory. This again confirms the dynamical quantum memory protection arising
from the flip-recovery jump and two-photon drive.

Figure S3: We initialize the dynamics with state
$\rho_{i}=|\psi\rangle\langle\psi|$, where
$|\psi\rangle=\frac{1}{\sqrt{5}}\ket{\downarrow}\ket{\alpha_{e}}+\frac{2e^{i\pi/4}}{\sqrt{5}}\ket{\downarrow}\ket{\alpha_{o}}$.
(a) Overlap between the initial and final states for $\kappa_{nn}=0$ [case
(i)] and $\kappa_{nn}/\kappa_{2}=0.3$ [case (ii)] as $N=\lambda/\kappa_{2}$
increases. Parameters: $\kappa_{d}/\kappa_{2}=0.1,\kappa_{1}/\kappa_{2}=0.1$.
(b) For case (ii), i.e. $\kappa_{nn}\neq 0$, the log scale plot shows that the
fidelity converges exponentially quickly to 1 as $N\to\infty$. Figure S4:
Repeating the same simulation as in Fig. S3, except the recovery (noiseless)
Lindblad evolution is done using $\kappa_{nn}=\kappa_{d}=\kappa_{1}=0$ in both
cases (i) and (ii) [not just case (i)]. The overlap between the initial and
final states is larger than $1/2$ when $\kappa_{nn}\neq 0$ during the noisy
dynamics [modified case (ii)], while the overlap saturates to $1/2$ when
$\kappa_{nn}=0$ [the original case (i)]. In the thermodynamic limit, the
overlap values agree with mean-field results (shown as horizontal dashed
lines).

Mean-field analysis of the toy model.—We use a mean field approach to show
that, despite the spin-boson coupling in our toy model, the $\mathbb{Z}_{2}$
symmetry-breaking phase diagram of the single cat qubit is reproduced. Given
an observable $\hat{O}$ and a Lindbladian term $\mathcal{L}$ generated by the
jump operator $L$, the expectation value obeys

$\Tr[\hat{O}\mathcal{L}\rho]=-\frac{1}{2}\Tr\left[[\hat{O},L^{{\dagger}}]L\rho+L^{{\dagger}}[L,\hat{O}]\rho\right].$
(S13)

Using this, we can derive a coupled set of mean-field equations of motion for
$\langle a\rangle$ and $\langle Z\rangle$:

$\displaystyle\frac{d}{dt}\langle a\rangle$ $\displaystyle=-i\lambda\langle
a^{{\dagger}}\rangle-\frac{1}{2}\left(\kappa_{1}+\kappa_{d}+\frac{\kappa_{nn}}{2}(1-\langle
Z\rangle)\right)\langle a\rangle-\kappa_{2}|\alpha|^{2}\langle a\rangle,$
(S14) $\displaystyle\frac{d}{dt}\langle Z\rangle$
$\displaystyle=-2\kappa_{1}|\alpha|^{2}\langle
Z\rangle+\kappa_{nn}|\alpha|^{2}(1-\langle Z\rangle).$ (S15)

This yields the mean-field fixed point solutions for both observables

$\displaystyle\langle Z\rangle$
$\displaystyle=\frac{\kappa_{nn}}{\kappa_{nn}+2\kappa_{1}},$ (S16)
$\displaystyle\kappa_{2}|\alpha|^{2}$
$\displaystyle=|\lambda|-\frac{1}{2}\left(\kappa_{1}+\kappa_{d}+\frac{\kappa_{1}\kappa_{nn}}{\kappa_{nn}+2\kappa_{1}}\right).$
(S17)

The expression closely matches the simulation in the thermodynamic limit (see
Fig. S4).

It is interesting to note that if $\kappa_{1}/\kappa_{2}$ is small enough,
then any non-zero $\kappa_{nn}$ can give rise to a stable memory ($\langle
Z\rangle,\langle a\rangle\neq 0$). On the other hand, if
$\kappa_{1}/\kappa_{2}$ is large, a large $\kappa_{nn}$ can destabilize the
memory, leading to $\langle a\rangle=0$.

##  III 3\. Mean-field solution for the 2D photonic-Ising model

In this section, we present the mean-field solution for the 2D photonic-Ising
model. The mean-field analysis shows the existence of two symmetry-breaking
transitions via two order parameters: $a^{2}$ and $Q\equiv e^{i\pi
a^{{\dagger}}a}$.

We consider a product-state mean-field ansatz
$\rho=\bigotimes^{M}_{x,y=1}\rho_{x,y}$. At each site, $\rho_{x,y}$ is a
density matrix for a two-level system in the basis of $\ket{\pm\alpha_{MF}}$
for some coherent parameter $\alpha_{MF}$. We first begin by deriving the
mean-field equation for $Q=e^{i\pi a^{{\dagger}}a}$. Note that all the terms
that commute with $Q$ do not contribute to the time evolution. We are
therefore left to consider only the single-photon loss term and the bit-flip
correction term. Using Eq. (S13), we obtain

$\displaystyle\frac{d\langle Q\rangle}{dt}=-2\left(\kappa_{1}\langle
a^{{\dagger}}aQ\rangle+\kappa_{nn}\langle
a^{{\dagger}}aQP_{\kappa_{nn}}\rangle+\tilde{\kappa}_{nn}\langle
a^{{\dagger}}aQP_{\tilde{\kappa}_{nn}}\rangle\right),$ (S18)

where $P_{\kappa_{nn}},P_{\tilde{\kappa}_{nn}}$ are sums of projectors onto
different parity configurations with rates $\kappa_{nn},\tilde{\kappa}_{nn}$,
as introduced in the main text. Within mean-field theory, we replace the
expectations by a product of expectations at each site, yielding

$\displaystyle-\frac{1}{2|\alpha|^{2}}\frac{d\langle
Q\rangle}{dt}=\frac{\kappa_{nn}-4\tilde{\kappa}_{nn}}{16}\langle
Q\rangle^{5}+\frac{\kappa_{nn}+4\tilde{\kappa}_{nn}}{8}\langle
Q\rangle^{3}-\left(\frac{3\kappa_{nn}+4\tilde{\kappa}_{nn}}{16}-\kappa_{1}\right)\langle
Q\rangle.$ (S19)

Similarly, we can derive the mean-field equation for $a^{2}$:

$\frac{d\langle a^{2}\rangle}{dt}=-\kappa_{2}(2\langle
a^{{\dagger}}aa^{2}\rangle+\langle a^{2}\rangle)-i\lambda(2\langle
a^{{\dagger}}a\rangle+1)-\kappa_{1}\langle a^{2}\rangle-2\kappa_{d}\langle
a^{2}\rangle-\kappa_{nn}\langle
a^{2}P_{\kappa_{nn}}\rangle-\tilde{\kappa}_{nn}\langle
a^{2}P_{\tilde{\kappa}_{nn}}\rangle.$ (S20)

With the mean-field ansatz, we may approximate $\langle
a^{{\dagger}}aa^{2}\rangle\approx|\alpha_{MF}|^{2}\langle a^{2}\rangle$. We
also have $\langle a^{2}P_{\tilde{\kappa}_{nn}}\rangle=\langle
a^{2}\rangle\langle P_{\tilde{\kappa}_{nn}}\rangle$ and $\langle
a^{2}P_{\kappa_{nn}}\rangle=\langle a^{2}\rangle\langle
P_{\kappa_{nn}}\rangle$. After some algebra, the mean-field fixed points at
the thermodynamic limit (e.g. $\kappa_{2}\to 0$) can be found to satisfy

$\displaystyle\langle Q\rangle^{2}$
$\displaystyle=\frac{2\sqrt{\kappa_{nn}^{2}-4\kappa_{1}(\kappa_{nn}-4\tilde{\kappa}_{nn})}-\kappa_{nn}-4\tilde{\kappa}_{nn}}{\kappa_{nn}-4\tilde{\kappa}_{nn}},$
(S21) $\displaystyle|\alpha_{MF}|^{2}$
$\displaystyle=\frac{2\lambda-\kappa_{1}-2\kappa_{d}-\gamma_{4}\langle
Q\rangle^{4}-\gamma_{2}\langle Q\rangle^{2}-\gamma_{0}}{2\kappa_{2}},$ (S22)

where $\gamma_{4}=(-3\kappa_{nn}+4\tilde{\kappa}_{nn})/16$,
$\gamma_{2}=(\kappa_{nn}-4\tilde{\kappa}_{nn})/8$, and
$\gamma_{0}=(\kappa_{nn}+4\tilde{\kappa}_{nn})/16$. In addition, $\langle
Q\rangle^{2}\neq 0$ is only possible when $|\alpha_{MF}|^{2}\neq 0$.
Intuitively, when $\langle a^{2}\rangle=0$, the cavity will lose coherence and
decay to the vacuum due to the noise. The logical states are no longer well-
defined in this case.

It is important to note that the mean-field solution suggests that the leakage
caused by both finite $\kappa_{1}$ and finite $\kappa_{d}$ is compensated by
the two-photon drive. The effect of this leakage amounts to a shift in the
steady state coherent parameter.

##  IV 4\. A microscopic derivation of the photonic-Ising dissipators

Here we establish an explicit connection between the Hamiltonian approach for
the photonic-Ising model and the microscopic Lindbladian approach. This
section (Section 4) provides an example where the proposed photonic-Ising
dissipators emerge naturally from a microscopic coupling, unlike the example
in the previous section. In the next section (Section 5), we will discuss an
experimental protocol that realizes the desired Hamiltonian coupling based on
superconducting circuits.

###  IV.1 The microscopic generators

Let us start by considering a microscopic Hamiltonian of both the system and
the bath:

$H=H_{S}+H_{B}+H_{SB},$ (S23)

where $H_{S}=-J\sum_{\langle i,j\rangle}Q_{i}Q_{j}$, and $H_{B},H_{SB}$ are
the bath and the system-bath coupling Hamiltonian, respectively. Notice that,
in contrast to the main text, we explicitly introduced an energy scale $J>0$
in $H_{S}$ to help us carry out the analysis. We consider
$H_{SB}=\sum_{i}(a_{i}+a^{{\dagger}}_{i})\otimes B_{i}$, where $a_{i}$ is the
annihilation operator on the photonic-Ising lattice and $B_{i}$ is some
Hermitian local operator on the bath. We assume that the bath is large and the
interaction $H_{SB}$ is weak such that effects of the coupling on the bath is
fast and can be neglected, i.e. the full density matrix approximately
factorizes into a product of a system density matrix and a bath density
matrix: $\rho(t)\approx\rho_{S}(t)\otimes\rho_{B},\ \forall t$. Provided that
the standard Born-Markov approximation (i.e. the smallness of the influence of
the system-bath coupling on the bath) and the rotating-wave approximation are
valid [62], we can derive the Master equation in the interaction picture as

$\frac{d\rho_{S}}{dt}=-i[H^{\prime},\rho_{S}]+\sum_{\omega}\sum_{i,j}\gamma_{i,j}(\omega)\left(A_{i}(\omega)\rho_{S}A_{j}^{{\dagger}}(\omega)-\frac{1}{2}\\{A^{{\dagger}}_{j}(\omega)A_{i}(\omega),\rho_{S}\\}\right).$
(S24)

Here $H^{\prime}$ is the Lamb-shift Hamiltonian which we will define below.
The operator $A_{i}(\omega)$ is defined as

$A_{i}(\omega)=\sum_{E}\Pi(E)(a_{i}+a^{{\dagger}}_{i})\Pi(E+\omega),$ (S25)

where $\Pi(E)$ is a projection the eigenstates of $H_{S}$ of energy E. So
$A_{i}(\omega)$ is a lowering operator: it is the part of
$a_{i}+a^{{\dagger}}_{i}$ that couples eigenstates of $H_{S}$ whose energies
differ by $\omega$. It is easy to verify that
$\sum_{\omega}A_{i}(\omega)=a_{i}+a^{{\dagger}}_{i}$. Note that
$a_{i}+a^{{\dagger}}_{i}$ can only create an energy difference of
$\omega=0,\pm 4J,\pm 8J$, and we can work out $A_{i}(\omega)$ for each case
explicitly.

It is straightforward to verify that $A_{i}(\omega)$ is geometrically local
and takes the form

$\displaystyle A_{i}(0)$
$\displaystyle=(a_{i}+a^{{\dagger}}_{i})\sum_{\sigma(2)}P_{\sigma(2)},\
A_{i}(0)=(a_{i}+a^{{\dagger}}_{i})\sum_{\sigma(2)}P_{\sigma(2)},$ (S26)
$\displaystyle A_{i}(+4J)$
$\displaystyle=(a_{i}+a^{{\dagger}}_{i})\sum_{\sigma(3)}P_{\sigma(3)},\
A_{i}(-4J)=(a_{i}+a^{{\dagger}}_{i})\sum_{\sigma(1)}P_{\sigma(1)},$ (S27)
$\displaystyle A_{i}(+8J)$
$\displaystyle=(a_{i}+a^{{\dagger}}_{i})\sum_{\sigma(4)}P_{\sigma(4)},\
A_{i}(-8J)=(a_{i}+a^{{\dagger}}_{i})\sum_{\sigma(0)}P_{\sigma(0)},$ (S28)

where $P_{\sigma(n)}$ denotes the projector onto different local
configurations $\sigma(n)$ around site $i$ with $n$ domain walls. The
following relationship is satisfied:
$A_{i}^{{\dagger}}(\omega)=A_{i}(-\omega)$.

In Eq. (S24), the Hamiltonian $H^{\prime}$ is the Lamb-shift Hamiltonian

$H^{\prime}=\sum_{\omega}\sum_{i,j}\chi_{i,j}(\omega)A^{{\dagger}}_{i}(\omega)A_{j}(\omega),$
(S29)

where the coupling $\chi_{i,j}(\omega)$ is given in terms of the Fourier
transform of the reservoir correlation functions [62]:

$\chi_{i,j}(\omega)=\Im{C_{i,j}(\omega)},\ \text{where
}C_{i,j}(\omega)\coloneqq\int^{\infty}_{0}dse^{i\omega s}\langle
B^{{\dagger}}_{i}(s)B_{j}(0)\rangle,$ (S30)

where $B^{{\dagger}}_{i}(s)$ denotes Heisenberg evolution under $H_{B}$ for
time $s$ and the expectation value is taken in the initial thermal state of
the bath. To simplify the expression, we make the assumption that $\langle
B^{{\dagger}}_{i}(s)B_{j}(0)\rangle\approx\langle
B^{{\dagger}}_{i}(s)\rangle\langle B_{j}(0)\rangle$ for $i\neq j$. In
particular, this is immediately satisfied if each site is coupled to its own
bath. Note that we also assumed $\langle B^{{\dagger}}_{i}(s)\rangle=\langle
B_{j}(0)\rangle=0$ in order to get to Eq. (S24) (see Ref. [62]). We find
$\chi_{i,j}(\omega)=0$ for $i\neq j$. Defining
$\chi_{i,i}(\omega)\coloneqq\chi_{i}(\omega)$, $H^{\prime}$ simplifies to

$\displaystyle H^{\prime}$
$\displaystyle=\sum_{\omega,i}\chi_{i}(\omega)A^{{\dagger}}_{i}(\omega)A_{i}(\omega)$
$\displaystyle=\sum_{i}\left(\sum_{\omega}\chi_{i}(\omega)\sum_{\sigma(n[\omega])}P_{\sigma(n[\omega])}\right)(a_{i}+a_{i}^{{\dagger}})^{2},$
(S31)

where $n[\omega]$ is the number of domain walls in the projected configuration
as in Eqs. (S26)-(S28). We therefore have a Hamiltonian contribution quadratic
in the bosonic operators. A Hamiltonian of this form will in general cause
dephasing of the quantum memory. However, such an effect will be exponentially
suppressed by the combination of two-photon drive and two-photon loss in Eq.
(8) of the main text.

In Eq. (S24), the dissipative rate $\gamma_{i,j}(\omega)$ is also expressed in
terms of the correlation functions as
$\gamma_{i,j}(\omega)=\real[C_{i,j}(\omega)]=\int^{\infty}_{-\infty}dse^{i\omega
s}\langle B^{{\dagger}}_{i}(s)B_{j}(0)\rangle$. We make use of the locality
assumption again so $\gamma_{i,j}=0$ for $i\neq j$. The master equation Eq.
(S24) simplifies to

$\displaystyle\frac{d\rho_{S}}{dt}=-i[H^{\prime},\rho_{S}]+\sum_{\omega\in\\{0,\pm
2J,\pm
4J\\}}\sum_{i}\gamma_{i}(\omega)\left(A_{i}(\omega)\rho_{S}A_{i}^{{\dagger}}(\omega)-\frac{1}{2}\\{A^{{\dagger}}_{i}(\omega)A_{i}(\omega),\rho_{S}\\}\right).$
(S32)

Imposing the Kubo-Martin-Schwinger (KMS) condition (i.e. assuming that the
bath is in thermal equilibrium [62]) leads to the detailed balance relation

$\gamma_{i}(\omega)=e^{\beta\omega}\gamma_{i}(-\omega).$ (S33)

for some temperature $\beta=1/(k_{B}T)$ set by the bath.

###  IV.2 Recovering the photonic-Ising dissipators

We now show that the photonic-Ising dissipators defined in Eq. (9) of the main
text are a special case of the dissipative part of Eq. (S32). Let us set
$\gamma(0)=\gamma(-4J)=\gamma(-8J)=\kappa_{1}/2$. This can be satisfied at low
temperature for a roughly constant density of states across this energy range.
Suppose we initialize the system $\rho_{S}(0)$ as a pure state corresponding
to a well-defined domain-wall configuration. Let us denote
$\mathcal{D}(L)[\rho_{S}]=L\rho_{S}L^{{\dagger}}-\frac{1}{2}\\{L^{{\dagger}}L,\rho_{S}\\}$,
then, for any site $i$ and time $t$, the density matrix $\rho_{S}(t)$
satisfies

$\mathcal{D}\left(a_{i}\sum_{\sigma}P_{\sigma}\right)[\rho_{S}]=\sum_{\sigma}\mathcal{D}(a_{i}P_{\sigma})[\rho_{S}],$
(S34)

where the sum is taken over some projector $P_{\sigma}$ that projects onto
distinct domain wall configuration $\sigma$. The same result holds when we
replace $a_{i}$ by $a_{i}^{{\dagger}}$ or $a_{i}+a_{i}^{{\dagger}}$. This is
because $\rho_{S}(t)$ evolved under Eq. (S32) is always a linear combination
of some $\outerproduct{\psi_{1}}{\psi_{2}}$, where
$\ket{\psi_{1}},\ket{\psi_{2}}$ might be different states but they are always
the same domain wall configuration defined by the parity misalignment. The
dissipative terms on site $i$ in Eq. (S32) now read

$\displaystyle\gamma_{i}(0)\mathcal{D}(A_{i}(0))[\rho_{S}]+\sum_{\omega=\pm
4J}\gamma_{i}(\omega)\mathcal{D}(A_{i}(\omega))[\rho_{S}]+\sum_{\omega=\pm
8J}\gamma_{i}(\omega)\mathcal{D}(A_{i}(\omega))[\rho_{S}]$
$\displaystyle=\frac{\kappa_{1}}{2}\sum_{\omega=0,\pm 4J,\pm
8J}\mathcal{D}(A_{i}(\omega))[\rho_{S}]+\left(\gamma_{i}(2J)-\frac{\kappa_{1}}{2}\right)\mathcal{D}(A_{i}(4J))[\rho_{S}]+\left(\gamma_{i}(8J)-\frac{\kappa_{1}}{2}\right)\mathcal{D}(A_{i}(8J))[\rho_{S}]$
$\displaystyle=\frac{\kappa_{1}}{2}\mathcal{D}\left(\sum_{\omega=0,\pm 4J,\pm
8J}A_{i}(\omega)\right)[\rho_{S}]+\frac{\kappa_{1}}{2}(e^{4\beta
J}-1)\mathcal{D}(A_{i}(4J))[\rho_{S}]+\frac{\kappa_{1}}{2}(e^{8\beta
J}-1)\mathcal{D}(A_{i}(8J))[\rho_{S}]$
$\displaystyle=\frac{\kappa_{1}}{2}\mathcal{D}(a_{i}+a_{i}^{{\dagger}})[\rho_{S}]+\frac{\kappa_{1}}{2}(e^{4\beta
J}-1)\sum_{\sigma(3)}\mathcal{D}((a_{i}+a^{{\dagger}}_{i})P_{\sigma(3)})[\rho_{S}]+\frac{\kappa_{1}}{2}(e^{8\beta
J}-1)\sum_{\sigma(4)}\mathcal{D}((a_{i}+a^{{\dagger}}_{i})P_{\sigma(4)})[\rho_{S}].$
(S35)

We can now choose
$\beta=\frac{1}{8J}\ln\left[\frac{\kappa_{nn}+\kappa_{1}}{\kappa_{1}}\right]$
for some $\kappa_{nn}$. Then we have

$\frac{\kappa_{1}}{2}\mathcal{D}(a_{i}+a_{i}^{{\dagger}})[\rho_{S}]+\frac{\tilde{\kappa}_{nn}}{2}\sum_{\sigma(3)}\mathcal{D}((a_{i}+a^{{\dagger}}_{i})P_{\sigma(3)})[\rho_{S}]+\frac{\kappa_{nn}}{2}\sum_{\sigma(4)}\mathcal{D}((a_{i}+a^{{\dagger}}_{i})P_{\sigma(4)})[\rho_{S}],$
(S36)

where
$\tilde{\kappa}_{nn}=\sqrt{\kappa_{1}\kappa_{nn}+\kappa_{1}^{2}}-\kappa_{1}$.
This is almost the same as the dissipators defined in Eq. (9) of the main
text. The first difference is that all jumps here are proportional to
$a_{i}+a_{i}^{\dagger}$, while in the main text all jumps are proportional to
just $a_{i}$. (We shall address this in the next paragraph.) The second
difference is that here $\kappa_{1}$ came solely from the bath that we
designed ourselves, representing the fluctuations associated with a non-zero-
temperature bath that disorders the state. We note that single-photon loss
processes can also occur due to spontaneous emission for a single cavity,
which arises due to processes that lower the energy within a single cavity
(not included in the analysis above for simplicity). Taking such terms into
account amounts to shifting up the effective temperature of the model.
However, as long as the composite system sits in the low-temperature part of
the phase diagram (i.e. if the “zero-temperature” part of the bath is
significantly stronger than the terms that generate bit flips), the system
will sit in the part of the phase diagram that hosts a quantum memory.

To recover the photonic-Ising dissipators, we take into account the two-photon
process on each photonic-Ising site. The two-photon process effectively
constrains the bosonic Hilbert space on each site to a two-dimensional
manifold, making the thermal equilibrium of $H_{S}$ well-defined. Without the
two-photon process, $H_{S}$ is infinitely degenerate even for a finite
photonic-Ising lattice, which leads to an ill-defined steady-state manifold
for the Lindbladian even when detailed balance is enforced.

To simplify the analysis, we assume that $N=|\alpha|^{2}\rightarrow\infty$ and
$\kappa_{2}$ is nonzero, in which case the system is at all times close to the
ideal manifold spanned by the two coherent states $\ket{\pm\alpha}=\ket{\pm
e^{-i\pi/4}\sqrt{N}}$ with $\langle\alpha|{-\alpha}\rangle=0$. So, within this
two-dimensional Hilbert space, both $a$ and $a^{\dagger}$ have vanishing
diagonal matrix elements, while their off-diagonal matrix elements are off by
only a phase, so we can write $a+a^{{\dagger}}\rightarrow(1+i)a$. This finally
yields

$\kappa_{1}\mathcal{D}(a_{i})[\rho_{S}]+\tilde{\kappa}_{nn}\sum_{\sigma(3)}\mathcal{D}(a_{i}P_{\sigma(3)})[\rho_{S}]+\kappa_{nn}\sum_{\sigma(4)}\mathcal{D}(a_{i}P_{\sigma(4)})[\rho_{S}],$
(S37)

which are the photonic-Ising dissipators [Eq. (9) in the main text] and the
single-photon loss noise at site $i$.

Note that, while we made many assumptions to arrive at the precise dissipators
from Eq. (9) in the main text, it is likely that nearly any low-temperature
thermal Markovian bath will result in a protected quantum memory.

##  V 5\. Engineering an Ising interaction between cavity modes

Figure S5:  (a) An $LC$ oscillator (representing a cavity) connected to a
Josephson junction. The junction generates a Hamiltonian term that is
proportional to the photon parity in the cavity. (b) Two $LC$ oscillators
connected via a Josephson junction. The junction generates an Ising-like
parity-parity interaction between the oscillators in a certain parameter
regime. (c) Schematic layout for a 2D quantum memory: Black circles represent
driven-dissipative resonator cavities (or $LC$ oscillators). Each resonator is
coupled to its neighbor via a Josephson junction which generates an Ising-like
interaction (in the appropriate limit).

We showed in the last section that the dissipators described in Eq. (9) most
naturally appear as thermal dissipators corresponding to processes that raise
and lower the energy of the Ising-like Hamiltonian $H=-J\sum_{\langle
ij\rangle}Q_{i}Q_{j}$, where $Q_{i}=e^{i\pi a_{i}^{\dagger}a_{i}}$. The most
straightforward approach to achieve a passively protected memory is therefore
to engineer an Ising-like interaction between nearest-neighbor cat qubits. At
low temperatures, the coupling of the system to its thermal environment can
lead to thermal processes that drive the system to a ferromagnetic state, thus
protecting the system against bit flips.

Previous studies have described how to achieve a parity-parity interaction
term between neighboring cat qubits [56, 55]. Here we first review the steps
needed to achieve a single-cavity Hamiltonian proportional to parity, closely
following Ref. [56]. We then discuss the generalization to a parity-parity
interaction which follows in a very similar manner.

Consider an LC oscillator (representing a cavity mode) connected to a
Josephson junction, as shown in Fig. S5(a), in the presence of a two-photon
drive on the cavity (not shown). The Hamiltonian of the system reads

$H=\frac{\hat{q}^{2}}{2C}+\frac{\hat{\phi}^{2}}{2L}-E_{J}\cos\left(\hat{\phi}/\Phi_{0}\right)+\lambda\left(\hat{a}^{2}e^{i\omega_{d}t}+(\hat{a}^{\dagger})^{2}e^{-i\omega_{d}t}\right),$
(S38)

where $\hat{q}$ is the charge of the capcitor $C$, $\hat{\phi}$ is the flux
through the inductor $L$, and $E_{J}$ is the Josephson energy,
$\Phi_{0}=\hbar/(2e)$ is the flux quantum, $\hat{a}$ annihilates an excitation
of the cavity, $\lambda$ is the drive strength, and $\omega_{d}$ is the drive
frequency. We note that the flux through the inductor is equal in magnitude to
the flux through the Josephson junction since Kirchoff’s law relates the
voltage drop across the two elements at any given time. We can rewrite the
Hamiltonian in terms of $a,a^{\dagger}$:

$H=\hbar\omega(a^{\dagger}a+1/2)-E_{J}\cos{[x(a+a^{\dagger})]}+\lambda\left(a^{2}e^{i\omega_{d}t}+(a^{\dagger})^{2}e^{-i\omega_{d}t}\right),$
(S39)

where $\phi=\sqrt{\hbar
Z/2}(a+a^{\dagger}),q=(1/i)\sqrt{\hbar/(2Z)}(a-a^{\dagger}),\omega=1/\sqrt{LC},Z=\sqrt{L/C},x=\Phi_{0}^{-1}\sqrt{\hbar
Z/2}$ and we drop the hats on operators henceforth.

Going to the interaction picture with respect to
$\hbar\omega(a^{\dagger}a+1/2)$, the Hamiltonian takes the form

$H_{int}=-E_{J}\cos{[x(ae^{-i\omega t}+a^{\dagger}e^{+i\omega
t})]}+\lambda(a^{2}+(a^{\dagger})^{2})=-\frac{E_{J}}{2}(D[\beta(t)]+D[-\beta(t)])+\lambda(a^{2}+(a^{\dagger})^{2}),$
(S40)

where we have assumed an on-resonant drive, $\omega_{d}=2\omega$, and defined
the displacement operator $D[\beta(t)]=e^{\beta(t)a^{\dagger}-\beta^{*}(t)a}$,
with $\beta(t)=ixe^{-i\omega t}$. We apply the rotating-wave approximation to
remove all time dependence in the Hamiltonian, which is valid in the limit:
$\omega\gg E_{J}$:

$H_{rw}=-E_{J}e^{-x^{2}/2}\sum_{n}\left(L_{n}(x^{2})|n\rangle\langle
n|\right)+\lambda(a^{2}+(a^{\dagger})^{2}),$ (S41)

where $L_{n}$ is the Laguerre polynomial of order $n$.

The two-photon drive (along with engineered two-photon loss) ensures that the
system is approximately confined to a two-dimensional manifold spanned by the
cat states $|\alpha_{e}\rangle$ and $|\alpha_{o}\rangle$. In this subspace,
the Hamiltonian is diagonal (since $H_{rw}$ is diagonal in the Fock basis).
Further, if we set $x=2|\alpha|$, then the Hamiltonian is exponentially close
to the parity operator:

$\displaystyle H_{rw}^{\text{cat}}$
$\displaystyle=-\frac{\hbar\Omega}{2}(|\alpha_{e}\rangle\langle\alpha_{e}|-|\alpha_{o}\rangle\langle\alpha_{o}|)+O(E_{J}e^{-|\alpha|^{2}/2})$
(S42)
$\displaystyle=-\left(\frac{\hbar\Omega}{2}\right)Q+O(E_{J}e^{-|\alpha|^{2}/2}),\qquad
Q=e^{i\pi a^{\dagger}a},$ (S43)

where $\Omega=E_{J}/(\hbar\sqrt{2\pi|\alpha|^{2}})$. The exponentially small
term is proportional to the identity in this subspace, which should not affect
the dynamics.

Let us briefly summarize the physical conditions required to achieve a parity
Hamiltonian of the form $H\sim Q$ in (S43) above. In the large
$|\alpha|^{2}\equiv N$ limit, we require that $E_{J}\sim\sqrt{N}$ such that
the coefficient $\Omega$ does not tend to zero at large $N$. In order for the
rotating wave approximation to be valid, we required that $\omega\gg E_{J}$.
This implies that $\omega\sim\sqrt{N}$, which can be achieved with a small
capacitance: $\omega\sim 1/\sqrt{LC}$ if $C\sim N^{-1}$. Finally, we required
$x\sim\sqrt{L/C}\sim\sqrt{N}$. This is again satisfied with the small
capacitance condition: $C\sim N^{-1}$.

So far, our discussion has focused on achieving a parity Hamiltonian for a
single cavity. A very similar setup will result in a parity-parity interaction
between neighboring cavity modes [56]. We briefly describe how this can be
done.

Consider two driven LC oscillators which are connected to a Josephson
junction, as shown in Fig. S5(b). (We again assume two-photon drives on each
cavity as before, but neglect to include them in the Hamiltonian since they
only serve the purpose of confining the state of the system to the cat state
subspace, as described in the previous paragraphs.) The Hamiltonian for the
system reads

$H=\left(\frac{q_{1}^{2}}{2C_{1}}+\frac{\phi_{1}^{2}}{2L_{1}}\right)+\left(\frac{q_{2}^{2}}{2C_{2}}+\frac{\phi_{2}^{2}}{2L_{2}}\right)-E_{J}\cos(\phi_{2}+\phi_{1}),$
(S44)

where $\phi_{1/2}$ are the node fluxes defined in Fig. S5(b). Moving to the
rotating frame of the cavity Hamiltonians leads to the interaction Hamiltonian

$H_{int}=-E_{J}\cos\left((a_{2}e^{-i\omega_{2}t}+a_{2}^{\dagger}e^{i\omega_{2}t})+(a_{1}e^{-i\omega_{1}t}+a_{1}^{\dagger}e^{i\omega_{1}t}))\right),$
(S45)

where $a_{1/2}$ are the annihilation operators associated with the two
cavities, and $\omega_{1/2}$ are the frequencies. Applying the rotating-wave
approximation leads to

$H_{rw}=-E_{J}e^{-(x_{1}^{2}+x_{2}^{2})/2}\sum_{n_{1},n_{2}}L_{n_{1}}(x_{1}^{2})L_{n_{2}}(x_{2}^{2})|n_{1},n_{2}\rangle\langle
n_{1},n_{2}|.$ (S46)

Note that we also require that the frequencies of the cavities should be
incommensurate in order for the terms above to be the only ones that are time
independent, i.e. $l_{1}\omega_{1}\neq l_{2}\omega_{2},\forall
l_{1},l_{2}\in\mathbb{Z}$. Specializing to the two-dimensional cat state
manifold leads to an interaction of the form

$H_{rw}^{\text{cat}}=-\left(\frac{\Omega_{1}\Omega_{2}}{4E_{J}}\right)Q_{1}Q_{2},\qquad
Q_{1/2}=e^{i\pi a_{1/2}^{\dagger}a_{1/2}},$ (S47)

where $\Omega_{1/2}=E_{J}/(\sqrt{2\pi|\alpha_{1/2}|^{2}})$ set the energy
scale of the coupling between the two cavities $1,2$.

Again, let us briefly summarize the physical conditions required to achieve
the parity-parity Hamiltonian of the form $H\sim Q_{1}Q_{2}$ in Eq. (S47)
above. In the large $|\alpha|^{2}\equiv N$ limit, we require that $E_{J}\sim
N$ such that the coefficient $\sim\Omega_{1}\Omega_{2}/E_{J}$ does not tend to
zero at large $N$. In order for the rotating wave approximation to be valid,
we required that $\omega_{1/2}\gg E_{J}$ and $|\omega_{1}-\omega_{2}|\gg
E_{J}$. This implies that $\omega\sim N$. We also require that
$x_{1/2}\sim\sqrt{L_{1/2}/C_{1/2}}\sim\sqrt{N}$. Both of these conditions can
be achieved for $C_{1/2}\sim N^{-3/2},L_{1/2}\sim N^{-1/2}$.

A schematic layout for a 2D quantum memory is provided in Fig. S5(c): Each
black dot represents a driven-dissipative resonator, connected to its neighbor
via a Josephson junction that results in an Ising interaction. The advantage
of the passive approach is the lack of ancilla qubits and of precise pulse
signals that are typically needed to make the measurements required for active
error correction. The passive approach also avoids the need for classical
communication with a decoder.

The scheme we have just described to achieve a parity-parity interaction
requires large tunability of the frequency and impedance of a cavity, and of
the Josephson tunneling coefficient of the junction. These conditions may
prove to be challenging to engineer in practice, and it is still unclear
whether such a protocol is the simplest way to arrive at the effective model
described in the main text. Any circuit QED proposal to arrive at our model
should involve nonlinear devices, since the desired parity-parity interaction
is highly nonlinear. In this section, we have explored a conceptually simple
route to achieve this nonlinearity via a single Josephson junction. It is
conceivable that adding other devices would relax the conditions required to
achieve a substantial Ising interaction, necessary for an exponentially
suppressed logical bit flip rate. We leave such investigations to future work.

##  VI 6\. Digital autonomous photonic-Ising local decoder

In this section, we describe a digital autonomous approach for realizing a
stochastic local error decoder inspired by the photonic-Ising dissipators
defined in the main text. This procedure involves implementing a sequence of
fault-tolerant local gates to correct the errors without the need of
measurements. The step is then iterated over time on the entire system as fast
as possible. Although this approach is different from directly realizing the
microscopic Lindbladian (which could be done by dividing the error-correcting
step infinitesimally [63, 64]), this digital approach is potentially easier to
realize experimentally, and we expect that it provides the same dynamical
protection of the quantum memory. Note that the protocol requires the rate of
the digital steps to scale linearly with the average photon number.

It suffices to consider the implementation of the local decoder at a single
time step. Consider a square lattice of photonic cavities. Over each lattice
cavity, we place an ancillary cavity (initialized in $\ket{\alpha_{o}}$).

We describe the implementation at a single site and the generalization to the
whole lattice follows straightforwardly. First, for a given site, we perform
an encoding unitary $U$. Depending on the state of the chosen lattice cavity
and its four neighbors, $U$ changes the state of the ancillary cavity via

$U=P\otimes(\outerproduct{\alpha_{e}}{\alpha_{o}}+\outerproduct{\alpha_{o}}{\alpha_{e}})+P^{\perp}\otimes(\outerproduct{\alpha_{o}}{\alpha_{o}}+\outerproduct{\alpha_{e}}{\alpha_{e}}),$
(S48)

where $P$ projects on a local configuration of domain walls (a specific
example was considered in the main text, where we project onto a configuration
with 3 or 4 misaligned neighboring lattice cavities) and $P^{\perp}$ is the
orthogonal subspace projector. We define the unitary $U$ such that it changes
the ancillary cavity from $\ket{\alpha_{o}}$ to $\ket{\alpha_{e}}$ if a local
error is detected; it does nothing otherwise. Note that $U$ can be implemented
using the fundamental set of bias-preserving gates in Ref. [54], where the
two-photon drive and two-photon loss can be kept on thus suppressing the
dephasing errors during the gate implementation.

Second, we apply a CNOT gate (described in Section IV.D in Ref. [54]) that is
controlled by the ancillary cavity and targets the corresponding lattice
cavity. (We use the convention that $\ket{\alpha_{o}}$ is $\ket{0}$ and
$\ket{\alpha_{e}}$ is $\ket{1}$).

Third, we use a strong dispersive coupling to a transmon to extract the
entropy from the ancillary cavity [41] and reset it back to the initial state
$\ket{\alpha_{o}}$. Due to the fault-tolerance of the cavity-cavity gates
[54], the phase errors stay suppressed when the Ising-type local decoder is
implemented autonomously.

The full procedure is achieved by implementing the encoding and the reset
across the entire lattice. To extend the single-site procedure to the entire
lattice, we note that the encoding operations $U$ on each site are local
around each lattice cavity and they commute across different lattice sites.
Therefore, the encoding $U$ can be implemented in parallel across all the
sites before a final reset, e.g. by dividing the lattice into bipartite
sublattices and operating on the cavities that belong to the same sublattice
in parallel.

For concreteness, let us consider an example where $U$ flips the ancillary
cavity if it identifies a corner formed by the domain walls, i.e. Toom’s rule.
(Digitally, this is simpler than the majority rule described in the main text,
but one could implemenent the majority rule approach as well.) Then $U$ for a
chosen orientation can be implemented by the following gate sequence

  1. 1.

Apply two CNOTs from center lattice cavity to its neighboring cavities on the
left and on the top.

  2. 2.

Apply a Toffoli gate (described in Section IV.E in Ref. [54]) controlled by
the two neighboring cavities and targeting the ancillary cavity.

  3. 3.

Repeat step 1 to invert the two CNOTs applied.

The digital procedure for the autonomous implementation of Toom’s rule on a
particular configuration is schematically depicted in Fig. S6.

Figure S6: Illustration of the digital autonomous photonic-Ising local decoder
when $U$ is implementing Toom’s rule. Purple color labels the ancillary
cavity. We denote $\ket{\alpha_{o}}$ by ($-$) and $\ket{\alpha_{e}}$ by ($+$).
The green arrows are CNOT gates, pointing from the control cavity to the
target cavity. A connected double green arrow is a Toffoli gate. The procedure
consists of the following steps: (i) two CNOTs from the central cavity to its
neighbors; (ii) a Toffoli gate from the neighbors to the ancillary cavity;
(iii) two CNOTs from the central cavity to its neighbors; (iv) a CNOT from the
ancillary cavity to the central cavity; (v) reset the ancillary cavity via
coupling to a transmon.

The autonomous approach above can be easily turned into an active error
correction protocol: instead of applying a CNOT from the ancillary cavity to
the central cavity and then resetting the ancillary cavity, we can measure the
ancillary cavity and flip the parity of the central cavity if the measurement
result is ($+$). Alternatively, when implementing active error correction, we
can place a syndrome cavity between each pair of neighboring cavities on the
lattice. By storing in the syndrome cavity the information regarding the
presence of a domain wall, we can implement a local decoder based on the rules
defined by $U$. This then becomes a 2D version of a repetition cat code.
Again, all the steps can be achieved with dephasing errors exponentially
suppressed. In contrast with the non-local processing of syndrome information
required by the 1D repetition cat code in Ref. [54], the 2D code allows for a
stochastic local decoding procedure. As mentioned at the beginning of the
section, to achieve an exponentially long memory time, both the autonomous and
the active error correction approaches require scaling the rate of the digital
step linearly with the average number of cavity photons $N$ because bit-flip
error rate scales with $N$.