Magnetic fragmentation and fractionalized Goldstone modes in a bilayer quantum spin liquid

We study the phase diagram of a bilayer quantum spin liquid model with Kitaev-type interactions on a square lattice. We show that the low energy limit is described by a $\pi$-flux Hubbard model with an enhanced SO(4) symmetry. The antiferromagnetic Mott transition of the Hubbard model signals a magnetic fragmentation transition for the spin and orbital degrees of freedom of the bilayer. The fragmented"N\'eel order"features a non-local string order parameter for an in-plane N\'eel component, in addition to an anisotropic local order parameter. The associated quantum order is characterized by an emergent $\mathbb{Z}_{2} \times \mathbb{Z}_{2}$ gauge field when the N\'eel vector is along the $\hat{z}$ direction, and a $\mathbb{Z}_2$ gauge field otherwise. We underpin these results with a perturbative calculation, which is consistent with the field theory analysis. We conclude with a discussion on the low energy collective excitations of these phases and show that the Goldstone boson of the $\mathbb{Z}_{2} \times \mathbb{Z}_{2}$ phase is fractionalized and non-local.

We study the phase diagram of a bilayer quantum spin liquid model with Kitaev-type interactions on a square lattice. We show that the low energy limit is described by a π-flux Hubbard model with an enhanced SO(4) symmetry. The antiferromagnetic Mott transition of the Hubbard model signals a magnetic fragmentation transition for the spin and orbital degrees of freedom of the bilayer. The fragmented "Néel order" features a non-local string order parameter for an in-plane Néel component, in addition to an anisotropic local order parameter. The associated quantum order is characterized by an emergent Z2 × Z2 gauge field when the Néel vector is along theẑ direction, and a Z2 gauge field otherwise. We underpin these results with a perturbative calculation, which is consistent with the field theory analysis. We conclude with a discussion on the low energy collective excitations of these phases and show that the Goldstone boson of the Z2 ×Z2 phase is fractionalized and non-local.
Quantum spin liquids (QSLs) are frustrated magnets that do not exhibit long range magnetic order down to zero temperature [1][2][3][4]. Quantum fluctuations in these systems give rise to exotic phenomena such as fractionalization and long-range entanglement, which now become the defining properties for QSLs [5][6][7]. The Kitaev model on the honeycomb lattice [8] is one of the few examples of an exactly solvable model with a QSL ground state (GS). In recent years, remarkable progress in identifying candidate materials with strong Kitaevtype interactions has been achieved, in such instances as the iridates [9,10] and α-RuCl 3 [11]. Kitaev interactions may also be strong in other van der Waals (vdW) materials [12]. Bilayers and moiré superlattices of vdW materials are new tunable quantum platforms for realizing a multitude of novel phases, with a variety of basic building blocks including graphene [13], semiconductors [14] and superconductors [15].
Motivated by these developments, we study the phase diagram of a bilayer QSL model with Kitaev-type interactions on a square lattice (see Fig 1(a)). First introduced in Ref. 16, the exact ground state of the monolayer model is an algebraic QSL featuring two flavors of Majorana fermions that are delocalized on the π-flux square lattice, and gapped π-flux (vison) excitations [17]. In the bilayer model Eq. (1), we add an Ising-type interlayer spin interaction, which commutes with the intra-layer flux operators and hence allows for controlled calculations. Our main results are summarized as follows: (i) Below the vison gap, we map the low-energy subspace of the bilayer model to a π-flux Hubbard model at half-filling, with an emergent SO(4) symmetry. Monte Carlo studies of this model show an antiferromagnetic (AFM) Mott transition at critical U c ∼ 6t [18]. (ii) The in-plane components of the AFM order parameter, i.e. n x,y of the Néel vector n, correspond to non-local order parameters whereas the out-of-plane component (n z ) is a local order parameter in terms of the spin and orbital degrees of freedom (DOF) of the bilayer system. (iii) The system features a Z 2 × Z 2 gauge field when the Néel vector points alongẑ, and a Z 2 gauge field otherwise. (iv) To complement the results of the Hubbard model, we perturbatively derive an effective Hamiltonian in the limit of large interlayer interactions. We confirm the magnetic fragmentation and topological degeneracy directly in terms of the original DOF, which are consistent with the Majorana fermion representation of the spin model. (v) We show that the Goldstone modes of the fragmented AFM order is fractionalized in the Z 2 × Z 2 phase, in comparison to the normal Goldstone modes in the Z 2 phase.
Microscopic model. One of the key conditions for the exact solution of the Kitaev model is the anticommutation relations of the Pauli matrices, {σ i , σ j } = 2δ ij . Since there are only three Pauli matrices, this method can only be applied to lattices with coordination number z = 3 such as honeycomb, hyperhoneycomb and hyperoctagon lattices. However, it is possible to extend Kitaev's method to Γ matrices that obey the Clifford algebra {Γ i , Γ j } = 2δ ij [19,20]. For instance, for a fourdimensional representation of the Clifford algebra, there are five Γ α operators along with ten Γ αβ = i 2 [Γ α , Γ β ] and an identity matrix, which span the local Hilbert space. Therefore, Kitaev's construction can be extended to lattices with coordination number up to z = 5 [19,20]. We adapt this representation and consider the intra-layer Hamiltonian [16], H K = − ij γ ,ν K ν (Γ γ νi Γ γ νj +Γ γ5 νi Γ γ5 νj ), where ν = 1, 2 is the layer index and γ is the type of the bond as depicted in Fig. 1. We also introduce an inter-layer Ising interaction H J = J i Γ 5 1i Γ 5 2i . The full Hamiltonian can be expressed in terms of spins (σ) and orbital (τ ) Pauli matrices using the relation Γ α = −σ y ⊗ τ α (α = x, y, z), Γ 4 = σ x ⊗ I 2 and Here τ γ = τ x , τ y , τ z , I for γ = 1, 2, 3, 4 respectively, corresponding to the 4 bonds incident on a vertex of the square lattice as shown in Fig. 1(a) and the sum is over all the γ bonds. Note that the γ = 4 (yellow) bond, which we refer as the 'identity' bond henceforth, has trivial orbital dependence. We consider K 1 = K 2 = K, unless specified otherwise. We identify two inequivalent intralayer flux plaquette operators W νp = σ z νk σ z νn τ x νi τ y νj τ x νk τ y νn and W νp = σ z νk σ z νn τ x νn τ y νk τ x νl τ y νm , each with ± 1 eigenvalues. Both types of plaquette operators commute with the Hamiltonian and the Hilbert space is divided into sectors of conserved fluxes. Note that the Ising form of the inter-layer exchange is crucial to preserve [W νp/p , H] = 0 [17]. The intra-layer Hamiltonian can be solved by using a Majorana fermion representation of the Γ matrices [16], where u γ ν,ij = ib γ νi b γ νj (see supplemental material (SM) for details). This representation is redundant and the physical states in each layer must be restricted to the eigenstates of D νj = ib 1 νj b 2 νj b 3 νj b 4 νj c x νj c y νj , with eigenvalues 1. As in the Kitaev model, these constraints are imposed by the projection operator P ν = i (1 + D νi )/2. The intra-layer bond operators, u γ ν,ij commute with H K and therefore are conserved with eigenvalues ±1. A Z 2 gauge transformation at site i for layer ν involves flipping the signs of the Majorana fermions and bond operators, c α νi → −c α νi ; u γ ν, ij → −u γ ν, ij . We combine the Majorana fermions on the two layers to form complex fermions, ]. According to Lieb's theorem [21], the GS manifold of H K lies in the π-flux sector and consequently the eigenvalue of W νp/p = p/p u γ νij is −1 in any GS configuration, for all square plaquettes. The spectrum is given by E K = ±4K cos 2 k x + sin 2 k y which includes two inequivalent Dirac points at (± π 2 , 0). Next, we represent the inter-layer interaction in terms of the Majorana fermions: . H J commutes with the intra-layer flux operators W p/p . However, the quartic form of the inter-layer exchange precludes the exact solvability of H, which can be expressed as where n νi = f † νi f νi . Enhanced emergent symmetry. The Hamiltonian in eq. 2 has a global U (1) symmetry in each layer (ν = 1, 2), e −iθ i σ z νi He iθ i σ z νi = H, a Z 2 layer exchange symmetry X , a particle-hole symmetry C, and a time reversal symmetry T . This results in a full symmetry group (1), as detailed in SM. After the Majoranization, the two U (1) rotations manifest themselves as: where κ = −1, 1 for ν = 1, 2 respectively. This can be viewed as "charge": U c = e iθ νi f † νi fνi and "pseudospin": U s = e iθ νi κf † νi fνi rotations where κ = +1(−1) for ν = 1(2). The particle-hole symmetry C and U (1) charge rotations form the O(2) c subgroup, while the layer exchange X and U (1) pseudo-spin rotations form the O(2) s subgroup. Next, we fix the gauge by choosing u γ ν,ij = u ij , for both the layers and pick the π-flux configuration as discussed above. The resulting lowenergy Hamiltonian (below the flux/vison gap) is a πflux Hubbard model at half-filling with a hopping amplitude t = 2K and interaction strength U = 4J. It is  (1) is spontaneously broken down to a different subgroup H ≤ G, with an order parameter manifold M = G/H. The gauge group for the associated topological order in each phase is also listed.
well established [22] that the Hubbard model on a bipartite lattice possesses an enhanced G = SO(4) The equivalence established above shows that our model also exhibits an enhanced SO(4) symmetry at the low energy sector. In fact, this emergent SO(4) symmetry exists in any subspace with a fixed flux configuration. Emergent symmetries can play a key role to describe the low energy physics of strongly correlated systems including cuprates [23] and iron pnictides [24].
Quantum Monte Carlo studies have shown that the repulsive (J > 0) π-flux Hubbard model displays a phase transition from Dirac semimetal to an AFM Mott insulator at J/K ∼ 3 (U/t ∼ 6) [25] (see Fig. 1(c)). Moreover, due to J → −J mapping in the Hubbard model, the phase diagram is symmetric for ferromagnetic (FM) and AFM inter-layer exchange, for which the Néel order maps to superconducting and charge density wave orders.
The Néel vector of the AFM order where N is the number of sites and r ix(y) is the x(y) coordinate of site i, can point along any direction on the Bloch sphere. Goldstone modes always arise since the emergent symmetry G = SO(4) × Z T 2 is spontaneously broken down to H = SU (2) c × U (1) s ZT 2 in the Néel order (see SM for details). However, as we show below, different orientations of the Néel vector correspond to distinct ground states with different symmetry and topological properties, as summarized in Table I. (i) The Néel vector points along the z-direction, n ẑ, In terms of Majorana fermions, eq. 5 takes the form , hence corresponding to a physical operator (−1) rix+riy (σ z 1i − σ z 2i ). In other words, n z is a local order parameter of a Landau-type long range order.
The Z 2 gauge fields for the two layers, u 1,ij and u 2,ij , are decoupled, leading to a Z 2 × Z 2 topolog-ical order described by 4-component Abelian Chern-Simons theory [26] characterized by matrix K = 0 2 2 0 ⊕ 0 2 2 0 . However, the Goldstone mode of the Néel order is not a gauge-invariant quantity, but instead an anyon obeying mutual semion statistics with the vison in each layer. More precisely, the above Goldstone mode carries the gauge charge for the Z 2 gauge field from each layer. Incorporating the gapless anyon b in (6) into the low energy description, the effective field theory for this algebraic spin liquid reads where A α=c,s µ label the charge and pseudo-spin external gauge fields, and are the charge and pseudo-spin vectors [26] for the Chern-Simons theory.
(ii) The Néel vector lies in-plane, e.g. n ⊥ẑ with Unlike n z , the in-plane components, n x and n y are not gauge invariant as the local gauge transformation maps n x(y) → −n x(y) . However, a non-local gauge invariant correlator can be defined [27].
where the gauge string for fermions, B(r, r ) = ij∈(r,r ) u 1ij u 2ij , connects operators at the end sites (r, r ). The value of C x(y) (r, r ) is the same in all gauge choices. Therefore, the ground state, symmetrized over all gauge configurations through the projection procedure, also has the same value of C x(y) (r, r ), signifying a string order parameter. Physically, the long-range string order corresponds to the condensation of anyon b in the field theory (7), hence breaking the gauge group down to Z 2 via the Higgs mechanism [28].
An alternative way to understand the gauge structure is to notice the following local order parameter for the in-plane Néel order for a pair of nearest neighbor sites i, j . Due to the mutual braiding phase of e iπ between a vison and a fermion in each layer, a vison from layer 1 (or 2) is nothing but a vortex for the above local order parameter, since S + i,j acquires a e ±2πi phase as it travels around a vison from layer 1 (2). The logarithmic confinement of vortices in the in-plane Neel phase suggest that the vison from layer 1 (or 2) is confined, therefore reducing the Z 2 × Z 2 gauge group down to Z 2 . A similar conclusion can be drawn if the Néel vector has both in-plane andẑ components.
In general, the ground state can have both non-zero out-of-plane (n z ) and in-plane (n x , n y ) components. The term, 'magnetic fragmentation' is coined for phases that display a coexistence of a local Landau-type order parameter and a non-local topological order [29][30][31][32][33][34]. Magnetic fragmentation is theoretically predicted [29] and experimentally observed [30] in spin ice materials such as Nd 2 Zr 2 O 7 where a local AFM order coexists with a spin liquid with FM correlations. The conclusions we draw from the Hubbard model rely on mean-field order parameters. Next, we underpin these results by a perturbative analysis.
Perturbative analysis in the limit of large interlayer exchange. We corroborate the results of the Hubbard model (Eq. 2) by considering the bilayer in the large-J limit, without reference to the Majorana representation (Eq. 1), on a torus. We introduce effective pseudo-spin and orbital DOF appropriate to this limit. We next derive effective models on the large-J GS manifold, to fourth order in the intra-layer coupling K. By analogy to the Hubbard model, we distinguish between cases with a) Z 2 and b) Z 2 × Z 2 topological order. For a), we show that the GS manifold is a state of uniform π flux, with a finite C x(y) correlator (Eq. 10). We also demonstrate that the GS manifold has four-fold topological degeneracy and that the visons are confined. For b), we also obtain a GS with uniform π flux, which has sixteen-fold topological degeneracy and deconfined vison excitations. These results naturally lead to the conclusion that the two phases are separated by a topological phase transition.
Effective degrees of freedom. We first introduce the effective pseudo-spin and orbital DOF. For K = 0 and finite FM inter-layer interactions (J < 0), the spins on overlapping sites form GS doublets |↑ 1 ↑ 2 , |↓ 1 ↓ 2 . These can be represented by a bilayer pseudo-spin obeying an SU(2) algebra. In addition, the orbital DOF for each pair of overlapping sites form a four dimensional Hilbert space, corresponding to one singlet and three triplet configurations. To represent these states, we introduce the inter-layer orbital operators q γ i = τ γ 1i τ γ 2i , γ = x, y, z, which mutually commute as [q α i , q β i ] = 0. The four orbital states can be labeled by the three eigenvalues q γ i = ±1, constrained to obey γ q γ i = −1. The Hilbert space thus includes all states of the form with the implicit local constraint. Pairs of nearestneighbor q γ i/j define the bond variables which take on values of ±1 for γ ∈ {x, y, z}, while they are trivially equal to 1 for additional identity bonds, labeled by ρ γ=I ≡ 1. To any {q γ i } configuration, we can associate a unique {ρ γ ij } bond configuration, while the converse is not true. Orbital states like |φ 0 = |∀ q γ i = −1 , which have uniform ρ γ ij = 1, play an important role in all subsequent discussions.
Defects in |φ 0 take the form of strings of negative bonds as shown in Fig. 2 (b). Defects in both pseudo-spin and ρ γ ij bonds are introduced by operating with the flux operators W νp/p . Each of these flips the pseudo-spin components along x/y as η x/y n,k → −η x/y n,k in the corresponding unit cell (2 (a)). Each also changes the signs of all six ρ γ bonds connected with sites n, k. Note that any string defect cannot be eliminated by application of W νp/p operators.
Effective Hamiltonian. The effective Hamiltonian H η−ρ , projected onto the K = 0 GS manifold, reads where is obtained at second order in K (g 2 = K 2 /4J). Note the distinction between NN, in-plane pseudo-spins connected via variable and trivial identity ρ γ ij bonds, respectively. We introduce small perturbations h x/z > 0 to explicitly break the continuous symmetry, enforcing the staggered pseudo-spin configurations along x/z, respectively. As it turns out, these respectively correspond to Z 2 and Z 2 ×Z 2 topological order. The fourth-order contribution is where g 4 = K 4 /J 3 . η p and η p include linear combinations of products of η and q γ operators around p/p plaquettes. H g2 commutes with all flux operators for h x → 0, while this always holds for H g4 . Technical details and derivations related to the effective Hamiltonian and the following sections are relegated to SM.
Z 2 topological order. We consider H η−ρ with h x > 0 and h z = 0, on a torus. Exact diagonalization calculations indicate that the GS manifold of H g2 includes |η xx , φ 0 , with η xx denoting a finite staggered pseudo-spin along x (see SM). Importantly, any configuration with string defects in the ρ γ bonds are gapped, with an energy cost which scales as the string length. We next consider the evolution of the GS manifold at H g4 level for h x g 4 g 2 . In SM, we show that exact diagonalization calculations indicate H g4 projects |η xx , φ 0 onto a state of uniform π flux per plaquette: Note thatη xx is a state of staggered pseudo-spins including corrections at both H g2 and H g4 levels. This is a state of definite π flux since W νp (1−W νp ) = −(1−W νp ). Moreover, all W νp/p commute with the operator Consequently, the latter has a finite expectation value in |Ψ GS for any pair of i, j, reflecting a locking of pseudospin and ρ γ ij bond configuration. Moreover, D x ij is equivalent to a gauge-invariant correlator of the Hubbard model corresponding to an in-plane Néel vector (Eq. 10).
Note that any state with an open string defect, obtained by first including strings in φ 0 , involves one or more visons on the plaquettes at each end (2 (b)). As the energy of this excitation depends on the string length and diverges for infinite vison separation, it follows that the latter are confined in an infinite system. case with Z 2 topological order, the GS manifold of H g2 now includes |η zz , φ 0 , where η zz indicates finite staggered pseudo-spins along z. Similarly, any open string defects are gapped. However, in contrast to the case for Z 2 topological order, the gap for these excitations remains finite for arbitrary string length, in the infinitesystem size. In the same limit, states with strings forming non-contractible loops become degenerate with |φ 0 for h z g 2/4 . The effect of H g4 is analogous to the case with Z 2 topological order. Consequently, the GS has the form in Eq. 18, with the replacement |η xx ; φ 0 → |η zz ; φ 0 , indicating a surviving pseudo-spin staggering. While the two-point correlator for the η z pseudo-spins is always finite, D x ij vanishes for infinite separation as in the Hubbard model with Z 2 × Z 2 topological order.
Since the energy cost of an open string remains finite even in an infinite-size system, the visons are deconfined. Similarly, as states with non-contractible loops of negative ρ γ bonds become degenerate with the GS, in the limit of infinite system-size, the latter acquires additional topological degeneracy. Consequently, this entails a sixteen-fold topological degeneracy on the torus, consistent with Z 2 × Z 2 topological order.
Conclusion and outlook. We have studied a bilayer adaptation of a QSL model on a square lattice with Kitaev-type interactions. We have shown that the low energy model exhibits an AFM Mott transition which corresponds to magnetic fragmentation in terms of the original DOF. We have corroborated these results by a perturbative calculation for the topological degeneracy, which is consistent with field theory analysis. The analysis we have presented here may be of particular value as a largely tractable yet highly non-trivial instance of magnetic fragmentation.
Interesting future directions include examining the role of fluctuations on the emergent symmetry which may reduce the ground state manifold via order by disorder [35][36][37]. Another direction is to generalize our mechanism for fractionalized Goldstone modes in bilayer systems to multilayer systems with larger emergent symmetries, such as SU (N ).
The study of moiré superlattices of QSLs is another intriguing direction [27,38]. Our work suggests that manifold new phenomena arising from a combination of emergent symmetry and strong interactions are awaiting discovery here, providing a new vista on strongly correlated magnetism.

