# Mixed-State Quantum Spin Liquids and Dynamical Anyon Condensations in Kitaev
Lindbladians

Kyusung Hwang  School of Physics, Korea Institute for Advanced Study, Seoul
02455, Korea<EMAIL_ADDRESS>
###### Abstract

Quantum spin liquids and anyons, used to be subjects of condensed matter
physics, now are realized in various platforms of qubits, offering
unprecedented opportunities to investigate fundamental physics of many-body
quantum entangled states. Qubits are inevitably exposed to environments’
effects such as decoherence and dissipation, which are believed to be
detrimental to many-body entanglement. Here, we argue that unlike the common
belief decoherence and dissipation can give rise to novel topological
phenomena in quantum spin liquids. We study open quantum systems of the Kitaev
spin liquid and the toric code via the Lindblad master equation approach. By
using exact solutions and numerical approaches, we show the dynamical
occurrence of anyon condensation by decoherence and dissipation, which results
in a topological transition from the initial state spin liquid to the steady
state spin liquid. The mechanism of the anyon condensation transition by the
Lindblad dynamics is elucidated. We also provide an insight into the
relationship between the Kitaev spin liquid and the toric code in the picture
of anyon condensation. Our work suggests open quantum systems to be a new
venue for topological phenomena of quantum spin liquids and anyons.

##  1 Introduction

Quantum spin liquids are exotic phases of matter featured with emergent anyon
quasiparticles and long range entanglement [1, 2, 3, 4, 5, 6, 7, 8, 9, 10,
11]. In addition to their importance in fundamental physics, anyons and
quantum spin liquids have been a central topic of recent studies due to their
potential applications in quantum memories and quantum computations [12, 2, 1,
13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31].
In the aspect of applications, it is essential to understand quantum devices
under the unavoidable influences of environments such as decoherence and
dissipation effects [32]. Studies of such open quantum systems also provide a
promising opportunity for exploring new phenomena fostered by many-body
quantum entanglement and the environment effects [33, 34, 35]. Nonetheless,
most studies of quantum spin liquids have focussed on closed systems so far.
Open quantum spin liquids remain largely unexplored due to the complexities of
the problems and the absence of solvable systems.

In this work, we introduce open quantum spin liquids that are exactly solved
in the steady state limit. As a primary example, we consider the Kitaev spin
liquid (KSL) [1] coupled with a Markovian environment via the Lindblad master
equation approach [36, 37, 38, 39, 40, 41, 42]. Under the Lindblad dynamics,
the pure state of the KSL evolves to a steady state that shows vanishing spin-
spin correlation, yet still preserving the zero-flux quantum number. The exact
form of the steady state is given by the maximally mixed state in the zero-
flux sector with equal weight. We call this state “mixed-state Kitaev spin
liquid”.

| Original system | Doubled system
---|---|---
State | $\hat{\rho}$ | $|\rho\rangle\\!\rangle$
Time evolution | $\frac{d}{dt}\hat{\rho}=\mathcal{L}(\hat{\rho})$ | $i\frac{d}{dt}|\rho\rangle\\!\rangle=\mathpzc{H}|\rho\rangle\\!\rangle$

Figure 1: Open quantum system of the Kitaev spin liquid in two equivalent
descriptions. (a) Original density matrix description. The Kitaev honeycomb
model is embedded in the Lindblad master equation with the jump operator,
$L_{\mu}=\sqrt{\gamma}\sigma_{j}^{\alpha}\sigma_{k}^{\alpha}$. (b) Doubled
state vector description. Vectorization of the Lindblad master equation leads
to a closed system defined on the bilayer honeycomb lattice with the non-
Hermitian interlayer interaction,
$i\gamma(\sigma_{j}^{\alpha}\sigma_{k}^{\alpha}\tau_{j}^{\alpha}\tau_{k}^{\alpha}-1)$.
The two descriptions are compared in the table.

The Lindblad dynamics of the Kitaev spin liquid presents rich physics with a
deep connection to anyon condensation, which is an important mechanism for
understanding topological transitions between distinct quantum spin liquids
[13, 43, 44, 45, 46, 47]. By mapping the Lindblad system to a doubled non-
Hermitian system [48, 49], we uncover that the time evolution from the pure
KSL to the mixed-state KSL is actually a dynamical transition from a KSL
bilayer product state to a toric code type $\mathbb{Z}_{2}$ spin liquid state
by anyon condensation.

Our study illuminates a new aspect of open quantum systems, i.e., a useful
platform for studying topological transitions in quantum spin liquids. In
addition to the Kitaev spin liquid, we study open quantum systems of the toric
code model [2]. Based on the concrete examples, we elucidate the mechanism of
dynamical anyon condensation transitions in quantum spin liquids induced by
the Lindblad dynamics.

##  2 Kitaev Lindbladian

We consider an open quantum system of the Kitaev spin liquid coupled to a
Markovian environment. The dynamics of the system is investigated by using the
Lindblad master equation [36, 37, 38, 39, 40, 41, 42],

$\frac{d\hat{\rho}}{dt}=\mathcal{L}(\hat{\rho})=-i[H,\hat{\rho}]+\sum_{\mu}(L_{\mu}\hat{\rho}L_{\mu}^{\dagger}-\frac{1}{2}\\{L_{\mu}^{\dagger}L_{\mu},\hat{\rho}\\}),$
(1)

where the system’s density matrix $\hat{\rho}(t)$ evolves in time by the
Lindbladian $\mathcal{L}$ composed of the Hamiltonian $H$ and the Lindblad
operators $\\{L_{\mu}\\}$ describing the environment effects. The Lindblad
equation has the completely positive and trace-preserving property, i.e.,
${\rm Tr}\hat{\rho}(t)={\rm const}$. For the Hamiltonian, we consider the
Kitaev honeycomb model [1],

$H=K\sum_{\langle jk\rangle_{x}}\sigma_{j}^{x}\sigma_{k}^{x}+K\sum_{\langle
jk\rangle_{y}}\sigma_{j}^{y}\sigma_{k}^{y}+K\sum_{\langle
jk\rangle_{z}}\sigma_{j}^{z}\sigma_{k}^{z},$ (2)

where $\sigma^{x,y,z}$ mean the Pauli matrices and $\langle jk\rangle_{x,y,z}$
denote $x,y,z$-bonds of the honeycomb lattice [Fig. 1(a)]. For the Lindblad
operators generating non-unitary dynamics, we assume the Kitaev bond
interactions,

$L_{\mu}=\sqrt{\gamma}\sigma_{j}^{\alpha}\sigma_{k}^{\alpha},$ (3)

with the dissipation strength $\gamma\leavevmode\nobreak\ (\geq 0)$. The
Lindblad operator generates decoherence and dissipation effects on the Kitaev
spin liquid state by creating a pair of fermion excitations at both sites of
the bond $\langle jk\rangle_{\alpha}$. The Hamiltonian $H$ and the Lindblad
operators $\\{L_{\mu}\\}$ both commute with the $\mathbb{Z}_{2}$ flux
operator,

$\hat{W}_{p}=\sigma_{1}^{z}\sigma_{2}^{y}\sigma_{3}^{x}\sigma_{4}^{z}\sigma_{5}^{y}\sigma_{6}^{x},$
(4)

i.e., $[H,\hat{W}_{p}]=[L_{\mu},\hat{W}_{p}]=0$. Hence, the Kitaev Lindbladian
$\mathcal{L}$ conserves the $\mathbb{Z}_{2}$ flux quantum number.

