Anyon condensation and confinement transition in a Kitaev spin liquid bilayer

Transitions between quantum spin liquids (QSLs) are fundamental problems lying beyond the Landau paradigm and requiring a deep understanding of the entanglement structures of QSLs called topological orders. The novel concept of anyon condensation has been proposed as a theoretical mechanism, predicting various possible transitions between topological orders, but it has long been elusive to confirm the mechanism in quantum spin systems. Here, we introduce a concrete spin model that incarnates the mechanism of anyon condensation transition. Our model harbors two topological QSLs in different parameter regions, a non-abelian Kitaev spin liquid (KSL) bilayer state and a resonating valence bond (RVB) state. The bilayer-KSL-to-RVB transition indeed occurs by the mechanism of anyon condensation, which we identify by using parton theories and exact diagonalization studies. Moreover, we observe"anyon confinement"phenomena in our numerical results, akin to the quark confinement in high energy physics. Namely, non-abelian Ising anyons of the bilayer KSL are confined in the transition to the RVB state. Implications and extensions of this study are discussed in various aspects such as (i) anyon-condensed multilayer construction of the Kitaev's sixteenfold way of anyon theories, (ii) additional vison condensation transition from the RVB to a valence bond solid (VBS) in the Kitaev bilayer system, (iii) dynamical anyon condensation in a non-Hermitian Kitaev bilayer, (iv) generalizations of our model to other lattice geometries, and (v) experimental realizations. This work puts together the two fascinating QSLs that are extensively studied in modern condensed matter and quantum physics into a concrete spin model, offering a comprehensive picture that unifies the anyon physics of the Kitaev spin liquids and the resonating valence bonds.