I. MAJORANA FERMION REPRESENTATION
The Hamiltonian of the bilayer is where H K and H J are the intra-and inter-layer terms, respectively. H K has the form [? ], where Γ γ matrices obey the Clifford algebra, where σ(τ ) Pauli matrices act on spin(orbital) degrees of freedom (DOF). Here τ γ = τ x , τ y , τ z , I and γ = 1, 2, 3, 4 correspond to the 4 inequivalent bonds at each site. We introduce a Majorana fermion representation via and relabel b 5 i → c x i and c i → c y i for convenience to obtain where u ν,ij = ib γ νi b γ νj . The inter-layer interaction Hamiltonian, H J = J i Γ 5 1i Γ 5 2i can be similarly expressed as Note that the fermion Hilbert space per site is 2 3 = 8-dimensional, twice the size of the 4-dimensional physical Hilbert space expanded by Pauli matrices σ and τ . To faithfully represent the physical system with the Majorana fermion degrees of freedom, we have to enforce the following onsite constraint with a dispersion E K = ±4K sin 2 k x + cos 2 k y . E K has two inequivalent Dirac points at (0, ± π 2 ) as shown in Supplementary figure 1.

A. Symmetry of the full Hamiltonian
We first consider the symmetry of the single layer Hamiltonian (2). The (spinless) time reversal symmetry operator T = K is simply implemented by the complex conjugation K. In addition, we have a unitary U (1) ν Z C 2 symmetry, generated by U (1) "charge" rotations in each layer ν and the particle-hole symmetry in layer ν where we can choose an arbitrary (but fixed) bond direction 1 ≤ γ ≤ 4. Note that all these symmetries commute with the onsite constraint (7). Next we take the inter layer interaction term (6) into account to analyze the symmetry of the full Hamiltonian (1). Note that the particle-hole symmetry (10) in each layer does not commute with the onsite constraint (7), and hence ceases to be a symmetry of the full Hamiltonian (1). Instead, the following particle-hole symmetry for both layers remain as a symmetry of the full Hamiltonian. Moreover, there is a layer exchange symmetry generated by As a result, the symmetry of the full Hamiltonian is given by The symmetry of the full Hamiltonian can be conveniently understood in the basis of complex fermions f ν,j defined as In this basis, the U (1) symmetries for each layer can be recombined into "charge" U (1) c and pseudo-spin (i.e. layer) U (1) s rotations: where κ = −1, 1 for ν = 1, 2 respectively and ν = ν + 1(mod2). If we represent the pseudo-spin-1/2 complex fermions using the following matrix the pseudo-spin and charge rotations correspond to left and right rotations on the matrix F j : where σ L/R ≡ (X L/R , Y L/R , Z L/R ) stands for Pauli matrices multiplied to the left/right of matrix F j in (17). Now it is clear that the full symmetry group can be written as where O(2) = U (1) Z 2 is the rotation group for either the charge or pseudo-spin degrees of freedom. , where the layer index ν plays the role of the psuedo-spin. We have used j = (j x , j y ) ∈ Z 2 to label a site on the square lattice. Therefore the emergent symmetry group for a fixed flux configuration is given by It is most convenient to represent the symmetry in the u ν,ij ≡ u ij = ±1, ∀ ν = 1, 2 gauge, where the SO(4) symmetry is simply implemented by an SO(4) unitary rotation on the Majorana basis since the fixed-flux Hamiltonian in this gauge is written aŝ Meanwhile, under the spinless time reversal symmetry operation, the Majorana fermions (24) transform as The SO(4) symmetry can also be understood in the basis of complex fermions f ν,j defined in (14), where SU (2) s pesudo-spin rotation and SU (2) c charge rotation are realized as left and right SU (2) rotations on the matrix in (17): Physically, at an energy scale much smaller than the vison gap, vison excitations are suppressed and the low energy subspace has a fixed flux configuration. As proven by Lieb[? ], a uniform π flux per plaquette has the lowest energy and the low energy physics is described by Majoranas ψ j hopping in a uniform background of π flux per plaquette. Therefore, the effective (or emergent) symmetry of the system below the vison excitation gap is exactly given by the symmetry of a fixed flux configuration as we described above.