Figure 2: Lindblad dynamics of the Kitaev spin liquid. (a) The trace of the
density matrix ${\rm Tr}\hat{\rho(t)}$ and the expectation value of the flux
operator $\langle\hat{W}_{p}\rangle$. (b) The Rényi entropy
$S_{\textup{R\'{e}nyi}}$. The dashed line marks the steady state value,
$\log(2^{N/2-1})$, where $N=24$ is the number of sites. (c) The spin-spin
correlator $\langle{\sigma}_{j}^{\alpha}{\sigma}_{k}^{\alpha}\rangle$ for the
nearest-neighbor bond $\langle jk\rangle_{\alpha}$. Different colors denote
the results of different dissipation strengths ($\gamma=0,\leavevmode\nobreak\
0.2,\leavevmode\nobreak\ 0.4,\leavevmode\nobreak\ 0.6,\leavevmode\nobreak\
0.8,\leavevmode\nobreak\ 1$). The unitary evolution of the closed system
($\gamma=0$) is shown together for comparison (black). In all the
calculations, the coupling constant is fixed by $K=-1$.

###  2.1 Lindblad dynamics of the Kitaev spin liquid

To understand the dynamics generated by many-body quantum entanglement and the
environment effects in an unbiased way, we numerically solve the Lindblad
equation by putting the system on a 24-site cluster with periodic boundary
condition [Fig. 1(a)]. First, we prepare the system in the pure state of the
Kitaev spin liquid,

$\hat{\rho}(t=0)=|\Psi_{\rm KSL}\rangle\langle\Psi_{\rm KSL}|,$ (5)

where $|\Psi_{\rm KSL}\rangle$ is the ground state of $H$ which we obtain by
exact diagonalization on the 24-site cluster (torus geometry). The model $H$
has threefold ground state degeneracy, among which we choose the ground state
$|\Psi_{\rm KSL}\rangle$ in the sector of the Wilson loop flux
$\mathcal{W}_{x,y}=-1$ in our calculations (see Appendix A). The time
evolution of the density matrix,
$\hat{\rho}(t)=e^{t\mathcal{L}}\hat{\rho}(0)$, is computed by using the Krylov
subspace methods [50, 51]. In the construction of the time evolution operator
$e^{t\mathcal{L}}$, we utilize the $\mathbb{Z}_{2}$ flux quantum numbers of
the Kitaev Lindbladian to reduce the size of the Hilbert space.

Figure 3: Time evolution of the spin structure factor $S({\bf q})$. The three
plots at $t=0,0.1,0.2$ are obtained for the case of $K=-1$ and $\gamma=1$. In
each plot, the two hexagons represent the first and second Brillouin zones in
momentum space.

Figure 2 shows the calculation results of the Lindblad dynamics. First, we
check the properties ${\rm Tr}\hat{\rho}(t)=1$ and
$\langle\hat{W}_{p}\rangle={\rm Tr}[\hat{\rho}(t)\hat{W}_{p}]=1$ in Fig. 2(a).
Next, we consider the Rényi entropy,

$S_{\textup{R\'{e}nyi}}=-{\rm log}[{\rm Tr}\hat{\rho}^{2}(t)],$ (6)

which tells us how the system’s purity changes by the environment effects. As
shown in Fig. 2(b), the Rényi entropy gradually increases converging to the
constant value: $S_{\textup{R\'{e}nyi}}=\log(2^{N/2-1})$ when
$t\rightarrow\infty$ ($N=24$ is the number of sites). The larger the
dissipation strength $\gamma$ is, the faster the convergence is reached. The
constant value reveals the size of the subspace ($2^{N/2-1}=2,048$) that the
system belongs to, i.e., the zero flux sector with the Wilson loop flux
$\mathcal{W}_{x,y}=-1$ (see Appendix A).

An interesting behavior appears in the spin-spin correlation
$\langle{\sigma}_{j}^{\alpha}{\sigma}_{k}^{\alpha}\rangle={\rm
Tr}[\hat{\rho}(t){\sigma}_{j}^{\alpha}{\sigma}_{k}^{\alpha}]$. It
monotonically decreases and vanishes in the long time limit:
$\langle{\sigma}_{j}^{\alpha}{\sigma}_{k}^{\alpha}\rangle=0$ when
$t\rightarrow\infty$ [Fig. 2(c)]. The system is still in a quantum spin liquid
phase since it stays in the zero flux sector. Nonetheless, the nearest
neighbor bond spin-spin correlation disappears in the long time limit. This
peculiar property is also reflected in the spin structure factor,

$S({\bf q})=\frac{1}{N}\sum_{j,k}e^{i{\bf q}\cdot({\bf r}_{j}-{\bf
r}_{k})}\langle{\bm{\sigma}}_{j}\cdot{\bm{\sigma}}_{k}\rangle,$ (7)

which becomes flat in momentum space in the long time limit: $S({\bf q})=3$
when $t\rightarrow\infty$ (see Fig. 3).

###  2.2 Exact solution of the steady state

The numerical approach establishes the existence of the steady state
characterized by the properties, $\langle\hat{W}_{p}\rangle=1$,
$S_{\textup{R\'{e}nyi}}=\log(2^{N/2-1})$, and
$\langle{\sigma}_{j}^{\alpha}{\sigma}_{k}^{\alpha}\rangle=0$. We find the
exact solution for the steady state:

$\displaystyle\hat{\rho}_{\rm MSKSL}$ $\displaystyle=$
$\displaystyle\frac{1}{2^{N+1}}\prod_{l=x,y}(1-\hat{\mathcal{W}}_{l})\prod_{p}(1+\hat{W}_{p})$
(8) $\displaystyle=$
$\displaystyle\sum_{n}\rho_{n}|\Psi_{n}\rangle\langle\Psi_{n}|.$

The steady state density matrix is simply given by the projection operator
into the zero flux sector with the Wilson loop flux $\mathcal{W}_{x,y}=-1$.
Hence, it is diagonal in the eigenbasis $\\{|\Psi_{n}\rangle\\}$ of $H$ and
$\hat{W}_{p}$ with the weight

$\rho_{n}=\left\\{\begin{array}[]{cc}\frac{1}{2^{N/2-1}}&(W_{p}=+1\leavevmode\nobreak\
\&\leavevmode\nobreak\ \mathcal{W}_{x,y}=-1)\\\ 0&({\rm
otherwise})\end{array}\right..$ (9)

This equal weight property is indeed identified in our numerical calculations
shown in Fig. 4. One can check that $\mathcal{L}(\hat{\rho}_{\rm MSKSL})=0$
along with the three properties mentioned earlier. This is the maximally mixed
state in the zero-flux sector with equal weight, which we call “mixed-state
Kitaev spin liquid (MSKSL)”.

##  3 Mapping to the non-Hermitian bilayer model

Lindblad system allows an exact mapping to a doubled system described by a
non-Hermitian Hamiltonian [48, 49]. The mapping is conducted by the
vectorization of the density matrix called the Choi-Jamiołkowski isomorphism
[52, 53]:

$\displaystyle\hat{\rho}=\sum_{m,n}\rho_{mn}|m\rangle\langle
n|\leavevmode\nobreak\ \leavevmode\nobreak\ \Rightarrow\leavevmode\nobreak\
\leavevmode\nobreak\
|\rho\rangle\\!\rangle=\sum_{m,n}\rho_{mn}|m\rangle\otimes|n\rangle,\leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ $ (10)

where the bra of the density matrix $\hat{\rho}$ is turned into a ket of the
resulting state vector $|\rho\rangle\\!\rangle$. The system is now effectively
doubled with an additional copy of the original Hilbert space. By the
vectorization, the Lindblad master equation takes the form of the Schrödinger
equation

$\displaystyle
i\frac{d|\rho\rangle\\!\rangle}{dt}=\mathpzc{H}|\rho\rangle\\!\rangle$ (11)

with the non-Hermitian Hamiltonian

$\displaystyle\mathpzc{H}$ $\displaystyle=$ $\displaystyle H\otimes I-I\otimes
H^{T}$ $\displaystyle+$ $\displaystyle i\sum_{\mu}\left[L_{\mu}\otimes
L_{\mu}^{*}-\frac{1}{2}(L_{\mu}^{\dagger}L_{\mu})\otimes
I-\frac{1}{2}I\otimes(L_{\mu}^{\dagger}L_{\mu})^{*}\right].$

In our case, the non-Hermitian Hamiltonian is given by

$\displaystyle\mathpzc{H}$ $\displaystyle=$ $\displaystyle K\sum_{\langle
jk\rangle_{\alpha}}^{\textup{upper
layer}}\sigma_{j}^{\alpha}\sigma_{k}^{\alpha}-K\sum_{\langle
jk\rangle_{\alpha}}^{\textup{lower layer}}\tau_{j}^{\alpha}\tau_{k}^{\alpha}$
(13) $\displaystyle+$ $\displaystyle i\gamma\sum_{\langle
jk\rangle_{\alpha}}^{\textup{inter-
layer}}(\sigma_{j}^{\alpha}\sigma_{k}^{\alpha}\tau_{j}^{\alpha}\tau_{k}^{\alpha}-1),$

where $\sigma$ & $\tau$ are Pauli spin operators acting on the original kets
and copied kets;
$\sigma_{j}^{\alpha}\tau_{k}^{\beta}(|m\rangle\otimes|n\rangle)=(\sigma_{j}^{\alpha}|m\rangle)\otimes(\tau_{k}^{\beta}|n\rangle)$.
We now have a bilayer spin model with the intralayer Kitaev interactions ($K$
& $-K$) and the interlayer interactions ($i\gamma$) [Fig. 1(b)]. Note that the
interlayer interactions are non-Hermitian originating from the environment
effects. This leads to non-unitary time evolution of the state vector,
$|\rho(t)\rangle\\!\rangle=e^{-it\mathpzc{H}}|\rho(0)\rangle\\!\rangle$, which
is identified from the relationship
$\langle\\!\langle\rho|\rho\rangle\\!\rangle=\exp(-S_{\textup{R\'{e}nyi}})$
and Fig. 2(b).

The non-Hermitian bilayer model has two types of $\mathbb{Z}_{2}$ flux
operators, $\hat{W}_{p}$ and
$\hat{Z}_{p}=\tau_{1}^{z}\tau_{2}^{y}\tau_{3}^{x}\tau_{4}^{z}\tau_{5}^{y}\tau_{6}^{x},$
commuting with the Hamiltonian
($[\mathpzc{H},\hat{W}_{p}]=[\mathpzc{H},\hat{Z}_{p}]=[\hat{W}_{p},\hat{Z}_{p^{\prime}}]=0$).
This model can be viewed as a non-Hermitian analog of the bilayer spin model
in Ref. [47].

Figure 4: Decomposition of the density matrix. Numerically obtained density
matrices at different times ($t=0.1,0.3,0.5,0.7$) are decomposed by
$\rho_{n}=\langle\Psi_{n}|\hat{\rho}(t)|\Psi_{n}\rangle$. At $t=0.7$, the
components become almost constant, $\rho_{n}=2^{-N/2+1}\simeq 4.9\times
10^{-4}$ ($N=24$). The results are obtained with the coupling constants,
$K=-1$ and $\gamma=1$.

###  3.1 Steady state degeneracy

The bilayer model approach enables us to discover exact steady state solutions
and their extensive degeneracy. We find the steady state condition,

$\sigma_{j}^{\alpha}\sigma_{k}^{\alpha}\tau_{j}^{\alpha}\tau_{k}^{\alpha}|\rho\rangle\\!\rangle=|\rho\rangle\\!\rangle\leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\
\Longleftrightarrow\leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \mathpzc{H}|\rho\rangle\\!\rangle=0.$ (14)

Any state that satisfies the above local constraint is basically a steady
state. For instance, we may consider the singlet product state,

$|\lambda\rangle\\!\rangle=\otimes_{j}|s\rangle_{j},$ (15)

where
$|s\rangle=\frac{1}{\sqrt{2}}(|\uparrow\rangle\otimes|\downarrow\rangle-|\downarrow\rangle\otimes|\uparrow\rangle)$
is the spin-singlet of $\sigma$ and $\tau$ spins. One can check that
$\sigma_{j}^{\alpha}\sigma_{k}^{\alpha}\tau_{j}^{\alpha}\tau_{k}^{\alpha}|\lambda\rangle\\!\rangle=|\lambda\rangle\\!\rangle$.
From the singlet product state, we may generate further steady states as
follows.

$|\rho_{\rm
ss}\\{n\\}\rangle\\!\rangle=\prod_{l=x,y}(\hat{\mathcal{W}}_{l})^{n_{l}}\prod_{p}(\hat{W}_{p})^{n_{p}}|\lambda\rangle\\!\rangle,$
(16)

where $n_{l}=0,1\leavevmode\nobreak\ (l={x,y})$, and $n_{p}=0,1$ at each
plaquette $p$. Instead of the states, we consider the $\mathbb{Z}_{2}$ flux
eigenstates,

$|\rho_{\rm
ss}\\{w\\}\rangle\\!\rangle=\mathcal{N}\prod_{l=x,y}(1+w_{l}\hat{\mathcal{W}}_{l})\prod_{p}(1+w_{p}\hat{W}_{p})|\lambda\rangle\\!\rangle,$
(17)

where $\mathcal{N}$ is a normalization constant, and $w_{p}(=\pm 1)$ &
$w_{l}(=\pm 1)$ specify the flux sector including the Wilson loop flux. By
counting all these states, we find the degeneracy of the steady state
manifold, or steady state degeneracy (SSD):

${\rm SSD}=2^{2}\times 2^{N/2-1}.$ (18)

The first factor accounts for the four Wilson loop flux sectors ($w_{x}=\pm
1\leavevmode\nobreak\ \&\leavevmode\nobreak\ w_{y}=\pm 1$) while the second
one counts the number of possible flux sectors ($w_{p}=\pm 1$) on torus
geometry.

By applying the inverse vectorization ($|s\rangle\Rightarrow
i\sigma^{y}/\sqrt{2}$) to $|\rho_{\rm ss}\\{w\\}\rangle\\!\rangle$, we obtain
the density matrices of the steady states:

$\displaystyle\hat{\rho}_{\rm
ss}\\{w\\}=\mathcal{N}\prod_{l=x,y}(1+w_{l}\hat{\mathcal{W}}_{l})\prod_{p}(1+w_{p}\hat{W}_{p})\prod_{j}\frac{i\sigma_{j}^{y}}{\sqrt{2}}.\leavevmode\nobreak\
\leavevmode\nobreak\ $ (19)

On torus geometry, $\prod_{j}\frac{i\sigma_{j}^{y}}{\sqrt{2}}$ is equivalent
to a product of flux operators thus can be absorbed into the term
$\prod_{p}(1+w_{p}\hat{W}_{p})$; see Appendix B. If we focus on the flux
sector ($w_{p}=+1,w_{l}=-1$), we obtain the steady state in Eq. (8).

Figure 5: Two equivalent representations of the $\mathbb{Z}_{2}$ spin liquid
steady state $|\rho(t\rightarrow\infty)\rangle\\!\rangle$. (a) Toric code
state representation. Black lines denote the $\mathbb{Z}_{2}$ variables
$\\{X_{jk}=-1\\}$. The sign-flipped $\mathbb{Z}_{2}$ variables
$\\{X_{jk}=+1\\}$ are depicted by red lines. (b) Anyon condensation
representation. The condensed anyons are fermion pairs across the layers (blue
balls connected by a dashed line).

###  3.2 Dynamical anyon condensation

A fascinating physics is hidden in the Lindblad dynamics of Fig. 2, which we
unveil by using the bilayer description. The initial state is the KSL bilayer
product state, $|\rho(t=0)\rangle\\!\rangle=|\Psi_{\rm
KSL}\rangle\otimes|\Psi_{\rm KSL}\rangle$, and the steady state is given by
$|\rho(t\rightarrow\infty)\rangle\\!\rangle=|\rho_{\rm
MSKSL}\rangle\\!\rangle$, i.e.,

$|\rho(t\rightarrow\infty)\rangle\\!\rangle=\frac{1}{2^{N/2+1}}\prod_{l=x,y}(1-\hat{\mathcal{W}}_{l})\prod_{p}(1+\hat{W}_{p})|\lambda\rangle\\!\rangle.$

We find that this steady state is nothing but a $\mathbb{Z}_{2}$ spin liquid.
Therefore, the time evolution from $|\rho(t=0)\rangle\\!\rangle$ to
$|\rho(t\rightarrow\infty)\rangle\\!\rangle$ is a dynamical transition from
the KSL$\times$KSL state to a $\mathbb{Z}_{2}$ spin liquid state.

How does the state $|\rho(t\rightarrow\infty)\rangle\\!\rangle$ represent a
$\mathbb{Z}_{2}$ spin liquid? First, we note that $|\lambda\rangle\\!\rangle$
satisfies the condition $\phi_{j}^{x}=\phi_{j}^{y}=\phi_{j}^{z}=-1$ at every
site, where $\phi_{j}^{\alpha}\equiv\sigma_{j}^{\alpha}\tau_{j}^{\alpha}$
[47]. If a flux operator $\hat{W}_{p}$ acts on $|\lambda\rangle\\!\rangle$, it
flips the sign of $\phi_{j}^{\alpha}$ ($-1\rightarrow+1$) along the bonds of
the plaquette $p$. Using these properties, we assign a $\mathbb{Z}_{2}$
variable $X_{jk}(\equiv\phi_{j}^{\alpha}=\phi_{k}^{\alpha}=\pm 1)$ to each
bond $\langle jk\rangle_{\alpha}$. Then, we immediately recognize that the
state $\prod_{p}(1+\hat{W}_{p})|\lambda\rangle\\!\rangle$ is equivalent with
the $\mathbb{Z}_{2}$ spin liquid state of the toric code model [2]; see Fig.
5(a).

As already shown in Eqs. (8) and (9), the steady state can be represented by

$|\rho(t\rightarrow\infty)\rangle\\!\rangle=\sum_{n}\rho_{n}|\Psi_{n}\rangle\otimes|\Psi_{n}\rangle,$

where $\\{|\Psi_{n}\rangle\\}$ are the states of the zero flux sector
including the KSL ground state and excited states. Note that each state has an
even number of fermion excitations (without any flux excitations). The state
$|\rho(t\rightarrow\infty)\rangle\\!\rangle$ is an equal weight superposition
of all the bilayer product states
$\\{|\Psi_{n}\rangle\otimes|\Psi_{n}\rangle\\}$. This implies that
$|\rho(t\rightarrow\infty)\rangle\\!\rangle$ is the $\mathbb{Z}_{2}$ spin
liquid state that emerges from the KSL$\times$KSL state by condensing fermion
pairs between the two layers [Fig. 5(b)].

Figure 6: Dynamical anyon condensation. Red points denote the expectation
value $\langle\\!\langle\hat{\mathcal{A}}\rangle\\!\rangle$ as a function of
time $t$. The flux operator expectation value
$\langle\\!\langle\hat{W}\rangle\\!\rangle$ is plotted together for comparison
(blue). In the calculations, the coupling constants are fixed by $K=-1$ and
$\gamma=1$.

To confirm the fermion pair condensation, we consider the loop operator

$\hat{\mathcal{A}}=(\tau_{1}^{y}\tau_{6}^{y})(\tau_{6}^{z}\tau_{5}^{z})(\tau_{5}^{x}\tau_{4}^{x})(\sigma_{4}^{y}\sigma_{3}^{y})(\sigma_{3}^{z}\sigma_{2}^{z})(\sigma_{2}^{x}\sigma_{1}^{x}),$
(20)

which measures whether fermion excitations of the KSL can move between the two
layers, i.e., the fermion pair condensation between the layers [47]. We indeed
find that the expectation value
$\langle\\!\langle\hat{\mathcal{A}}\rangle\\!\rangle\equiv\langle\\!\langle{\rho}|\hat{\mathcal{A}}|{\rho}\rangle\\!\rangle/\langle\\!\langle{\rho}|{\rho}\rangle\\!\rangle$
increases from almost zero ($t=0$) to one ($t\rightarrow\infty$); see Fig. 6.
Furthermore, we explicitly check the steady state condition,
$\langle\\!\langle\sigma_{j}^{\alpha}\sigma_{k}^{\alpha}\tau_{j}^{\alpha}\tau_{k}^{\alpha}\rangle\\!\rangle=1$,
and entanglement entropy between the two layers in the steady state limit
(Appendix C). All these results confirm the dynamical anyon condensation from
the KSL$\times$KSL state to the $\mathbb{Z}_{2}$ spin liquid state.

##  4 Discussion

This work establishes that the $\mathbb{Z}_{2}$ toric code (TC) spin liquid is
equivalent with an anyon condensed phase of the Kitaev spin liquid bilayer.
Here, the condensed anyons are the fermion pairs $\psi_{1}\boxtimes\psi_{2}$
($\psi_{1},\psi_{2}$: fermion excitations in the upper and lower layer KSLs,
respectively).

$\mathcal{L}_{\rm KSL}:{\rm KSL}\times{\rm
KSL}\xrightarrow[]{\langle\psi_{1}\boxtimes\psi_{2}\rangle\neq
0}\mathbb{Z}_{2}{\rm TC}$ (21)

We showed this by using the Lindblad dynamics of the KSL and analyzing the
exact solution of the mixed-state KSL via the Choi-Jamiołkowski isomorphism.
The results summarized in Eq. (8) and Fig. 5 elucidate the anyon condensation
transition in a simple, exact, but nontrivial way.

###  4.1 Mechanism of the dynamical anyon condensation

The open quantum system of the Kitaev spin liquid reveals a general mechanism
for inducing anyon condensations in quantum spin liquids by the Lindblad
dynamics. A quantum spin liquid is usually realized by preserving a certain
type of local constraints (e.g., flux quantum number in the Kitaev spin
liquid). The quantum spin liquid is embedded in a Lindbladian that preserves
the local constraints. Then, the steady state should appear in the form of a
maximally mixed state within the subspace defined by the constraints. In the
bilayer description, the maximally mixed state corresponds to a new quantum
spin liquid emerging from the initial-state spin liquid bilayer by anyon
condensation. The key point is that the maximal mixing in the original system
corresponds to the anyon condensation in the doubled system.

###  4.2 Toric code open quantum systems

The mechanism applies to general quantum spin liquids. We take one more
example: an open quantum system of the toric code model [2]. The toric code
system $H_{\rm TC}=-J_{A}\sum_{s}A_{s}-J_{B}\sum_{p}B_{p}$ (where $J_{A,B}>0$,
$A_{s}=\prod_{j\in s}\sigma_{j}^{x}$ and $B_{p}=\prod_{j\in p}\sigma_{j}^{z}$)
is embedded in the Lindbladian with the jump operator
$L_{\mu}=\sqrt{\gamma}\sigma_{j}^{x}$. Note that the jump operator creates a
pair of $m$-particles ($B_{p}=-1$) without touching the $e$-particle sector.
In this case, the Lindblad equation [Eq. (1)] is exactly solved in the steady
state limit. The $\mathbb{Z}_{2}$ toric code spin liquid state
$|A_{s}=1\leavevmode\nobreak\ \&\leavevmode\nobreak\ B_{p}=1\rangle$ evolves
in time to the exact steady state

$\hat{\rho}_{\textup{ss}}\propto\prod_{s}\frac{1+A_{s}}{2},$ (22)

which is the maximally mixed state of all possible $m$-particle excitations
(with no $e$-particles). In the doubled system description, the maximally
mixedness corresponds to the condensation of $m$-particle pairs,
$m_{1}\boxtimes m_{2}$, yielding another $\mathbb{Z}_{2}$ toric code spin
liquid in the steady state $|\rho_{\textup{ss}}\rangle\\!\rangle$.

$\mathcal{L}_{\rm TC}:\mathbb{Z}_{2}{\rm TC}\times\mathbb{Z}_{2}{\rm
TC}\xrightarrow[]{\langle m_{1}\boxtimes m_{2}\rangle\neq 0}\mathbb{Z}_{2}{\rm
TC}$ (23)

One can condense $e$-particles ($A_{s}=-1$) instead by considering the jump
operator $L_{\mu}=\sqrt{\gamma}\sigma_{j}^{z}$. Qualitatively same results are
obtained in this case due to the duality between the $e$\- and $m$-particles.
See Appendix D for more details.

###  4.3 Summary

In this work, we proposed a new platform for anyon condensation topological
transitions, i.e., open quantum spin liquids. The main idea is to induce a
mixed-state quantum spin liquid in the steady state limit by using the
Lindblad dynamics, which gives rise to an anyon-condensed phase in the
corresponding doubled system. This framework can be applicable to more generic
systems with topological orders.

## Acknowledgements

I thank Jong Yeon Lee for valuable discussions and also for introducing
related studies [54, 55]. This work was supported by Individual Grants (No.
PG071402 & PG071403) of Korea Institute for Advanced Study (KIAS).
Computations were performed on clusters at the Center for Advanced Computation
(CAC) of KIAS.

Figure 7: The 24-site cluster and the Wilson loop operators

##  Appendix A Wilson loop operators and the reduced Hilbert space

To reduce the size of the Hilbert space in the actual calculations, we
consider two different Wilson loop flux operators $\hat{\mathcal{W}}_{x,y}$
(in addition to the hexagon flux operators $\\{\hat{W}_{p}\\}$). The Wilson
loop flux operators are defined by generalizing the hexagon flux operators to
non-contractible loops of the system. The 24-site cluster used in our
calculations is shown in Fig. 7(a); a periodic boundary condition is imposed
across the gray hexagon. Along the non-contractible loop in Fig. 7(b),
$\hat{\mathcal{W}}_{x}$ is defined by

$\displaystyle\hat{\mathcal{W}}_{x}$ $\displaystyle=$
$\displaystyle(\sigma_{12}^{z}\sigma_{5}^{z})(\sigma_{5}^{x}\sigma_{6}^{x})(\sigma_{6}^{y}\sigma_{7}^{y})(\sigma_{7}^{x}\sigma_{8}^{x})$
(24)
$\displaystyle(\sigma_{8}^{z}\sigma_{2}^{z})(\sigma_{2}^{x}\sigma_{3}^{x})(\sigma_{3}^{y}\sigma_{4}^{y})(\sigma_{4}^{x}\sigma_{12}^{x})$
$\displaystyle=$
$\displaystyle-\sigma_{12}^{y}\sigma_{5}^{y}\sigma_{6}^{z}\sigma_{7}^{z}\sigma_{8}^{y}\sigma_{2}^{y}\sigma_{3}^{z}\sigma_{4}^{z}.$
(25)

This is nothing but the product of the Kitaev terms along the non-contractible
loop. Repeating for other possible non-contractible loops [Fig. 7(c),(d)], we
can define other Wilson loop flux operators:

$\displaystyle\hat{\mathcal{W}}_{y}=-\sigma_{0}^{z}\sigma_{1}^{z}\sigma_{2}^{x}\sigma_{8}^{x}\sigma_{9}^{z}\sigma_{10}^{z}\sigma_{11}^{x}\sigma_{18}^{x},$
(26)
$\displaystyle\hat{\mathcal{W}}_{z}=-\sigma_{1}^{y}\sigma_{0}^{y}\sigma_{6}^{x}\sigma_{7}^{x}\sigma_{14}^{y}\sigma_{13}^{y}\sigma_{19}^{x}\sigma_{20}^{x}.$
(27)

All these Wilson loop flux operators commute with the Lindbladian:
$[H,\hat{\mathcal{W}}_{l}]=[L_{\mu},\hat{\mathcal{W}}_{l}]=0$ ($l=x,y,z$).
Note that the three operators are not independent of each other in the zero
flux sector ($W_{p}=+1$) due to the relationship
$\hat{\mathcal{W}}_{x}\hat{\mathcal{W}}_{y}\hat{\mathcal{W}}_{z}=+1$. Among
the three operators, we use $\hat{\mathcal{W}}_{x}$ and
$\hat{\mathcal{W}}_{y}$. Thanks to the flux quantum numbers of
$\hat{\mathcal{W}}_{x,y}$ and $\\{\hat{W}_{p}\\}$, the size of the Hilbert
space in each flux sector is significantly reduced to

$\frac{2^{N}}{2^{2}\cdot 2^{N/2-1}}=2^{N/2-1}.$ (28)

By solving the Hamiltonian $H$ by exact diagonalization, we obtain threefold
degenerate ground states $\\{|\psi_{\rm KSL}^{\rm I}\rangle,|\psi_{\rm
KSL}^{\rm II}\rangle,|\psi_{\rm KSL}^{\rm III}\rangle\\}$ which are
distinguished by different eigenvalues of the Wilson loop flux operators.

$\begin{array}[]{c|ccc|c}\hline\cr&\mathcal{W}_{x}&\mathcal{W}_{y}&\mathcal{W}_{z}&W_{p}\\\
\hline\cr|\psi_{\rm KSL}^{\rm I}\rangle&+1&-1&-1&+1\\\ |\psi_{\rm KSL}^{\rm
II}\rangle&-1&+1&-1&+1\\\ |\psi_{\rm KSL}^{\rm III}\rangle&-1&-1&+1&+1\\\
\hline\cr\end{array}$ (29)

In our calculations, the state $|\psi_{\rm KSL}^{\rm III}\rangle$ is chosen
for the initial state of the Lindblad dynamics [Eq. (5)].

##  Appendix B Exact steady state solution

Here we show the relationship,

$\prod_{j}\frac{i\sigma_{j}^{y}}{\sqrt{2}}\propto\prod_{p\in\mathcal{P}}\hat{W}_{p},$
(30)

where $\mathcal{P}$ means the set of the hexagon plaquettes marked in Fig.
7(a). The hexagon plaquettes cover a half of the cluster in the stripe pattern
perpendicular to the $y$-bond direction. In this case, one can find that

$\prod_{p\in\mathcal{P}}\hat{W}_{p}\propto\prod_{\langle
jk\rangle_{x}}\sigma_{j}^{x}\sigma_{k}^{x}\prod_{\langle
jk\rangle_{z}}\sigma_{j}^{z}\sigma_{k}^{z}\propto\prod_{j}i\sigma_{j}^{y}.$
(31)

From this, we identify the relationship, Eq. (30). If we apply this to Eq.
(19) together with $\mathcal{W}_{x,y}=-1$, $W_{p}=+1$ and appropriate
normalization, we finally obtain the solution in Eq. (8).

The relationship of Eq. (30) is valid on generic clusters if the periodic
boundary condition allows covering the hexagon plaquettes in the stripy
pattern as shown in Fig. 7(a).

Figure 8: Numerical results of the non-Hermitian Kitaev bilayer. (a) The
steady state condition
$\sigma_{j}^{\alpha}\sigma_{k}^{\alpha}\tau_{j}^{\alpha}\tau_{k}^{\alpha}=1$.
(b) The entanglement entropy $S_{\rm layer}$.

##  Appendix C Non-Hermitian Kitaev bilayer

To get convinced about the nature of the steady state, we have examined other
quantities of the non-Hermitian Kitaev bilayer. We checked the steady state
condition [Eq. (14)] by computing the expectation value,

$\langle\\!\langle\sigma_{j}^{\alpha}\sigma_{k}^{\alpha}\tau_{j}^{\alpha}\tau_{k}^{\alpha}\rangle\\!\rangle=\frac{\langle\\!\langle{\rho}|\sigma_{j}^{\alpha}\sigma_{k}^{\alpha}\tau_{j}^{\alpha}\tau_{k}^{\alpha}|{\rho}\rangle\\!\rangle}{\langle\\!\langle{\rho}|{\rho}\rangle\\!\rangle}.$
(32)

Figure 8(a) shows the result. In the initial Kitaev bilayer product state
($t=0$), we have
$\langle\\!\langle\sigma_{j}^{\alpha}\sigma_{k}^{\alpha}\tau_{j}^{\alpha}\tau_{k}^{\alpha}\rangle\\!\rangle=\langle\\!\langle\sigma_{j}^{\alpha}\sigma_{k}^{\alpha}\rangle\\!\rangle\langle\\!\langle\tau_{j}^{\alpha}\tau_{k}^{\alpha}\rangle\\!\rangle$.
This quantity gradually increases and converges to the expected value
$\langle\\!\langle\sigma_{j}^{\alpha}\sigma_{k}^{\alpha}\tau_{j}^{\alpha}\tau_{k}^{\alpha}\rangle\\!\rangle=1$,
i.e.,
$\langle\\!\langle\phi_{j}^{\alpha}\phi_{k}^{\alpha}\rangle\\!\rangle=1$. This
confirms that the steady state has the property
$\phi_{j}^{\alpha}=\phi_{k}^{\alpha}(=\pm 1)$ at each bond $\langle
jk\rangle_{\alpha}$, allowing us to define the $\mathbb{Z}_{2}$ variable
$X_{jk}(\equiv\phi_{j}^{\alpha}=\phi_{k}^{\alpha})$ to characterize the steady
state.

We also calculated the entanglement entropy between the two layers,

$S_{\rm layer}=-\log({\rm Tr}_{\\{\sigma\\}}\rho_{\rm layer}^{2}),$ (33)

where $\rho_{\rm layer}={\rm
Tr}_{\\{\tau\\}}|\tilde{\rho}\rangle\\!\rangle\langle\\!\langle\tilde{\rho}|$
is the reduced density matrix obtained by tracing out the $\tau$-spins in the
normalized wave function
$|\tilde{\rho}\rangle\\!\rangle=|\rho\rangle\\!\rangle/\sqrt{\langle\\!\langle\rho|\rho\rangle\\!\rangle}$.
Figure 8(b) shows the calculated entanglement entropy. Starting from zero in
the initial Kitaev bilayer product state, the entanglement entropy reaches the
maximum value $S_{\rm layer}=2^{N/2-1}$ in the steady state. The steady state
entanglement entropy is in agreement with the number of all possible states
satisfying the steady state condition $\phi_{j}^{\alpha}\phi_{k}^{\alpha}=1$
[Eq. (14)].

##  Appendix D Toric code open quantum systems

The open quantum system of the toric code model can be exactly solved in the
steady state limit. We find the exact steady state solution by taking the
mapping to the non-Hermitian bilayer system,

$\displaystyle\mathpzc{H}_{\rm TC}$ $\displaystyle=$ $\displaystyle-
J_{A}\sum_{s}A_{s}^{(\sigma)}-J_{B}\sum_{p}B_{p}^{(\sigma)}$ (34)
$\displaystyle+J_{A}\sum_{s}A_{s}^{(\tau)}+J_{B}\sum_{p}B_{p}^{(\tau)}$
$\displaystyle+i\gamma\sum_{j}(\sigma_{j}^{x}\tau_{j}^{x}-1),$

where the upper layer and lower layer terms are denoted by the superscripts
$(\sigma,\tau)$, respectively, and

$\displaystyle A_{s}^{(\sigma)}=\prod_{j\in
s}\sigma_{j}^{x},\leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ B_{p}^{(\sigma)}=\prod_{j\in p}\sigma_{j}^{z},$ (35)
$\displaystyle A_{s}^{(\tau)}=\prod_{j\in s}\tau_{j}^{x},\leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ B_{p}^{(\tau)}=\prod_{j\in
p}\tau_{j}^{z}.$ (36)

We are going to use the basis shown in Table 1. In this basis, spin operators
are given by

$\displaystyle\sigma^{x}$ $\displaystyle=$ $\displaystyle|s\rangle\langle
t_{x}|+|t_{x}\rangle\langle s|-i|t_{y}\rangle\langle
t_{z}|+i|t_{z}\rangle\langle t_{y}|,$ (37) $\displaystyle\tau^{x}$
$\displaystyle=$ $\displaystyle-|s\rangle\langle t_{x}|-|t_{x}\rangle\langle
s|-i|t_{y}\rangle\langle t_{z}|+i|t_{z}\rangle\langle t_{y}|.$ (38)

Other spin operators are obtained by cyclic permutations of the above.

$\sigma^{z}\otimes\tau^{z}$ basis |  $|\phi^{x},\phi^{y},\phi^{z}\rangle$ basis
---|---
$|s\rangle=\frac{1}{\sqrt{2}}(|\uparrow\downarrow\rangle-|\downarrow\uparrow\rangle)$ | $|-1,-1,-1\rangle$
$|t_{x}\rangle=\frac{-1}{\sqrt{2}}(|\uparrow\uparrow\rangle-|\downarrow\downarrow\rangle)$ | $|-1,+1,+1\rangle$
$|t_{y}\rangle=\frac{i}{\sqrt{2}}(|\uparrow\uparrow\rangle+|\downarrow\downarrow\rangle)$ | $|+1,-1,+1\rangle$
$|t_{z}\rangle=\frac{1}{\sqrt{2}}(|\uparrow\downarrow\rangle+|\downarrow\uparrow\rangle)$ | $|+1,+1,-1\rangle$
Table 1: Four local states of $\sigma$ and $\tau$ spins. The four states
$\\{|s\rangle,|t_{x}\rangle,|t_{y}\rangle,|t_{z}\rangle\\}$ are eigenstates of
$\phi^{x,y,z}$ (where $\phi^{\alpha}=\sigma^{\alpha}\tau^{\alpha}$).

The steady state solution of the non-Hermitian Hamiltonian $\mathpzc{H}_{\rm
TC}$ takes the following form:

$|\rho_{\rm
ss}\rangle\\!\rangle=\mathcal{N}\left(\prod_{s}\frac{1+A_{s}^{(\sigma)}}{2}\right)\otimes_{l}|t_{y}\rangle_{l},$
(39)

where $\mathcal{N}$ is a normalization constant. Now we show that $|\rho_{\rm
ss}\rangle\\!\rangle$ is indeed the steady state ($\mathpzc{H}_{\rm
TC}|\rho_{\rm ss}\rangle\\!\rangle=0$). First, we note the property

$\sigma_{j}^{x}\tau_{j}^{x}|\rho_{\rm ss}\rangle\\!\rangle=|\rho_{\rm
ss}\rangle\\!\rangle,$ (40)

which is easily checked from the identity
$\sigma^{x}\tau^{x}|t_{y}\rangle=\phi^{x}|t_{y}\rangle=|t_{y}\rangle$. From
the above property, we find that

$A_{s}^{(\sigma)}|\rho_{\rm ss}\rangle\\!\rangle=A_{s}^{(\tau)}|\rho_{\rm
ss}\rangle\\!\rangle=|\rho_{\rm ss}\rangle\\!\rangle.$ (41)

Lastly, we note that

$B_{p}^{(\sigma)}|\rho_{\rm ss}\rangle\\!\rangle=B_{p}^{(\tau)}|\rho_{\rm
ss}\rangle\\!\rangle,$ (42)

which is derived from the identity
$\sigma^{z}|t_{y}\rangle=-i|t_{x}\rangle=\tau^{z}|t_{y}\rangle$. The three
properties of Eqs. (40),(41),(42) prove that $\mathpzc{H}_{\rm TC}|\rho_{\rm
ss}\rangle\\!\rangle=0$.

We go back to the original system by applying the inverse vectorization
$|t_{y}\rangle\Rightarrow i\mathbb{I}_{2}/\sqrt{2}$ ($\mathbb{I}_{2}$:
2$\times$2 identity matrix) to the state $|\rho_{\rm ss}\rangle\\!\rangle$.
Then, we obtain the steady state density matrix shown in Eq. (22):

$\hat{\rho}_{\rm
ss}=\mathcal{N}\left(\frac{i}{\sqrt{2}}\right)^{N}\left(\prod_{s}\frac{1+A_{s}}{2}\right),$
(43)

where $N$ is the number of sites. This is the projection operator into the
subspace of $\\{A_{s}=1\\}$, leading to another representation

$\hat{\rho}_{\rm
ss}\propto\sum_{n}|\Psi_{n}^{(m)}\rangle\langle\Psi_{n}^{(m)}|,$ (44)

i.e., the maximally mixed state of all possible $m$-particle excitations (with
no $e$-particles) with equal weight. This implies that the state $|\rho_{\rm
ss}\rangle\\!\rangle\propto\sum_{n}|\Psi_{n}^{(m)}\rangle\otimes|\Psi_{n}^{(m)}\rangle$
is an $m$-particle-condensed state of the $\mathbb{Z}_{2}\times\mathbb{Z}_{2}$
toric code bilayer.

On the other hand, the state $|\rho_{\rm ss}\rangle\\!\rangle$ can be
represented in the form of a $\mathbb{Z}_{2}$ toric code state. To show this,
we define a new set of Pauli matrices in the space of $|t_{y}\rangle$ and
$|t_{z}\rangle$.

$X\equiv|t_{y}\rangle\langle\tilde{t}_{z}|+|\tilde{t}_{z}\rangle\langle
t_{y}|\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\&\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
Z\equiv|t_{y}\rangle\langle
t_{y}|-|\tilde{t}_{z}\rangle\langle\tilde{t}_{z}|,$ (45)

where we used $|\tilde{t}_{z}\rangle=i|t_{z}\rangle$ instead of
$|t_{z}\rangle$. Using the new operators, we can rewrite the state as

$|\rho_{\rm
ss}\rangle\\!\rangle=\mathcal{N}\left(\prod_{s}\frac{1+A_{s}^{(X)}}{2}\right)\otimes_{l}|Z=1\rangle_{l}.$
(46)

We may also write it as

$|\rho_{\rm
ss}\rangle\\!\rangle=\mathcal{N}\left(\prod_{p}\frac{1+B_{p}^{(Z)}}{2}\right)\otimes_{l}|X=1\rangle_{l}.$
(47)

This is exactly the $\mathbb{Z}_{2}$ toric code state of
$\\{A_{s}=1\leavevmode\nobreak\ \&\leavevmode\nobreak\ B_{p}=1\\}$. All these
results confirm the anyon condensation transition summarized in Eq. (23).

## References

  * [1] A. Kitaev.  “Anyons in an exactly solved model and beyond”.  Ann. Phys. 321, 2–111 (2006).
  * [2] A.Yu. Kitaev.  “Fault-tolerant quantum computation by anyons”.  Ann. Phys. 303, 2–30 (2003).
  * [3] X.-G. Wen.  “Quantum Field Theory of Many-Body Systems: From the Origin of Sound to an Origin of Light and Electrons”.  Oxford University Press. (2007).
  * [4] S. Sachdev.  “Quantum Phases of Matter”.  Cambridge University Press. (2023).
  * [5] L. Savary and L. Balents.  “Quantum spin liquids: a review”.  Rep. Prog. Phys. 80, 016502 (2016).
  * [6] C. Broholm, R. J. Cava, S. A. Kivelson, D. G. Nocera, M. R. Norman, and T. Senthil.  “Quantum spin liquids”.  Science 367, eaay0668 (2020).
  * [7] J. Knolle and R. Moessner.  “A Field Guide to Spin Liquids”.  Annu. Rev. Condens. Matter Phys. 10, 451–472 (2019).
  * [8] Y. Kasahara, T. Ohnishi, Y. Mizukami, O. Tanaka, S. Ma, K. Sugii, N. Kurita, H. Tanaka, J. Nasu, Y. Motome, et al.  “Majorana quantization and half-integer thermal quantum Hall effect in a Kitaev spin liquid”.  Nature 559, 227–231 (2018).
  * [9] H. Takagi, T. Takayama, G. Jackeli, G. Khaliullin, and S. E Nagler.  “Concept and realization of Kitaev quantum spin liquids”.  Nat. Rev. Phys. 1, 264–280 (2019).
  * [10] Y. Motome and J. Nasu.  “Hunting Majorana Fermions in Kitaev Magnets”.  J. Phys. Soc. Jpn. 89, 012002 (2019).
  * [11] S. Trebst and C. Hickey.  “Kitaev materials”.  Physics Reports 950, 1–37 (2022).
  * [12] E. Dennis, A. Kitaev, A. Landahl, and J. Preskill.  “Topological quantum memory”.  Journal of Mathematical Physics 43, 4452–4505 (2002).
  * [13] H. Bombin and M. A. Martin-Delgado.  “Family of non-Abelian Kitaev models on a lattice: Topological condensation and confinement”.  Phys. Rev. B 78, 115421 (2008).
  * [14] H. Bombin.  “Topological Order with a Twist: Ising Anyons from an Abelian Model”.  Phys. Rev. Lett. 105, 030403 (2010).
  * [15] J. K. Pachos.  “Introduction to Topological Quantum Computation”.  Cambridge University Press. Cambridge (2012).
  * [16] B. M. Terhal.  “Quantum error correction for quantum memories”.  Rev. Mod. Phys. 87, 307–346 (2015).
  * [17] H.-N. Dai, B. Yang, A. Reingruber, H. Sun, X.-F. Xu, Y.-A. Chen, Z.-S. Yuan, and J.-W. Pan.  “Four-body ring-exchange interactions and anyonic statistics within a minimal toric-code Hamiltonian”.  Nature Physics 13, 1195–1200 (2017).
  * [18] K. J. Satzinger, Y.-J. Liu, A. Smith, C. Knapp, M. Newman, C. Jones, Z. Chen, C. Quintana, X. Mi, A. Dunsworth, et al.  “Realizing topologically ordered states on a quantum processor”.  Science 374, 1237–1241 (2021).
  * [19] G. Semeghini, H. Levine, A. Keesling, S. Ebadi, T. T. Wang, D. Bluvstein, R. Verresen, H. Pichler, M. Kalinowski, R. Samajdar, et al.  “Probing topological spin liquids on a programmable quantum simulator”.  Science 374, 1242–1247 (2021).
  * [20] R. Samajdar, W. W. Ho, H. Pichler, M. D. Lukin, and S. Sachdev.  “Quantum phases of Rydberg atoms on a kagome lattice”.  Proc. Natl. Acad. Sci. U.S.A. 118, e2015785118 (2021).
  * [21] R. Verresen, M. D. Lukin, and A. Vishwanath.  “Prediction of Toric Code Topological Order from Rydberg Blockade”.  Phys. Rev. X 11, 031005 (2021).
  * [22] T. D. Ellison, Y.-A. Chen, A. Dua, W. Shirley, N. Tantivasadakarn, and D. J. Williamson.  “Pauli Stabilizer Models of Twisted Quantum Doubles”.  PRX Quantum 3, 010353 (2022).
  * [23] Y.-J. Liu, K. Shtengel, A. Smith, and F. Pollmann.  “Methods for Simulating String-Net States and Anyons on a Digital Quantum Computer”.  PRX Quantum 3, 040315 (2022).
  * [24] T.-C. Lu, L. A. Lessa, I. H. Kim, and T. H. Hsieh.  “Measurement as a Shortcut to Long-Range Entangled Quantum Matter”.  PRX Quantum 3, 040337 (2022).
  * [25] N. Tantivasadakarn, A. Vishwanath, and R. Verresen.  “Hierarchy of Topological Order From Finite-Depth Unitaries, Measurement, and Feedforward”.  PRX Quantum 4, 020339 (2023).
  * [26] Google Quantum AI and Collaborators.  “Non-Abelian braiding of graph vertices in a superconducting processor”.  Nature 618, 264–269 (2023).
  * [27] R. Verresen, N. Tantivasadakarn, and A. Vishwanath.  “Efficiently preparing Schrödinger’s cat, fractons and non-Abelian topological order in quantum devices” (2022).  arXiv:2112.03061.
  * [28] S. Bravyi, I. Kim, A. Kliesch, and R. Koenig.  “Adaptive constant-depth circuits for manipulating non-abelian anyons” (2022).  arXiv:2205.01933.
  * [29] N. Tantivasadakarn, R. Verresen, and A. Vishwanath.  “The shortest route to non-abelian topological order on a quantum processor” (2022).  arXiv:2209.03964.
  * [30] M. Foss-Feig, A. Tikku, T.-C. Lu, K. Mayer, M. Iqbal, T. M. Gatterman, J. A. Gerber, K. Gilmore, D. Gresh, A. Hankin, N. Hewitt, C. V. Horst, M. Matheny, T. Mengle, B. Neyenhuis, H. Dreyer, D. Hayes, T. H. Hsieh, and I. H. Kim.  “Experimental demonstration of the advantage of adaptive quantum circuits” (2023).  arXiv:2302.03029.
  * [31] M. Iqbal, N. Tantivasadakarn, R. Verresen, S. L. Campbell, J. M. Dreiling, C. Figgatt, J. P. Gaebler, J. Johansen, M. Mills, S. A. Moses, J. M. Pino, A. Ransford, M. Rowe, P. Siegfried, R. P. Stutz, M. Foss-Feig, A. Vishwanath, and H. Dreyer.  “Creation of Non-Abelian Topological Order and Anyons on a Trapped-Ion Processor” (2023).  arXiv:2305.03766.
  * [32] J. Preskill.  “Quantum Computing in the NISQ era and beyond”.  Quantum 2, 79 (2018).
  * [33] B. Kraus, H. P. Büchler, S. Diehl, A. Kantian, A. Micheli, and P. Zoller.  “Preparation of entangled states by quantum Markov processes”.  Phys. Rev. A 78, 042307 (2008).
  * [34] S. Diehl, A. Micheli, A. Kantian, B. Kraus, H. P. Büchler, and P. Zoller.  “Quantum states and phases in driven open quantum systems with cold atoms”.  Nat. Phys. 4, 878–883 (2008).
  * [35] S. Diehl, E. Rico, M. A. Baranov, and P. Zoller.  “Topology by dissipation in atomic quantum wires”.  Nat. Phys. 7, 971–977 (2011).
  * [36] G. Lindblad.  “On the generators of quantum dynamical semigroups”.  Communications in Mathematical Physics 48, 119–130 (1976).
  * [37] V. Gorini, A. Kossakowski, and E. C. G. Sudarshan.  “Completely positive dynamical semigroups of N-level systems”.  Journal of Mathematical Physics 17, 821–825 (1976).
  * [38] H.-P. Breuer and F. Petruccione.  “The Theory of Open Quantum Systems”.  Oxford University Press. (2007).
  * [39] V. V. Albert and L. Jiang.  “Symmetries and conserved quantities in Lindblad master equations”.  Phys. Rev. A 89, 022118 (2014).
  * [40] V. V. Albert, B. Bradlyn, M. Fraas, and L. Jiang.  “Geometry and Response of Lindbladians”.  Phys. Rev. X 6, 041031 (2016).
  * [41] D. A. Lidar.  “Lecture Notes on the Theory of Open Quantum Systems” (2020).  arXiv:1902.00967.
  * [42] P. M. Harrington, E. J. Mueller, and K. W. Murch.  “Engineered dissipation for quantum information science”.  Nat. Rev. Phys. 4, 660–671 (2022).
  * [43] F. A. Bais and J. K. Slingerland.  “Condensate-induced transitions between topologically ordered phases”.  Phys. Rev. B 79, 045316 (2009).
  * [44] F. J. Burnell, S. H. Simon, and J. K. Slingerland.  “Condensation of achiral simple currents in topological lattice models: Hamiltonian study of topological symmetry breaking”.  Phys. Rev. B 84, 125434 (2011).
  * [45] F. J. Burnell.  “Anyon Condensation and Its Applications”.  Annu. Rev. Condens. Matter Phys. 9, 307–327 (2018).
  * [46] R. Wiedmann, L. Lenke, M. R. Walther, M. Mühlhauser, and K. P. Schmidt.  “Quantum critical phase transition between two topologically ordered phases in the Ising toric code bilayer”.  Phys. Rev. B 102, 214422 (2020).
  * [47] K. Hwang.  “Topological Quantum Dimers Emerging from Kitaev Spin Liquid Bilayer: Anyon Condensation Transition” (2023).  arXiv:2301.05721.
  * [48] N. Shibata and H. Katsura.  “Dissipative spin chain as a non-Hermitian Kitaev ladder”.  Phys. Rev. B 99, 174303 (2019).
  * [49] N. Shibata and H. Katsura.  “Dissipative quantum Ising chain as a non-Hermitian Ashkin-Teller model”.  Phys. Rev. B 99, 224432 (2019).
  * [50] Y. Saad.  “Analysis of Some Krylov Subspace Approximations to the Matrix Exponential Operator”.  SIAM Journal on Numerical Analysis 29, 209–228 (1992).
  * [51] R. B. Sidje.  “Expokit: A software package for computing matrix exponentials”.  ACM Transactions on Mathematical Software (TOMS) 24, 130–156 (1998).
  * [52] A. Jamiołkowski.  “Linear transformations which preserve trace and positive semidefiniteness of operators”.  Reports on Mathematical Physics 3, 275–278 (1972).
  * [53] M.-D. Choi.  “Completely positive linear maps on complex matrices”.  Linear Algebra and its Applications 10, 285–290 (1975).
  * [54] Y. Bao, R. Fan, A. Vishwanath, and E. Altman.  “Mixed-state topological order and the errorfield double formulation of decoherence-induced transitions” (2023).  arXiv:2301.05687.
  * [55] J. Y. Lee, C.-M. Jian, and C. Xu.  “Quantum criticality under decoherence or weak measurement” (2023).  arXiv:2301.05238.