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

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 environment 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.

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".
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 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.

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 Lind- blad master equation [36][37][38][39][40][41][42], where the system's density matrix ρ(t) evolves in time by the Lindbladian L composed of the Hamiltonian H and the Lindblad operators {L µ } describing the environment effects.The Lindblad equation has the completely positive and trace-preserving property, i.e., Trρ(t) = const.For the Hamiltonian, we consider the Kitaev honeycomb model [1], where σ x,y,z mean the Pauli matrices and ⟨jk⟩ 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, with the dissipation strength γ (≥ 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 ⟨jk⟩ α .The Hamiltonian H and the Lindblad operators {L µ } both commute with the Z 2 flux operator, i.e., [H, Ŵp ] = [L µ , Ŵp ] = 0. Hence, the Kitaev Lindbladian L conserves the Z 2 flux quantum number.

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, where |Ψ KSL ⟩ 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 |Ψ KSL ⟩ in the sector of the Wilson loop flux W x,y = −1 in our calculations (see Appendix A).The time evolution of the density matrix, ρ(t) = e tL ρ(0), is computed by using the Krylov subspace methods [50,51].In the construction of the time evolution operator e tL , we utilize the Z 2 flux quantum numbers of the Kitaev Lindbladian to reduce the size of the Hilbert space.
Figure 2 shows the calculation results of the Lindblad dynamics.First, we check the properties Trρ(t) = 1 and ⟨ Ŵp ⟩ = Tr[ρ(t) Ŵp ] = 1 in Fig. 2(a).Next, we consider the Rényi entropy, 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 Rényi = log(2 N/2−1 ) when t → ∞ (N = 24 is the number of sites).The larger the dissipation strength γ is, the faster the convergence is reached.The constant value reveals the size of the subspace (2  when t → ∞ [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 spinspin correlation disappears in the long time limit.This peculiar property is also reflected in the spin structure factor, which becomes flat in momentum space in the long time limit: S(q) = 3 when t → ∞ (see Fig. 3).

Exact solution of the steady state
The numerical approach establishes the existence of the steady state characterized by the properties, ⟨ Ŵp ⟩ = 1, S Rényi = log(2 N/2−1 ), and ⟨σ α j σ α k ⟩ = 0. We find the exact solution for the steady state: The steady state density matrix is simply given by the projection operator into the zero flux sector with the Wilson loop flux W x,y = −1.Hence, it is diagonal in the eigenbasis {|Ψ n ⟩} of H and Ŵp with the weight This equal weight property is indeed identified in our numerical calculations shown in Fig. 4. One can check that L(ρ 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]: where the bra of the density matrix ρ is turned into a ket of the resulting state vector |ρ⟩⟩.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 with the non-Hermitian Hamiltonian In our case, the non-Hermitian Hamiltonian is given by where σ & τ are Pauli spin operators acting on the original kets and copied kets; We now have a bilayer spin model with the intralayer Kitaev interactions (K & −K) and the interlayer interactions (iγ) [Fig.1(b)].Note that the interlayer interactions are non-Hermitian originating from the environment effects.This leads to nonunitary time evolution of the state vector, |ρ(t)⟩⟩ = e −itH |ρ(0)⟩⟩, which is identified from the relationship ⟨⟨ρ|ρ⟩⟩ = exp(−S Rényi ) and Fig. 2(b).

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, Any state that satisfies the above local constraint is basically a steady state.For instance, we may consider the singlet product state, From the singlet product state, we may generate further steady states as follows.
where n l = 0, 1 (l = x, y), and n p = 0, 1 at each plaquette p. Instead of the states, we consider the Z 2 flux eigenstates, where N is a normalization constant, and w p (= ±1) & w l (= ±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): The first factor accounts for the four Wilson loop flux sectors (w x = ±1 & w y = ±1) while the second one counts the number of possible flux sectors (w p = ±1) on torus geometry.By applying the inverse vectorization (|s⟩ ⇒ iσ y / √ 2) to |ρ ss {w}⟩⟩, we obtain the density matrices of the steady states: On torus geometry, j iσ y j √ 2 is equivalent to a product of flux operators thus can be absorbed into the term p (1 + w p Ŵ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).

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, |ρ(t = 0)⟩⟩ = |Ψ KSL ⟩ ⊗ |Ψ KSL ⟩, and the steady state is given by |ρ(t → ∞)⟩⟩ = |ρ MSKSL ⟩⟩, i.e., We find that this steady state is nothing but a Z 2 spin liquid.Therefore, the time evolution from |ρ(t = 0)⟩⟩ to |ρ(t → ∞)⟩⟩ is a dynamical transition from the KSL×KSL state to a Z 2 spin liquid state.
As already shown in Eqs. ( 8) and ( 9), the steady state can be represented by where {|Ψ n ⟩} 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 |ρ(t → ∞)⟩⟩ is an equal weight superposition of all the bilayer product states {|Ψ n ⟩ ⊗ |Ψ n ⟩}.This implies that |ρ(t → ∞)⟩⟩ is the Z 2 spin liquid state that emerges from the KSL×KSL state by condensing fermion pairs between the two layers [Fig.5(b)].
To confirm the fermion pair condensation, we consider the loop operator 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 ⟨⟨ Â⟩⟩ ≡ ⟨⟨ρ| Â|ρ⟩⟩/⟨⟨ρ|ρ⟩⟩ increases from almost zero (t = 0) to one (t → ∞); see Fig. 6.Furthermore, we explicitly check the steady state condition, ⟨⟨σ α j σ α k τ α j τ α k ⟩⟩ = 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×KSL state to the Z 2 spin liquid state.

Discussion
This work establishes that the 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 ψ 1 ⊠ ψ 2 (ψ 1 , ψ 2 : fermion excitations in the upper and lower layer KSLs, respectively).
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.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.
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 TC = −J A s A s − J B p B p (where J A,B > 0, A s = j∈s σ x j and B p = j∈p σ z j ) is embedded in the Lindbladian with the jump operator L µ = √ γσ x j .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 Z 2 toric code spin liquid state |A s = 1 & B p = 1⟩ evolves in time to the exact steady state which is the maximally mixed state of all possible mparticle excitations (with no e-particles).In the doubled system description, the maximally mixedness corresponds to the condensation of m-particle pairs, m 1 ⊠m 2 , yielding another Z 2 toric code spin liquid in the steady state |ρ ss ⟩⟩.
One can condense e-particles (A s = −1) instead by considering the jump operator L µ = √ γσ z j .Qualitatively same results are obtained in this case due to the duality between the e-and m-particles.See Appendix D for more details.
We make some comments on related works to our study.Lindblad open quantum systems of the KSL have been studied in a different setup from ours in Ref. [54]; they assumed the jump operator L µ = σ α j + σ α k and no occurrence of quantum jump [accordingly, they ignored the term L µ ρL † µ in the master equation, Eq. ( 1)].
This work shows that the non-Hermitian physics without quantum jump leads to interesting spin liquids featured with exceptional points, fermi arcs, and non-Hermitian skin effects.Such exceptional spin liquids are only maintained in the absence of any quantum jump event.By contrast, our mixed-state spin liquids appear in the steady state limit as an outcome of quantum jump events.
Recently, decoherence-induced transitions and mixed-state topological orders are widely investigated in various aspects.The same anyon condensation physics in Eq. ( 23) was shown to arise in the toric code subject to local decoherence channels [55,56].For the intrinsic characterizations of mixed-state topological orders, several information-theoretical diagnostics have been proposed and examined such as quantum relative entropy, coherent information, and topological entanglement negativity [57,58].General frameworks to characterize/classify mixed-state topological orders are also being developed very recently [59,60].Our mixed-state spin liquids appearing from the Lindblad dynamics provide concrete realizations of mixed-state topological orders.Exact solutions of our mixed-state spin liquids are expected to help to sharpen our understanding on mixed-state topological orders.
To summarize, 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 anyoncondensed phase in the corresponding doubled system.This framework can be applicable to more generic systems with topological orders.

A Wilson loop operators and the reduced Hilbert space
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: All these Wilson loop flux operators commute with the Lindbladian: [H, Ŵl ] = [L µ , Ŵ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 Ŵx Ŵy Ŵz = +1.Among the three operators, we use Ŵx and Ŵy .Thanks to the flux quantum numbers of Ŵx,y and { Ŵp }, the size of the Hilbert space in each flux sector is significantly reduced to By solving the Hamiltonian H by exact diagonalization, we obtain threefold degenerate ground states {|ψ I KSL ⟩, |ψ II KSL ⟩, |ψ III KSL ⟩} which are distinguished by different eigenvalues of the Wilson loop flux operators.
In our calculations, the state |ψ III KSL ⟩ is chosen for the initial state of the Lindblad dynamics [Eq.( 5)].

B Exact steady state solution
Here we show the relationship, where 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 From this, we identify the relationship, Eq. ( 30).If we apply this to Eq. ( 19) together with 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).
We also calculated the entanglement entropy between the two layers, where ρ layer = Tr {τ } |ρ⟩⟩⟨⟨ρ| is the reduced density matrix obtained by tracing out the τ -spins in the normalized wave function |ρ⟩⟩ = |ρ⟩⟩/ ⟨⟨ρ|ρ⟩⟩.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 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 ϕ α j ϕ α k = 1 [Eq.( 14)].

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, where the upper layer and lower layer terms are denoted by the superscripts (σ, τ ), respectively, and We are going to use the basis shown in Table 1.In this basis, spin operators are given by  Other spin operators are obtained by cyclic permutations of the above.
The steady state solution of the non-Hermitian Hamiltonian H TC takes the following form: where N is a normalization constant.Now we show that |ρ ss ⟩⟩ is indeed the steady state (H TC |ρ ss ⟩⟩ = 0).First, we note the property which is easily checked from the identity σ x τ x |t y ⟩ = ϕ x |t y ⟩ = |t y ⟩.From the above property, we find that Lastly, we note that which is derived from the identity σ z |t y ⟩ = −i|t x ⟩ = τ z |t y ⟩.The three properties of Eqs. ( 40),( 41), (42) prove that H TC |ρ ss ⟩⟩ = 0.
We go back to the original system by applying the inverse vectorization |t y ⟩ ⇒ iI 2 / √ 2 (I 2 : 2×2 identity matrix) to the state |ρ ss ⟩⟩.Then, we obtain the steady state density matrix shown in Eq. (22): where N is the number of sites.This is the projection operator into the subspace of {A s = 1}, leading to another representation i.e., the maximally mixed state of all possible m-particle excitations (with no e-particles) with equal weight.This implies that the state |ρ ss ⟩⟩ ∝ n |Ψ We may also write it as This is exactly the Z 2 toric code state of {A s =

Figure 1 :
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µ = √ γσ α j σ α k .(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γ(σ α j σ α k τ α j τ α k − 1).The two descriptions are compared in the table.

Figure 2 :Figure 3 :
Figure 2: Lindblad dynamics of the Kitaev spin liquid.(a) The trace of the density matrix Tr ρ (t) and the expectation value of the flux operator ⟨ Ŵp⟩.(b) The Rényi entropy S Ré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 ⟨σ α j σ α k ⟩ for the nearest-neighbor bond ⟨jk⟩α.Different colors denote the results of different dissipation strengths (γ = 0, 0.2, 0.4, 0.6, 0.8, 1).The unitary evolution of the closed system (γ = 0) is shown together for comparison (black).In all the calculations, the coupling constant is fixed by K = −1.

Figure 6 :
Figure 6: Dynamical anyon condensation.Red points denote the expectation value ⟨⟨ Â⟩⟩ as a function of time t.The flux operator expectation value ⟨⟨ Ŵ ⟩⟩ is plotted together for comparison (blue).In the calculations, the coupling constants are fixed by K = −1 and γ = 1.
To reduce the size of the Hilbert space in the actual calculations, we consider two different Wilson loop flux operators Ŵx,y (in addition to the hexagon flux operators { Ŵp }).The Wilson loop flux operators are defined by generalizing the hexagon flux operators to non-

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

2 ⊗
n ⟩ is an m-particle-condensed state of the Z 2 × Z 2 toric code bilayer.On the other hand, the state |ρ ss ⟩⟩ can be represented in the form of a Z 2 toric code state.To show this, we define a new set of Pauli matrices in the space of |t y ⟩ and |t z ⟩.X ≡ |t y ⟩⟨ tz | + | tz ⟩⟨t y | & Z ≡ |t y ⟩⟨t y | − | tz ⟩⟨ tz |, (45) where we used | tz ⟩ = i|t z ⟩ instead of |t z ⟩.Using the new operators, we can rewrite the state as |ρ ss ⟩⟩ = N s l |Z = 1⟩ l .(