C. Remnant symmetry of the Néel order
In a Hubbard model on the square lattice, as the onsite repulsive Hubbard interaction strength J/K increases beyond a critical value J c /K, the system enters a Néel ordered phase with the long-range antiferromagnetic order. Now that we have obtained the symmetry of the full Hamiltonian (1) and of any fixed flux configuration, next we analyze the remnant (unbroken) symmetry of the Néel ordered ground state in our model (1).
The Néel order in the Hubbard model is characterized by the staggered magnetization as its order parameter, represented by the Néel vector: Note that the Néel vector always changes sign when transformed under the spinless time reversal symmetry T = K, therefore it always spontaneously breaks the antiunitary time reversal symmetry T . From (28) it is clear that the Néel vector is invariant under charge rotations M c ∈ SU (2) c in (14), therefore it always preserves the O(2) c subgroup of "charge" symmetries of the full Hamiltonian (1), or the SU (2) c subgroup of the emergent G = SO(4) × Z T 2 symmetry for a fixed flux configuration. In fact, independent of the Néel vector orientation, in a fixed flux configuration, the Néel order always spontaneously breaks the emergent symmetry group G = SO(4) × Z T 2 down to a subgroup where in the complex fermion basis (17) we have defined Herem ẑ × n is an in-plane unit vector perpendicular to the unit Néel vectorn. Below we analyze three cases with different orientations of the Néel vector n, where the symmetry group G in (20) of the full Hamiltonian is spontaneously broken down to different subgroups.
(1) If the Néel vector n is parallel to theẑ axis, the U (1) s symmetry is preserved while the layer exchange X is broken. As a result, the remnant symmetry group is The associated order parameter manifold is given by (2) If the Néel vector n lies in the plane perpendicular to theẑ axis, the U (1) s spin rotation is broken while a layer exchange symmetry X n defined by in the complex fermion basis (17) is preserved. As a result, the remnant symmetry group is where the preserved anti-unitary symmetry is defined as As a result, the associated order parameter manifold is (3) If the Néel vector has both in-plane andẑ components, both U (1) s rotation and layer exchange X are broken, leaving anti-unitary T n defined in (30) an unbroken symmetry. Therefore the remnant symmetry group is H(n z = 0, n x n y = 0) = O(2) c × Z Tn 2 (37) and the associated order parameter manifold is A summary of the long range orders and associated topological defects of the local order parameters is shown in Table I.