Transitions between QSLs with anyons are less studied, yet outstanding problems in modern condensed mat-ter and quantum physics that promise a profound understanding of quantum entanglement and anyons.A theoretical mechanism for such topological transitions has been formulated using the concept of anyon condensation in their seminal work by Bais and Slingerland [48].In this mechanism, anyon condensation reconstructs manybody quantum entanglement and anyons leading to a new topological phase and new anyons, where the braiding statistics with the condensed anyon plays a key role for the determination of the fate of each anyon, i.e., "confined" vs. "deconfined".The anyon condensation mechanism predicts a variety of transitions among topological phases with anyons, providing global insights on a wide class of topological quantum matter [48][49][50][51][52][53][54][55][56][57][58][59][60][61][62][63][64][65][66].In the field of quantum information, anyon condensation has been an important concept that underlies in designing and manipulating topological quantum codes for faulttolerant quantum computations [67][68][69][70][71]. Unfortunately, it has been elusive to confirm the mechanism in quantum spin systems, i.e., QSL-to-QSL transitions by anyon condensation, because of the scarcity of appropriate microscopic models and the difficulty with defining an order parameter for anyon condensation in terms of local spin operators.
In quantum magnetism, studies of anyon physics have been focused on symmetry breaking transitions from QSLs to long-range orders in frustrated magnets including kagome and triangular lattice antiferromagnets.Z 2 spin liquids such as short-ranged resonating valence bonds and the toric code may be the simplest cases, described by Z 2 gauge theories and four types of anyon excitations: trivial boson with no fractionalization (1), bosonic spinon (e), bosonic vison (m), and fermionic spinon (ϵ = e × m) [5,.Mutual statistics exist between different anyon species.The bosonic spinon carries a spin-1/2 quantum number, thus condensing the bosonic spinon leads to a long-range magnetic order [73,77,87,94].Unlike the spinon, the vison is a spin-0 excitation, preserving time-reversal symmetry.Vison condensation breaks lattice symmetries yielding a crystalline order of spin-singlet dimers, i.e., valence bond solid (VBS) [43,84,90,93,[97][98][99][100].In such transitions to VBS orders (or magnetic orders), the system becomes trivial in the sense that every nontrivial anyon gets either condensed or confined.Essential features of the transitions can be summarized as follows.In this paper, we extend the scope of anyon condensation physics to quantum phase transitions between QSLs.Especially, we construct a spin model that establishes a transition from a QSL (having non-abelian anyons) to another QSL (having only abelian anyons) via anyon condensation.We find such a topological transition in a bilayer system of Kitaev spin liquids.To be specific, our spin model harbors a non-abelian Kitaev spin liquid (KSL) bilayer state and a resonating valence bond (RVB) state in different parameter regions.The RVB state is induced by entangling the KSL bilayer with bonddependent inter-layer interactions similar to the Kitaev interactions (Fig. 1).
The nature of the transition between the KSL bilayer state and the RVB state is identified by computing the condensed anyon and order parameters.In the KSL bilayer state, each layer supports three different anyons: trivial boson (1), non-abelian Ising anyon (σ) and fermion (ψ).There is nontrivial braiding between σ and ψ; a ψparticle sees a σ-particle as a Z 2 flux.In the bilayer-KSLto-RVB transition, the fermion-pair (ψ I ⊠ ψ II ) between layer-I and layer-II is condensed by inter-layer interactions.The anyon physics realized by our spin model is outlined as follows.
  The two points, (i) ψ I ⊠ ψ II is the condensed anyon and (ii) nontrivial braiding exists between ψ and σ in each layer, are the key elements to understand the anyon physics of the bilayer-KSL-to-RVB transition as we shall see in Sec.VI.
Remarkably, our model allows to study the phenomena of "anyon confinement" akin to the quark confinement in high energy physics.Namely, non-abelian Ising anyons (σ I & σ II ) of the KSL bilayer get confined, i.e., they cannot be directly observed in the low energy physics of the RVB state.We confirm the confinement physics by our numerical exact diagonalization calculations (Fig. 2).
Unlike transitions from Z 2 spin liquids to VBS orders, the bilayer-KSL-to-RVB transition does not break any symmetries of the system.The anyon condensation only reconstructs the underlying entanglement structure and supported anyons (i.e., topological order).
Our main results are displayed in Figs. 1, 2, and 9. Figure 1 shows the phase diagram of the model obtained by numerical exact diagonalization.In particular, Fig. 2 highlights the confinement of Ising anyons in the RVB state by demonstrating extremely large energy costs for the excitations.In Fig. 9, we provide a numerical evidence on the existence of the anyon condensation transition.
This work builds a bridge between the two intriguing QSLs extensively studied over the fields of topological quantum matter, quantum magnetism, and quantum information, and offers a comprehensive picture unifying the anyon physics of the Kitaev spin liquids and the resonating valence bonds.

A. Outline of the paper
The paper is organized as follows.In Sec.II, we introduce the bilayer spin model with an emphasis on the conserved quantities.The RVB state and the KSL×KSL state arise in the strong coupling limit and the weak coupling limit of inter-layer interactions, respectively.In Sec.III, an effective quantum dimer model is developed for the RVB state by conducting a sixth order degenerate perturbation theory.We introduce a dimer representation for the original spin model and the hardcore dimer constraint defining the dimer Hilbert space, which plays a crucial role in understanding the anyon condensation and anyon confinement phenomena from our numerical results later.In Sec.IV, we construct a Majorana meanfield theory for the KSL×KSL state and show that timereversal can be spontaneously broken in each layer due to inter-layer interactions.It is shown that both layers have finite energy gaps and nonzero Chern numbers with opposite signs, keeping the whole system achiral.In Sec.V, we investigate the full phase diagram of the model by numerical exact diagonalization methods on a (24+24)-site bilayer cluster.The RVB state and the KSL×KSL state are characterized and distinguished by computing various quantities including entanglement entropy, hardcore dimer constraint, chirality structure factor, and topological degeneracy.Importantly, we find that the numerical results are all consistent with the effective theories in Secs.III and IV.
Section VI is the most important part, where we discuss our main results on anyon condensation and anyon confinement.We start by articulating the non-abelian Ising×Ising topological order of the KSL×KSL state and the abelian Z 2 toric code topological order of the RVB state.We review the anyon condensation mechanism for the transition between the two topological orders.Then, we present numerical evidences for the anyon condensation transition in our model.We clarify the condensed anyon by explicitly calculating associated order parameters.Furthermore, we elucidate the confinement phenomena of Ising anyons by numerically investigating excitation energies of Ising anyons.We will show a simple picture that enables to intuitively understand the anyon confinement from the condensed anyon.
In Sec.VII, we explore generic parameter regimes of the model and show that the emergence of the RVB state or Z 2 spin liquid is determined by the sign structure of the coupling constants of the model.
Section VIII is another important part, where we discuss implications of this work and promising future directions.To name a few, we will discuss (i) the Kitaev's sixteenfold way of anyon theories in the perspective of anyon condensation, (ii) vison condensation transition from the RVB state to a VBS state in our Kitaev bilayer system, (iii) a non-Hermitian Kitaev bilayer and dynamical anyon condensation, (iv) generalizations of our model to other lattice geometries, and (v) experimental realizations.
Appendices A and B provide details of the sixth order degenerate perturbation theory and the quantum dimer model for the RVB state.

II. BILAYER MODEL
We place two copies of the Kitaev honeycomb model [6] on a bilayer geometry of AA stacking as illustrated in Fig. 1(a).In this bilayer setup, our model Hamiltonian consists of three parts: where σ γ and τ γ (γ = x, y, z) are Pauli spins on the upper and lower layers coupled by the bond-dependent Kitaev interactions (K σ , K τ ) and also the inter-layer interactions (G).We label upper and lower layer spins with same site-indices (j, k), and the inter-layer interactions are nothing but the products of adjacent upper-layer and lower-layer Kitaev interactions, The Hamiltonian commutes with the hexagon plaquette operators defined on the upper and lower layers, i.e., [H, Ŵp ] = [H, Ẑp ] = [ Ŵp , Ẑp ′ ] = 0; see Fig. 3(a) for the site convention within a hexagon plaquette p.

III. RVB STATE: STRONG COUPLING LIMIT
To understand the strong coupling limit in an intuitive way, we employ a dimer mapping to a dual kagome lattice.Using the dimer mapping, we derive an effective quantum dimer model and the RVB state.

A. Dimer mapping to a dual kagome lattice
The bilayer model can be regarded as a single layer honeycomb model with four states per site.Each site may have either spin-singlet state |s⟩ or spin-triplet state |t x,y,z ⟩ as shown in the table of Fig. 3.We take a mapping from the honeycomb lattice to a dual kagome lattice.In Fig. 3(a), the dual kagome lattice is constructed by connecting the mid-points of the bonds of the honeycomb lattice.Sites of the honeycomb lattice are now replaced by triangles of the kagome lattice.Interestingly, the kagome lattice has three different bond directions, which are perpendicular to the x, y, z-bond directions of the honeycomb lattice.By using this property, we assign a bond character (x, y, z) to each bond of the kagome lattice [denoted by different colors in Fig. 3 3(b) and 3(c).This kind of dimer mapping was introduced in a recent study of a spin-3/2 transverse field Ising model [47].
By this dimer mapping, spin operators take the following representations.
Other spin operators are obtained by cyclic permutations of (x, y, z) in the above.In this dimer picture, σ γ and τ γ spin operators do monomer pair-creation/annihilation or monomer hopping at the γ-bond of the kagome lattice (γ = x, y, z).

B. Emergence of a quantum dimer model and resonating valence bonds
We develop a degenerate perturbation theory for the strong coupling limit of H.In this case, the unperturbed Hamiltonian is the inter-layer interaction part, Here we employed the composite spin operator, which is Z 2 -valued [(ϕ γ j ) 2 = 1] and satisfies the commutation relation, [ϕ γ j , ϕ λ k ] = 0, and the local constraint, ϕ x j ϕ y j ϕ z j = −1.Interestingly, the four states {|s⟩, |t x ⟩, |t y ⟩, |t z ⟩} are the eigenstates of the composite spin operators (ϕ x , ϕ y , ϕ z ) as summarized in the table of Fig. 3. Hence, the composite spin operators are diagonalized in the dimer representation, e.g., Note that −ϕ z measures the dimer parity across the x, y-bonds, and similarly for −ϕ x,y .The unperturbed Hamiltonian H 0 is readily diagonalized in the dimer representation.We find that the ground state manifold of H 0 is extensively degenerate satisfying the condition, at every bond ⟨jk⟩ γ of the honeycomb lattice.This condition implies the "hardcore dimer" constraint on the kagome lattice, i.e., each site of the kagome lattice is occupied by a single dimer.Specifically, the condition ϕ γ j ϕ γ k = −1 only allows the eight configurations, Note that z measures the dimer parity across the x, ybonds, and similarly for x,y .
The unperturbed Hamiltonian H 0 is readily diagonalized in the dimer representation.We find that the ground state manifold of H 0 is extensively degenerate satisfying the condition, at every bond hjki of the honeycomb lattice.This condition implies the "hardcore dimer" constraint on the kagome lattice, i.e., each site of the kagome lattice is occupied by a single dimer.Specifically, the condition j k = 1 only allows the eight configurations, (10) in every adjacent two triangles.The central site is always occupied by a single dimer because the condition j k = 1 requires the odd dimer parity over the four bonds connected to the central site.Hence, the ground state manifold of H 0 consists of the dimer states that respect the hardcore dimer constraint at every site.The ground state energy is , where N is the number of sites of the honeycomb lattice, and the degeneracy is given by 2 N/2 1 in periodic boundary condition.It must be noted that the ground state manifold of H 0 satisfies the flux condition, lence bo should u guage.A dimer re The full factor, f note tha tor) is i guich, S ble 2 N/2 by the p weight.have an and the state wit later). IV.
The w turbatio tain the , (10) in every adjacent two triangles.The central site is always occupied by a single dimer because the condition ϕ γ j ϕ γ k = −1 requires the odd dimer parity over the four bonds connected to the central site.Hence, the ground state manifold of H 0 consists of the dimer states that respect the hardcore dimer constraint at every site.The ground state energy is , where N is the number of sites of the honeycomb lattice, and the degeneracy is given by 2 N/2−1 in periodic boundary conditions.It must be noted that the ground state manifold of H 0 satisfies the flux condition, The Kitaev terms in H − H 0 generate quantum dimer motions in the ground state manifold of H 0 .By conducting perturbative expansions within the manifold, we obtain the effective Hamiltonian, where Ŵp is the hexagon plaquette operator in Eq. ( 2), and where N is a normalization constant and |Φ⟩ is an arbitrary state in the ground state manifold of H 0 .This state is massively quantum superposed with the uniform zero-flux (W p = Z p = +1).Also, it has the property, ⟨σ γ j σ γ k ⟩ = ⟨τ γ j τ γ k ⟩ = 0, because the Kitaev terms violate the hardcore dimer constraint, Eq. ( 9); see Fig. 13.
The ground state |Ψ⟩ is nothing but a resonating valence bond state on the kagome lattice.To see this, one should understand the effects of Ŵp in the dimer language.Acting on hardcore dimer states, Ŵp gives rise to dimer resonance motions within the 12-site David star: The full list of dimer motions and the associated sign factor, f (D) = ±1, are provided in Table III (also see Appendix B).We note that the above description (apart from the sign factor) is identical to the quantum dimer model by Misguich, Serban, and Pasquier [44].In Eq. ( 13), all possible 2 N/2−1 hardcore dimer configurations are generated by the plaquette operators, and superposed with equal weight.Therefore, in the strong coupling limit of H, we have an effective quantum dimer model H eff (G ≫ K), and the ground state |Ψ⟩ is a resonating valence bond state whose anyon properties are characterized by the Z 2 toric code topological order (discussed later).

IV. KSL×KSL STATE: WEAK COUPLING LIMIT
We construct a degenerate perturbation theory for the weak coupling limit of H. Second order perturbations create effective spin interactions that lead to spontaneous symmetry breaking of time-reversal in the KSL×KSL state.We show this by a Majorana mean-field theory.

A. Effective Hamiltonian
To construct an effective theory for the weak coupling limit, we arrange the Hamiltonian into the form H = H 0 + H 1 (where , and conduct a second order degenerate perturbation theory.The resulting effective Hamiltonian is given by where ∆E(∝ |K σ |) means the energy difference between the ground state and an intermediate excited state of H 0 .Nontrivial effects are generated by the second order term, (16) Among the various combinations of spin operators, we are particularly interested in the combinations that are defined on connected two bonds.For instance, (σ ) on adjacent two bonds, ⟨12⟩ x and ⟨23⟩ y .This eight-spin operator can be simplified to a sixspin operator: By collecting this kind of terms, we arrange the effective Hamiltonian into the following form.

B. Majorana mean-field theory
The effective Hamiltonian is solved by a mean-field theory.To the six-spin operators, we apply the mean-field decoupling, Similar mean-field decoupling schemes may be applied to the other terms that are not explicitly shown in Eq. (18).Yet, those terms merely renormalize the Kitaev term H 0 .Hence, we ignore those terms and just focus on the meanfield Hamiltonian, The three-spin terms, σ α i σ β j σ γ k & τ α i τ β j τ γ k , break timereversal symmetry and create a finite energy gap in the fermion excitations of the Kitaev spin liquid state on each layer [6].
This can be shown by using the Majorana representation for spin-1/2 operator, where the subscript (I,II) means the upper/lower layer.This representation leads to the Majorana mean-field Hamiltonian, where we have chosen the simplest gauge for the uniform zero-flux (u I,ij = ib α I,i b α I,j = +1 and u II,ij = ib α II,i b α II,j = +1 with site i in A sublattice and j in B sublattice), and introduced the mean-field parameters, In this mean-field theory, the two layers are coupled only by the mean-field parameters, M I and M II .Figure 4 shows the result of the self-consistent meanfield calculations.We find that the mean-field parameters have a same magnitude but opposite signs (M I = −M II ).This implies that the KSL state of each layer has a finite energy gap (∆) and a nonzero Chern number (ν) [6]: More importantly, the two layers have opposite signs for the Chern number (e.g., ν I = +1 & ν II = −1).This theory indicates that the KSL×KSL state has the nonabelian Ising × Ising topological order [6,31,55].
On the time-reversal symmetry breaking in the KSL×KSL state, we provide further evidences from our numerical calculations shown in the next section.We will discuss the spontaneous symmetry breaking and the gap-openning problem more carefully by comparing our results with a field theory argument about gapless Dirac fermions in Sec.VIII.

V. EXACT DIAGONALIZATION
Now we investigate the phase diagram of H by exact diagonalization (ED) on the (24+24)-site bilayer cluster in Fig. 1(a).We impose periodic boundary conditions and utilize the flux quantum numbers, W p = ±1 & Z p = ±1, to reduce the size of the Hilbert space in our ED calculations.
Figure 1(c) displays the resulting phase diagram.We find that the ground state appears in the uniform zeroflux sector (W p = Z p = +1).The KSL×KSL and RVB states, which we considered in the weak and strong coupling limits, are stabilized over large regions separated by a single transition at θ c ≃ 0.26π (shown by a small peak in the derivative of the ground state energy, −∂ 2 E gs /∂θ 2 ).We characterize the two QSL states by investigating entanglement entropy, chirality structure factor, hardcore dimer constraint, and topological degeneracy below.

A. Entanglement entropy
Quantum entanglement between the two layers can be measured by the entanglement entropy, where of KSL.By contrast, the RVB state shows strong interlayer entanglement and the dimension of the dimer Hilbert space manifests through the entanglement entropy S layer = log(2 N/2−1 ) in the strong coupling limit [Fig.5(a)].

B. Chirality order
Time-reversal symmetry breaking and the opposite Chern numbers discussed in our mean-field theory for the KSL×KSL state can be checked in our ED via the chirality structure factor, where I, J(= σ, τ ) are layer indices and p, q are plaquette indices.The chiral spin operators, χσ p & χτ p , are defined at each plaquette p as with α ̸ = β ̸ = γ.To be specific, the chiral spin operators can be written as where we followed the site convention shown in Fig. 3(a).The operators χσ p & χτ p are nothing but the (time-reversal odd) gap-opening terms considered by Kitaev [6], which we already have seen in Eq. (20).
Namely, in thermodynamics limit, σ-and τ -layers have a k = 0 order of opposite chiralities, e.g., ⟨ χσ ⟩ > 0 and ⟨ χτ ⟩ < 0. See Figs.This result, together with the entanglement entropy in Fig. 5(a), suggests that in thermodynamic limit the two layers of the KSL×KSL state can be individually described by the effective Hamiltonians [102], where the coupling constants λ σ & λ τ have the opposite signs (λ σ λ τ < 0) as required by the inter-layer antiferrochirality.Notice that the above effective Hamiltonians are exactly same as what we have derived in our Majorana mean-field theory [Eq.(20)].Therefore, timereversal symmetry breaking in the KSL×KSL state is unequivocally shown not only from the mean-field theory but also from the numerical calculations of chirality structure factor.
Our approach of using the chirality-chirality correlation to detect the broken time-reversal in the KSL×KSL state is analogous to the conventional approach of using the spin-spin correlation to distinguish between magnetic orders and spin liquids (see Table I).

C. Hardcore dimer constraint
The hardcore dimer constraint in Eq. ( 9) is another good measure to distinguish between the KSL×KSL and RVB states.Based on our analytical approaches, the two states are expected to show distinct behaviors: ⟨ϕ γ j ϕ γ k ⟩ ≈ ⟨σ γ j σ γ k ⟩⟨τ γ j τ γ k ⟩ in the KSL×KSL state near θ = 0 and ⟨ϕ γ j ϕ γ k ⟩ ≈ −1 in the RVB state around θ = π/2.We confirm the distinct behaviors,

D. Topological degeneracy
Topological degeneracy, i.e., the ground state degeneracy on torus geometry, is a useful probe to detect a topological order in the system.Namely, topological degeneracy tells about the number of anyon types allowed in the system [103,104].We identify "ninefold" degeneracy in the KSL×KSL state and "fourfold" degeneracy in the RVB state (see Fig. 7).
For the analysis of topological degeneracy, we employ four Wilson loop operators, {W σ x , W σ y , W τ x , W τ y }, commuting with the Hamiltonian H.The Hilbert space can then be partitioned into sixteen topological sectors distinguished by the eigenvalues of the Wilson loop operators (W σ,τ x,y = ±1).
The Wilson loop operators are defined by generalizing the hexagon plaquette operators ( Ŵp , Ẑp ) to noncontractible loops of the system.Figure 8(a) shows the cluster used in our exact diagonalization.Periodic boundary conditions are imposed at the cluster boundary (gray hexagon).Along the non-contractible loop in Fig. 8(b), we define the Wilson operator W σ x as follows.
Note that the three operators satisfy the relationship, W σ x W σ y W σ z = +1, in the ground state manifold due to the uniform zero-flux property (W p = Z p = +1).Among the three operators, we use W σ x & W σ y .Similarly, we repeat the same procedure for τ -spins and obtain the following Wilson loop operators.
Again, we have the relationship, W τ x W τ y W τ z = +1, and only use W τ x & W τ y .The Wilson loop operators, {W σ x , W σ y , W τ x , W τ y }, commute with themselves, all the hexagon plaquette operators, and the Hamiltonian H.The sixteen topological sectors of the Wilson loop operators (#1 ∼ 16) are specified in the table of Fig. 7.
In the KSL×KSL state, the ground states appear in the topological sectors having −1 in two of W σ x,y,z and also in two of W τ x,y,z (ninefold degeneracy).In the RVB state, the ground states stay in the topological sectors satisfying the condition, W σ l = W τ l (l = x, y, z) (fourfold degeneracy).

VI. TOPOLOGICAL TRANSITION BY ANYON CONDENSATION AND CONFINEMENT
The topological degeneracies provide useful hints about the underlying topological orders of the two QSLs.The ninefold degeneracy in the KSL×KSL state is consistent with the Ising ×Ising topological order with the nine different anyon-pairs, {1 I , σ I , ψ I } ⊠ {1 II , σ II , ψ II }.Here the subscript means the layer index (I: upper layer, II: lower layer), and each layer has the three anyon sectors: trivial boson (1), Ising anyon (σ), and fermion (ψ).Note that there is nontrivial braiding between σ and ψ; moving a ψ-particle around a σ-particle changes the overall sign of the wave function (|Ψ⟩ → e iπ |Ψ⟩).These anyons satisfy the fusion rules, where the fusion outcome of two σ-particles has two possibilities due to the non-abelian nature of the Ising anyon σ (quantum dimension: √ 2) [6,26,29,31,32].This Ising anyon topological order emerges in each layer of the KSL×KSL state, hence the Ising × Ising topological order for the whole system (bar indicates that the two layers are time-reversal partners) [6,31,32,48,55].
On the other hand, the fourfold degeneracy in the RVB state implies the Z 2 toric code topological order with the four different anyons, {1, e, m, ϵ}; the e-and m-particles are self-bosons with mutual statistics, and the ϵ-particle is a self-fermion composed of the e-and m-particles.All these particles are abelian anyons satisfying the fusion rules, which define the Z 2 toric code topological order [5,6].
The transition between the non-abelian Ising × Ising topological order and the abelian Z 2 topological order can be understood by the mechanism of anyon condensation transition [48,51,52,55].Suppose we condense the fermion-pair ψ I ⊠ ψ II .Then, ψ I ⊠ 1 II and 1 I ⊠ ψ II become indistinguishable and identified as a same type of anyons in the condensed phase: ψ I ⊠ 1 II = 1 I ⊠ ψ II .Moreover, anyons having nontrivial braiding with the fermion-pair are confined, thus only ψ I ⊠ 1 II , 1 I ⊠ ψ II , and σ I ⊠ σ II remain deconfined.The following table highlights the identical anyons (colored) and the confined anyons (crossed) resulting from the anyon condensation.
We note that (44) where we have two trivial bosons (1 I ⊠ 1 II ,ψ I ⊠ ψ II ) and two fermions (ψ I ⊠ 1 II ,1 I ⊠ ψ II ).In order to reproduce the abelian Z 2 topological order in the condensed phase, σ I ⊠ σ II (which has the quantum dimension 2) should be split into two abelian anyons.Therefore, we may assume that Then, Eq. ( 44) becomes identical to the toric code's fusion rules shown in Eq. ( 42) [32,48,55].

A. Anyon condensation
Now we confirm the mechanism of the anyon condensation induced transition in our microscopic model.First, we check if ψ I ⊠ 1 II and 1 I ⊠ ψ II become indeed indistinguishable identical anyons (ψ I ⊠ 1 II = 1 I ⊠ ψ II ).To this end, we consider two different loop operators.The first one is the usual hexagon plaquette operator, which is defined within a single layer [Fig.9(b)].If we apply the Majorana representation, σ γ j σ γ k = −u jk (ic j c k ) where u jk = ib γ j b γ k , the action of W is to create a pair of c-fermions, move one of the fermions along the hexagon, and finally annihilate the fermion-pair [6].Notice that the c-fermions correspond to ψ I fermions.The second one is defined over the two layers: This operator moves a fermion (ψ I ) around an upper "half" hexagon, and moves another fermion (ψ II ) around a lower "half" hexagon [Fig.9(c)].This type of loop operator has been considered in a recent study on anyon condensation in kagome quantum spin liquids [61].Figure 9(a) shows the expectation value ⟨L⟩.We find that ⟨L⟩ is small in the KSL×KSL phase (close to zero near θ = 0), meaning that ψ I and ψ II cannot move across the two layers because they are distinct anyons supported on different layers.By contrast, substantially large values of ⟨L⟩ are observed in the RVB phase (reaching one at θ = π/2).This implies that ψ I and ψ II become identical particles, thus the fermions (ψ I = ψ II ) can complete the hexagon loop motion.It is only possible when we have the condensation of ψ I ⊠ ψ II .Therefore, ⟨L⟩ plays a role of an order parameter for the anyon condensation.

B. Anyon confinement
Next, we investigate the confinement of the Ising anyons, To allow Ising anyon excitations, we consider flux sectors with W = −1 or Z = −1.Under broken time-reversal symmetry, flux excitations with W = −1 or Z = −1 trap localized Majorana zero modes, behaving as the nonabelian Ising anyons [6,31,36].By comparing energy costs of different flux patterns, we can identify in which flux patterns the confinement occurs.
Figure 2 shows the lowest excitation energy profiles in various flux patterns.Depending on the phase, flux excitations have different behaviors in the energy cost.In the KSL×KSL phase, the excitation energy with respect to the (zero-flux) ground state is roughly proportional to the number of excited fluxes.In sharp contrast, in the RVB phase, such simple counting does not work.To be specific, when fluxes are excited only on a single layer, the excitation energy dramatically increases across the transition point θ c ≃ 0.26π [see the flux pattern #1 in Figs.2(a The large energy costs due to unpaired fluxes, shown in the flux patterns #1,3,4, can be explained by the hardcore dimer constraint of the RVB phase (⟨ϕ γ j ϕ γ k ⟩ ≈ −1).Any unpaired flux (W p Z p = −1) necessarily violates the hardcore dimer constraint leading to an energy increase in the order of G.In contrast, paired fluxes between the two layers at same locations (W p = Z p = −1) can be consistent with the hardcore dimer constraint.The lowest excitations of such paired fluxes are captured by the effective quantum dimer model which we derived in large-G limit [Eq.(12)].Notice that the excitation energy scale is λ (∝ K 6 σ /G 5 = cos 6 θ/ sin 5 θ), which explains the vanishingly decreasing excitation energies of the flux patterns #2,5-8 [Figs.2(e) and 2(f)].Remarkably, unpaired fluxes occur at an extremely high energy scale (G) compare to paired fluxes described by the quantum dimer model (λ).In the low energy physics of the RVB phase, unpaired fluxes never appear, essentially confined due to the extremely large energy cost.
The paired fluxes in Figs.The fact that only paired fluxes appear as deconfined anyons in the RVB phase can be understood in a more intuitive way.In order for ψ I & ψ II to be identical particles, ψ I & ψ II have to see a same flux pattern on the upper and lower layers.Otherwise, ψ I & ψ II cannot be identical.This simple fact results in the confinement of unpaired fluxes (σ I & σ II ), only allowing paired fluxes (σ I ⊠ σ II ) as deconfined anyons.

C. Topological degeneracy revisited
The fourfold topological degeneracy of the RVB state is also easily understood by the same picture.Fermions, ψ I & ψ II , must see a same pattern of Wilson loop fluxes on the upper and lower layers.In other words, the topological sectors allowed by the ψ I ⊠ ψ II -condensation must satisfy the condition, W σ l = W τ l (l = x, y, z), which leads to exactly the four sectors (#1, 6, 11, 16) in Fig. 7.

VII. GENERAL CASES
We investigate other parameter regions of the model beyond Eq.( 3).We still constrain ourselves into the cases of |K σ | = |K τ | since more generic cases can be easily understood from our discussion here.We find that the signs of the coupling constants {G, K σ , K τ } determine the emergence of a RVB-type Z 2 spin liquid.Table II lists six possible cases for the sign structure.Three cases (#1,5,6) stabilize a Z 2 spin liquid (Z 2 SL) whereas the other three cases do not.The case #1 corresponds to the parameter choice in Eq. ( 3).
Here we discuss the cases #5 and #6 (G < 0, K σ = K τ ). Figure 10 shows ED results for the two cases with the parametrization, The results are almost same as in the case #1 in terms of the ground state energy E gs , the topological degeneracy, the entanglement entropy S layer , and the four-spin correlator ⟨ϕ γ j ϕ γ k ⟩.Yet, there is a difference between #5,6 and #1 in the sign of ⟨ϕ γ j ϕ γ k ⟩, which is positive in the former cases but negative in the latter case.In the strong coupling limit (|G| ≫ |K σ,τ |), we see distinct behaviors in the associated Z 2 spin liquids: The property ϕ γ j ϕ γ k = −1 of #1 allows us to construct the effective quantum dimer model on the dual kagome lattice as we have seen in Sec.III.However, the property ϕ γ j ϕ γ k = +1 of #5,6 leads to states violating the hardcore dimer constraint, so the dimer representation is no longer useful for the two cases.
Instead, a Z 2 gauge theory can be more conveniently constructed on the honeycomb lattice for the strong coupling limits of the cases #5,6.Utilizing the property ϕ γ j ϕ γ k = +1 at each bond ⟨jk⟩ γ , we may define a Z 2 link variable, Such X-variables are subject to the Gauss law constraint, which is simply the product of X-variables on the three links sharing site j.The Gauss law constraint arises as a combined effect of (i) the property ϕ γ j ϕ γ k = +1 at each bond ⟨jk⟩ γ and (ii) the constraint ϕ x j ϕ y j ϕ z j = −1 at each site j.In other words, the property ϕ γ j ϕ γ k = +1 defines the Hilbert space of the Z 2 link variables {X jk } subject to the Gauss law constraint.
We repeat a sixth order degenerate perturbation theory for the cases #5,6 and obtain an effective Hamiltonian that is exactly in the same form of Eq. ( 12): The plaquette operator Ŵp flips the signs of X-variables at plaquette p as denoted by the Pauli Z-operators in the above.Here we see a Z 2 gauge theory on the honeycomb lattice.
In the other cases (#2,3,4), the effects of the K σ,τ terms completely cancel each other, failing to create any Z 2 spin liquid state in the strong coupling limit.For instance, in the cases #2,3 (G > 0 and K σ = K τ ), we have the relationship, where |Φ⟩ is an arbitrary ground state (ϕ γ j ϕ γ k = −1) of the interlayer interactions.Then, one can easily see that Due to this cancellation, the Kitaev interactions do not create any effect on the dimer Hilbert space and any Z 2 spin liquid.

VIII. DISCUSSION AND OUTLOOK
The anyon condensation transition between the nonabelian KSL × KSL state and the RVB state was identified in our model via the loop operator L. Interestingly, there is an intimate relationship between the anyon condensation and the hardcore dimer constraint [Eq.( 9)].
To see this, we consider the constraint in a slightly different fashion: σ γ j σ γ k = −τ γ j τ γ k .This tells us that moving the fermions, ψ I & ψ II , has a same effect, meaning that the RVB state does not distinguish between ψ I and ψ II ; they are essentially same particles.To be more precise, if we apply the hardcore dimer constraint to Eq. ( 47), we obtain the relationship L = W = Z.Therefore, the hardcore dimer constraint itself means the ψ I ⊠ψ II anyon condensation.
We emphasize that the inter-layer interactions (σ γ j σ γ k τ γ j τ γ k = ϕ γ j ϕ γ k ) are crucial for realizing the two topological spin liquids in the bilayer system and also for understanding the anyon condensation transition.In the weak coupling regime, a chirality order is developed in each layer due to the inter-layer interactions, yielding the non-abelian KSL×KSL state.In the strong coupling regime, the inter-layer interactions construct the Hilbert space of hardcore dimers for the RVB state.The hardcore dimer constraint established by the inter-layer interactions plays a role of an order parameter for the anyon condensation transition between the two topological spin liquids.Furthermore, it was due to the inter-layer interactions conserving the flux quantum number that we could easily trace the Ising anyons and their confinement across the transition (Fig. 2).
This work demonstrates an anyon-condensed multilayer construction for the Kitaev's sixteenfold way of anyon theories.In his original work, Kitaev provided a simple idea to construct sixteen different topological orders [6].Namely, in the system of Majorana fermions coupled with Z 2 fluxes, the spectral Chern number (ν) of the Majorana fermions determines the anyon properties of the Z 2 fluxes.Changing the Chern number results in sixteen distinct types of topological orders determined by ν mod 16.Such sixteenfold way constructions have been discussed in several solvable models recently [105][106][107][108].Here we discuss a different approach using anyon condensation in mutilayer systems.First, we note that the non-abelian KSL with ν = 1 realizes the Ising anyon topological order.Then, the bilayer system of KSLs with the opposite Chern numbers ν I = +1 & ν II = −1 (having zero for the net Chern number: ν = ν I + ν II = 0) may construct the Z 2 toric code topological order.Our bilayer model demonstrates this via the RVB state that emerges from the KSL×KSL state by anyon condensation.Generalizing this idea, all the sixteen topological orders could be constructed by stacking non-abelian KSLs with ν = ±1 and inducing anyon condensation over the multilayers, i.e., anyon-condensed multilayer constructions of the Kitaev's sixteenfold way (Fig. 11).
In our study, the KSL×KSL state was shown to open a finite energy gap by developing the chirality order that breaks time-reversal symmetry.From renormalization group analyses for interacting Dirac fermions, any shortranged four-fermion interaction is, however, perturbatively irrelevant in the systems with time-reversal and lattice symmetries.Hence, the gaplessness of the Dirac fermions is maintained until some finite strengths of fourfermion interactions [109].The same physics applies to the Kitaev spin liquid where the low energy physics is described by the Majorana version of the field theory.Such gaplessness is expected around the weak coupling limit of our system since the inter-layer interactions are essentially four-fermion interactions in the Majorana representation [σ γ j σ γ k τ γ j τ γ k = u I,jk u II,jk (ic I,j c I,k )(ic II,j c II,k ), where I & II are layer indices].Nonetheless, our meanfield theory and finite-size ED calculations suggest that the bilayer system can spontaneously break time-reversal symmetry and open a finite energy gap before the transition to the RVB state.Combining our results with the field theory argument, we anticipate that in the thermodynamic limit the bilayer system has a gapless region from the θ = 0 point, then opens up a gap at some finite value θ gap , and finally undergoes the anyon condensation transition to the RVB state at θ c (0 < θ gap < θ c ).
Low energy field theory for the anyon condensation transition (at θ c ) is an interesting problem lying beyond the Landau-Ginzburg-Wilson paradigm.To our best knowledge, a general framework for such transitions is currently unknown.Nonetheless, the 3D Ising universality has been proposed for the ψ I ⊠ψ II condensate induced transition between the Ising × Ising and Z 2 topological orders by Burnell, Simon, and Slingerland [51,52].Our system is expected to be in the same 3D Ising universality class.
We compare our work with the recent spin-3/2 transverse field Ising model by Verresen and Vishwanath [47].1)] extended by the interactions in Eq. ( 55).Here we assume Kσ = 1 & Kτ = −1.In large-G limit, it is expected that the RVB-to-VBS transition occurs around D ∼ λ(∝ K 6 σ /G 5 ) and the resulting VBS state is the 12-site pinwheel VBS state that is well known from kagome antiferromagnetism [9,12].
Both cases lead to emergent quantum dimer models with similar structures.The hardcore dimer constraint for the dimer Hilbert space is energetically implemented, and the dimer resonance is induced by anyon fluctuations-by the Kitaev interactions in our case, and by the transverse field in their case.Despite the similarities, the two works have different interests.Ref. [47] focuses on the possible quantum liquids induced by different types of anyon fluctuations.In our study, we are mainly interested in the anyon condensation transition and its identification.
Another anyon condensation transition can be generated in our bilayer system if we add the interactions, where next-nearest-neighbor composite spins are coupled in a bond-dependent way.If the composite spins are connected by α-and β-bonds, then only the γ-components are coupled (α ̸ = γ ̸ = β).In the dimer Hilbert space (appearing in large-G limit), such interactions become the dimer interactions, and the values of the energy coefficient E D are listed in Table III.The emergence of the dimer interactions implies that there should be a transition from the RVB state to a valence bond solid (VBS) state, i.e., a crystalline order of dimers.The RVB-to-VBS transition is driven by the condensation of m-particles (or visons) of the RVB state [89,90].By the vison condensation transition, the topological spin liquid (RVB state) becomes trivial loosing all nontrivial anyons.Namely, e-and ϵ-particles are confined due to their nontrival braiding with m-particles.
An anticipated phase diagram is schematically drawn in Fig. 12.
Our model can be generalized to other tri-coordinated lattices beyond the honeycomb lattice as in the original Kitaev model [110].
In particular, the star lattice (left figure) would be an interesting case where the pure Kitaev model itself allows a chiral spin liquid state as shown by Yao and Kivelson [111].Such spontaneous time-reversal breaking may have some novel effects on the RVB side and the supported anyon types.Generalization of our work to other lattice geometries would be an interesting direction for future studies.
The bilayer model has an interesting connection to a non-Hermitian Kitaev system.Motivated by inevitable environment effects such as decoherence and dissipation on quantum many-body states prepared in quantum devices, there has been growing interest in open quantum systems recently [112][113][114][115][116][117].Such environment effects can be investigated by using the so called Choi-Jamio lkowski isomorphism, which transforms the system's density matrix (ρ) into a doubled state vector (|ρ⟩⟩) that is governed by a non-Hermitian Schrödinger equation: The problem of a Kitaev spin liquid coupled to an environment becomes identical to a bilayer Kitaev system featured with non-Hermitian inter-layer interactions, i.e., a non-Hermitian analog of our bilayer model [117].More interestingly, anyon condensation phenomena dynamically occur in the non-Hermitian Kitaev system by the non-Hermitian inter-layer interactions [117].Quantum spin liquids coupled to environments (open quantum spin liquids) are an interesting setup which allows novel anyon physics via non-unitary dynamics.
Lastly, we comment on experimental realizations of our system.There have been several concrete proposals upon realizing the Kitaev's honeycomb model and toric code model by using ultracold atoms or Rydberg atoms together with Floquet engineering [118][119][120][121][122].Moreover, non-abelian topological orders and anyons have been successfully simulated in superconducting-qubit and trapped-ion quantum processors recently [19][20][21].We expect that such quantum simulators with high controllability, especially the reconfigurable atom arrays [23]  of our bilayer system in future.

FIG. 1 .
FIG. 1. Bilayer spin model.(a) Honeycomb lattice bilayer with the AA stacking.Pauli spins, σ and τ , reside on the upper and lower layers coupled by the intra-layer Kitaev interactions (Kσ and Kτ ) and the four-spin inter-layer interactions (G).The x, y, z-bonds are denoted by different colors of red, green, and blue.The figure depicts the (24 + 24)-site bilayer cluster used in the numerical exact diagonalization.(b) Illustration of the four-spin inter-layer interactions.(c) Phase diagram of the model.The KSL×KSL and RVB states are connected by a continuous transition at θc ≃ 0.26π; indicated by a small peak in the derivative of the ground state energy, −∂ 2 Egs/∂θ 2 .

FIG. 3 .
FIG. 3. Dimer mapping to the dual kagome lattice.(a) Mapping from the honeycomb lattice to the kagome lattice.Different colors denote the three bond characters, x(red), y(green), and z(blue).The numbers (1 ∼ 6) indicate the site convention within a hexagon plaquette p. (b),(c) Illustrations of several dimer states.The four local states {|s⟩, |tx⟩, |ty⟩, |tz⟩} and the dimer mapping are defined in the bottom table.

FIG. 4 .
FIG. 4. Result of the Majorana mean-field theory.The meanfield parameters (MI,II) are obtained by solving Eq. (23) selfconsistently.The Kitaev term is fixed by |Kσ| = |Kτ | = 1.The result does not depend on the relative sign of the coupling constants.
Flat structure factors shown in the RVB state indicate short-ranged correlations but essentially no long-range order in the chirality.See Figs.6(c), 6(d), 6(g), and 6(h).

FIG. 9 .
FIG. 9. Condensation of the fermion-pair.(a) The expectation values of the loop operators, ⟨W ⟩ and ⟨L⟩.In the RVB state, the substantially large values of ⟨L⟩ indicate the condensation of the fermion-pair, ψI ⊠ ψII.(b),(c) Visualizations of the W and L operators.

ρ
mn |m⟩⟨n| ⇒ |ρ⟩⟩ = m,n ρ mn |m⟩ ⊗ |n⟩.(57) (a)].In this mapping, the four local states {|s⟩, |t x ⟩, |t y ⟩, |t z ⟩} on the honeycomb lattice are represented by four dimer states {| ⟩, | ⟩, | ⟩, | ⟩} on the dual kagome lattice.The state |s⟩ is the no-dimer state at the corresponding local triangle of the kagome lattice while the state |t γ ⟩ is the dimer state occupying the γ-bond of the local triangle; see examples shown in Figs.

TABLE II .
Possible sign structures of the coupling constants {G, Kσ, Kτ } and the emergence of a RVB-type Z2 spin liquid state.
, could realize and further extend the rich anyon physics

TABLE III .
(56)r motion graphs and dimer interaction energies.Top: each graph depicts a dimer state (red) and the closed path (light blue) of the dimer motion by Ŵp.Middle: the sign factor f (D) in Eq. (14).Bottom: the energy coefficient ED of Vp [Eq.(56)].Other cases related by symmetry are dropped for simplicity.