Néel vector
Unbroken subgroup H G/H Topological defects Gauge group  (25) of interacting fermions with bipartite hopping and particle-hole symmetry, a uniform π flux per plaquette has the lowest energy. As a result, in the low-energy manifold below the vison gap, we can consider a π-flux Hubbard model of fermions {f ν,j } on the square lattice. To understand the topological nature of the ground state and its excitations, we need to understand the wavefunctions of the bilayer spin liquid model.
In the low-energy manifold, since the flux (21) in every plaquette is fixed to F ν,ijkl ≡ −1, ∀ ν, we can choose a gauge so that where u ij is fixed as shown in Fig. 1. In the fixed gauge, let's denote the eigenstates for the Hubbard model (25) of fermions {f ν,j } (or equivalently of {ψ j }) as |ψ . Then the mean-field ansatz |M F of all fermions {c x,y ν,j , b γ ν,j } can be written as = u ij denotes the unique eigenstate of fermions {b γ ν,j } that satisfies the fixed gauge (39). The physical eigenstate of the bilayer spin liquid model (1) is then obtained by enforcing the onsite constraint (7) on the mean-field ansatz: via the projection operator 1+Dν,j 2 on every site. In the Néel order of the square-lattice Hubbard model, there is an off-diagonal long range order for the Néel vector (28) in the mean-field ansatz: However, this does not guarantee the same local order parameter to exist in the physical ground state (41). More precisely, while theẑ-component of the Néel vector where κ = −1, 1 for ν = 1, 2 respectively. n z j is a local operator and, the in-plane components do not commute with the onsite constraint (7). As a result, the off-diagonal long range order of local order parameters in the mean-field anstaz will translate into a string order parameter in the physical wavefunction: where (j 1 , j 2 , · · · , j M ) is a string of consecutive sites that connect j 1 = i and j M = j. It is straightforward to check that the string order parameter above commutes with the onsite constraint (7). Although the in-plane component n + of the Néel vector is not a local operator, there is a local order parameter for the in-plane Néel order: where i, j is a pair of nearest neighbor sites. And one can check that the in-plane Néel order in the mean-field ansatz |M F implies the following long range order of the physical wavefunction |QSL :

E. Topological orders in the Néel ordered phase
For a small value of J/K, the π-flux Hubbard model (25) has a gapless spectrum with massless Dirac fermion excitations. Each layer hosts one type of vison (i.e. π-flux) excitations, giving rise to a Z 2 × Z 2 gauge field. The low energy physics is therefore described by the two decoupled layers of gapless Z 2 spin liquids, each layer described by two branches of massless Dirac fermions coupled to a Z 2 gauge field.
In the large J/K case, a long-range Néel order develops in the ground states of the π-flux Hubbard model (25), where the Dirac fermions acquires a mass, giving rise to a gapped spectrum for all fermions. Different orientations of the Néel vector (28) not only give rise to distinct long range orders with different order parameters, but also distinct topological orders with different gauge groups. Below we analyze the 3 different orientations of the Néel vector summarized in Table I.

n ẑ
As previously discussed, when the Néel vector aligns alongẑ-axis, theẑ-component of the Néel vector is a local order parameter for the Ising-type Néel order, with a order parameter manifold of G/H = Z 2 . Since the fermions are all gapped, the ground state exhibits a Z 2 × Z 2 topological order, due to deconfined vison excitations from each layer.
The only gapless excitations in this case are the Goldstone modes of the Néel order, described by the spin flip operator in (44). However, the operator n + j ∼ (−1) jx+jy f † 2,j f 1,j is not a gauge invariant operator as it violates the onsite constraint (7). The gapless Goldstone modes of the Hubbard model lead to a power-law decaying correlation of the string order parameter (45) in this case.
If the system is gapped, it will naively be an Abelian Z 2 × Z 2 topological order described by a 4-component Chern-Simons theory with the following K matrix However, the nature of the gapless "Goldstone modes" here, plays a crucial role to interpret the physical ground state (41) after the projection 1+Dν,j 2 on each site. Notice, the two real Goldstone modes together form a complex bosonic mode which is gapless, with a linear dispersion. This mode carries a unit gauge charge of both Z 2 gauge fields, one from each layer, and is hence not a local excitation. Therefore, this ground state should be understood as an algebraic Z 2 × Z 2 spin liquid, where the bound state (50) (with quasiparticle vector (1, 1, −1, 1) T in the Abelian Chern-Simons theory) of two types of fermionic spinons f 2 (with quasiparticle vector (1, 1, 0, 0) T ) and f 1 (with quasiparticle vector (0, 0, −1, 1) T ) becomes gapless with an algebraic correlation in (45). The mode b obeys bosonic self statistics, but is really an anyon excitation with semionic mutual braiding with fermions f ν,j from either layer. Therefore the system exhibits fractionalized Goldstone modes that obey anyonic statistics. After incorporating the gapless anyon b in (50) into the low energy description, the effective field theory for this algebraic spin liquid writes where b denotes the critical anyon mode out of the Goldstone modes of the broken emergent symmetry. Typically, the robust presence of a gapless anyon excitation in a system requires a mechanism to protect it. In this case, it is the emergent pseudo-spin SU (2) s symmetry that protects the criticality of this anyon mode.

n ⊥ẑ
It's straightforward to show that a nonzero in-plane Néel component b = 0 corresponds to condensing the anyon characterized by vector (1, 1, −1, 1) in the Chern-Simons theory, and drives the system into a usual Z 2 topological order described by a 2-component Abelian Chern-Simons theory with K = 0 2 2 0 .
This can be understood by considering the topological defects of the local order parameter (46): since the order parameter manifold for an in-plane Néel order (n ⊥ẑ) is G/H = U (1) S 1 , the long-range order supports point defects (i.e. vortices) classified by an integer-valued vorticity π 1 (G/H) = Z. The vortex of local order parameter S + = S x + i S y in (46) is nothing but a π flux from either layer, characterized by quasiparticle vector (1, 0, 0, 0) or (0, 0, 1, 0) in the 4-component Abelian Chern-Simons theory. Since the vortices are logarithmically confined in two spatial dimensions, the π flux from either layer is confined with an in-plane Néel order. As a result, the ground state features a Z 2 topological order, whose vison excitation is the bound state of a π flux from both layers.
Next we describe the Goldstone modes in the in-plane Néel order. Without loss of generality, we assume the Néel vector to point alongx axis in the ground state, i.e. n x = 0 in the ground state of the gauge-fixed π-flux Hubbard model (25). In this case, we can write down the following local operator that creates the Goldstone mode n y : whereê γ is the unit lattice vector along bond γ direction. Due to the long-range string order (45) in the ground state, it is straightforward to show that When the Néel vector has both in-plane andẑ components, the ground state also exhibits a Z 2 topological order, for the same reason as discussed above. With an order parameter manifold of G/H = O(2), the topological defects in this case include both domain walls classified by π 1 (O(2)) = Z 2 , and vortices classified by π 1 (O(2)) = Z.

IV. EFFECTIVE η − ρ HAMILTONIAN IN THE LARGE-J LIMIT
In this section, we provide the details of the derivation for the effective Hamiltonian to fourth-order in |K/J| 1, which we label as an "η − ρ" model. For K = 0, the ground state (GS) is doubly degenerate, |↑ 1 ↑ 2 , |↓ 1 ↓ 2 at each site. These doublets can be represented by the pseudospins |⇑ , |⇓ with associated operators which obey an SU(2) algebra. Note that η z i is a dipolar operator whereas η x i and η y i , are quadrupole operators. Apart from the pseudospin doublet, there is an additional four-fold degeneracy per site, arising from the orbital DOF. We introduce an operator which projects out the states |↑ 1 ↓ 2 , |↓ 1 ↑ 2 at each site.
In the following, we derive the effective Hamiltonian perturbatively to fourth-order in K. We also note that all odd-ordered terms vanish since H K only connects states inside the GS manifold to states outside of it.
A. Second-order effective Hamiltonian Hg 2 The effective Hamiltonian to second order is and H J is the inter-layer coupling. In order to evaluate eq. 57, we introduce the operators and re-write H K = ij γ X ij + X ji . Since P 0 X 2 ij P 0 = 0, we obtain Furthermore, we use the identities to express where g 2 = K 2 /4J. We perform rotations by π about the z-axis on the η pseudo-spins at every other site to obtain the effective Hamiltonian at g 2 level Note that τ γ 1i τ γ 1j τ γ 2i τ γ 2j = ρ γ ij , where ρ γ ij = q γ i q γ j .
B. Fourth-order effective Hamiltonian Hg 4 The general expression to fourth order reads We focus on the first term, which, as described below, includes coupled pseudo-spin and orbital DOF around the p/p plaquettes. These play a critical role in lifting the extensive degeneracy of the GS of the leading H g2 , due to the orbital DOF, as implied by ρ γ ij = 1 ∀ γ, ij (Supplementary section VI). By contrast, the remaining terms in eq. 64 are either independent of the orbital DOF or involve two pairs of NN pseudo-spins coupled to ρ γ ij ρ γ jk , where γ = γ for the two connected bonds. These involve subleading corrections to the pseudo-spin configuration, but do not lift the GS degeneracy of H g2 . Consequently, we ignore these contributions and focus instead on the plaquette terms in the following. To illustrate, we consider a single p plaquette (Supplementary figure 1) in the first term as where . . . includes all remaining p/p plaquettes. We simplify the expressions by using the identities in eq. 61 and by rotating the pseudo-spins such that η ± i → −η ± i at every other site. Carrying out similar operations for the remaining plaquettes, we obtain the effective Hamiltonian at fourth order where Note that these are due to the terms in eq. 65.

V. THE EFFECTS OF τν,i AND W νp/p OPERATORS IN THE η − ρ BASIS
We first consider τ ν,i operating on any local orbital state determined by the set {q x i , q y i , q z i } subject to the constraint q x i q y i q z i = −1: In summary, τ x 1/2,i preserve q x i while changing the signs of q y/z i , and so on for all permutations. In terms of the ρ γ ij = q γ i q γ j bonds , the same operator changes the signs of ρ y ij and ρ z ik , where ij and ik are NN pairs. The flux operators are defined by W νp =σ z νk σ z νn τ x νi τ y νj τ x νk τ y νn (68) where the index convention is illustrated in Supplementary figure 1. These operators act both on the spin (σ), as well as on the orbital (τ ) DOF. The labels p, p label the plaquettes, as illustrated in Supplementary figure 1. By applying eq. 68 around a single p plaquette, we see that W 1p flips the q γ 's at four of the six sites of the unit cell. Consequently, six ρ γ ij bonds change signs, as shown in Supplementary figure 2 (a). The effect of W 1p is similar (Supplementary figure 2 (b)) although the end state differs in the q γ values at each of the six sites. Due to the presence of the σ z k and σ z n operators, both W 1p and W 1p also rotate the in-plane pseudo-spin components at sites k and n as η x,y k → −η x,y k and η x,y n → −η x,y n .

VI. EXACT DIAGONALIZATION OF Hg 2
The effective Hamiltonian couples the in-plane pseudo-spin components η x/y i with the orbital bonds ρ γ ij = q γ i q γ j . The latter are preserved and H g2 can be solved for the pseudo-spins in any given bond configuration which is consistent with the local constraint Here, we consider H g2 on a 4 × 4 torus with both h x/z = 0. We determine the GS energy via exact diagonalization in the pseudo-spin sector for various bond configurations. We find that the uniform ρ γ ij = 1 configuration has the lowest energy. In this sector, H g2 maps onto a 2D Heisenberg model with periodic boundary conditions. In contrast, all states with one or more strings of ρ γ ij = −1 bonds have higher energy, as shown in Supplementary figure 3.
Note that exceptions to these conclusions occur for h z > 0, h x = 0, or when the pseudo-spin spontaneously order along the z direction only, in the infinite-size system. As shown in Supplementary section VIII, a non-contractible string of ρ γ ij = −1 bonds becomes degenerate with the uniform ρ γ ij = 1 state.

VII. EXACT DIAGONALIZATION OF Hg 4
As discussed in Supplementary section VI, a state with uniform ρ γ ij = 1 bonds and zero net pseudo-spin, as in a 2D Heisenberg model, has the lowest energy at H g2 level. For h x = 0, H g2 commutes with the flux operators W νp/p . Therefore, the states where the ket is the GS at H g2 level, are degenerate with the latter. The products can include any combinations of W νp/p operators and this implies an extensive degeneracy of the GS. For h x > 0, H g2 only commutes with W ν,p W ν,p , where p/p are in the same unit cell, such that the degeneracy is partially lifted. We shall assume a small h x in the following and ignore this additional complication. For a uniform ρ γ ij = 1 configuration, this reduces to a 2D AFM Heisenberg model which has finite staggered pseudospins in the GS for any h x = 0. A non-contractible loop implies ABC forη z i = −η z i+N 1/2 . At a semi-classical level, the presence of the staggered field h x implies that such a twist costs an amount of energy which scales with N 2/1 . Moreover, this energy cost diverges for an infinite-size system. This also holds in the subsequent h x → 0 limit.
The same conclusion can be reached using a perturbative analysis. It is convenient to re-writẽ where When considering matrix elements, we shall label the expectation value in the GS ofH 0 as 00 . 01 indicates matrix elements between GS and an excited state with a NN pair of flipped pseudo-spins, and similarly 12 for the matrix elements between two excited states with 1 and 2 pairs of flipped pseudo-spins, respectively. The GS has staggered pseudo-spins to zeroth order in g 2 . To first order, the correction to the GS energy is where we explicitly indicated the number of contributing configurations and their weights. These are due to the Ising exchange terms along the (rotated) z axis, and amount to taking into account all of the bonds inH String andH Coupling . ρ stands for the two cases with either all bonds positive or negative on the string. The remaining contributions from H Background are independent of ρ. The corrections to second order are due exclusively to the (rotated) in-plane terms, which introduce pseudo-spin flips for NNs. Importantly, each flip carries a (1 + ρ γ ij ) factor as where ∆E 01 = −2h x is the energy cost of an excited state with a single pair of flipped pseudo-spins. Also note that V γ ⊥,ij = 2 −2 (1 + ρ γ ij ). To third order, we obtain These corrections are due to the energy cost introduced by V γ z terms in the presence of a pair of flipped pseudo-spins in the string and connecting bonds. To fourth order, the corrections read The first line includes all pseudo-spin flips around p/p plaquettes. The next three lines amount to flips on single bonds and on pairs of connected bonds, while the last line contains corrections similar to those at third order. Their effective contribution is The energy cost of the string to fourth order is It is apparent that this diverges as the system size, and thus n, grows to infinity. Note that this also implies that the cost of an open string of infinite extent likewise diverges. An open string of length (n − 2) is obtained from the non-contractible loop by changing the sign of two neighboring bonds. This is a local perturbation, which costs a finite amount of energy. The energy of the string therefore diverges as that of the parent state. As discussed in the main text, this implies that the visons are confined.
B. hz > 0, hx = 0 In contrast to the case with h x > 0, h z = 0, here the BCs associated with the loop do not imply a finite energy cost in a semi-classical approximation. This is due to the fact that the pseudo-spins are staggered along z, and η z does not couple to the bonds. We further illustrate by considering the Holstein-Primakoff representation for the pseudo-spins [? ].
Here, b i , b † i are bosons, n b,i = b † i b i , and η n b,i is the size of the pseudo-spin. The BCs for the pseudo-spins can be implemented by imposing ABCs for the bosons as b i = −b i+N2 . The remaining analysis follows the usual semi-classical approximation. In particular, as h z → 0, the zero-point energy of the GS involves a summation over the reciprocal unit cell as [? ].
where α k is the usual form factor for a square lattice. ABC amount to a trivial translation of the entire reciprocal unit cell in the extended BZ by G 1/2 /2N 1/2 , where G 2 is a reciprocal lattice vector. The correction to the classical GS energy is invariant. Consequently, at this level of approximation, the string is degenerate with the uniform ρ γ ij = 1 state.
We reach the same conclusion using a perturbative analysis. In contrast to the h x > 0 case, the first-order corrections are independent of ρ γ ij : since the Ising interaction along z does not couple to the bonds. The second-order contribution, due to pairs of NN pseudo-spin flips, reads Note the important difference w.r.t. the h x > 0 case, where each flip of NN pseudo-spins was multiplied by a (1 + ρ γ ij ) factor. Here, the corrections are independent of bond configuration. A similar conclusion holds at third order. At fourth order, we obtain the plaquette terms where we omitted contributions similar to those at third order, which likewise are independent of the bonds. The terms shown above involve products of ρ's around p/p plaquettes. Since each plaquette contains at least two bonds on the string, these are also trivial. However, similar corrections involving pseudo-spin flips around closed loops, which intersect the string an odd number of times, are non-trivial. It is straightforward to check that only non-contractible loops, which wrap around the torus at least once, can contribute. These occur at leading nth order and scale with a factor of g n 2 /h n−1 z . Consequently, each is exponentially suppressed as the system size, and therefore n, goes to infinity. The string becomes degenerate with the uniform ρ γ ij = 1 state in this limit. This also implies an upper bound on the energy cost of any open string of length n − 2. The latter can be obtained by changing the sign of any two NN bonds, which is a local perturbation, implying a finite energy cost. This also holds as the system size goes to infinity. As discussed in the main text, it follows that visons are deconfined.

IX. CONNECTION BETWEEN HUBBARD AND η − ρ MODELS
In this section, we show that the GSs of the Hubbard model and H η−ρ models are equivalent. We start with the Hubbard model (eq. 3 in the main text). In general, the GS for large J is characterized by a finite expectation value of the gauge-invariant correlators where are staggered, pseudo-spins along x, which are connected via a gauge string, and similarly for y components. We write the correlator in terms of Majorana fermions as We next consider a similar correlator in the π-flux GS of H η−ρ , (eq. 19 in the main text) The correlator commutes with W νp/p . It follows from the form of the GS of H η−ρ that D x ij is finite for h x > 0 and any i, j. We next perform a rotation on all of the spin (σ) DOF of the second layer as We introduce a Majorana representation for the rotated pseudo-spins as as well as for the string i j ∈Cij where we used ρ γ ij = ν τ γ νi τ γ νj and τ γ νi = b 4 νi b γ νi c x νi c y νi . Substituting these in eq. 102, we obtain D x ij = i j ∈(i,j) (−1) ix+iy+i x +i y (c x 1i c x 2i + c y 1i c y 2i )u 1i j u 2i j (c x 1j c x 2j + c y 1j c y 2j ) ΨGS .
This expression is formally identical to C x ij .

X. TOPOLOGICAL DEGENERACY OF THE GROUND STATE OF Hη−ρ
In this section we show that the GS of H η−ρ exhibits Z 2 topological order for h x > 0, h z = 0 and Z 2 ×Z 2 topological order for h z > 0, h x = 0.

A. Z2 topological order
For h x > 0, h z = 0, the GS of H η−ρij has the form where |η xx ; φ 0 is the GS of the 2D Heisenberg model for finite h x , while φ 0 is an orbital configuration with all q γ i = −1, which implies that it has uniform ρ γ ij = 1 bonds. Small corrections to this state occur because h x does not commute with the flux operators W νp/p . The projection of |η xx ; φ 0 onto the π-flux state in eq. 107 contains factors of the form W νp W νp , where p/p are in the same unit cell. When operating on |η xx ; φ 0 these preserve all pseudo-spin and ρ γ ij bonds. However, they change the sign of all pairs of NN pairs of q γ 's around the unit cell. To illustrate, the NN pair q z 1 = −1, q z j = −1 is mapped onto q z 1 = 1, q z j = 1, with ρ γ ij = 1 in both cases, and similarly for the remaining NN pairs. However, since these states are obtained from the φ 0 configuration via local changes in the set of q's, they are not topologically distinct. The W νp W νp operators are in this sense the analogs of plaquette operators in the toric code, while the signs of NN q γ 's can be used to define effective bonds. We construct topologically distinct configurations by operating with Here, the orbital operators τ x/y 1i are applied on every site along non-contractible loops C x/y , as illustrated in Supplementary figure 5. F x/y commutes with both H η−ρ and W ν,p/p , and has the effect of changing the signs of all NN q's along the loop. Note that τ x/y 2i would have the same effect. Any state thus obtained is degenerate with |Ψ GS . In order to capture the topology of the GS manifold, we introduce the Wilson loop operators where denotes every other site on C y/x (Supplementary figure 5). Since I x/y , H η−ρ = 0 and I x/y , W νp/p = 0, |Ψ GS can be labeled by λ x = λ y = 1. More generally, the states belong to four distinct topological sectors. The four-fold degeneracy is characteristic of Z 2 topological order. This procedure exhausts the degeneracy of the GS manifold, since any non-local changes in the pseudo-spin configurations cost finite energy for h x > 0. Similarly, changes in the ρ γ bonds are equivalent to the strings discussed in Supplementary section VIII and are likewise gapped for h x > 0.
B. Z2 × Z2 topological order The previous arguments for h x > 0, h z = 0 also hold for h z > 0, h x = 0. The Z 2 degeneracy due to configurations of the local orbital states q γ consistent with a given set of ρ γ bonds still holds. However, as shown in Supplementary section VIII, non-contractible loops of flipped (ρ γ ij → −ρ γ ij ) bonds becomes degenerate with |Ψ GS in the infinite-size limit. In order to account for this additional degeneracy, we can define the Wilson loop operators with eigenvalues λ A,x/y = ±1. These commute with H η−ρ and W ν,p/p . For any eigenvalues of I x/y (eq. 109), we can therefore distinguish four additional topologically distinct states The operators act along every other site on path C x/y to create non-contractible loops of flipped ρ γ bonds. The GS thus exhibits a sixteen-fold degeneracy on the torus, as consistent with a Z 2 × Z 2 topological order.