Competing Spin Liquid Phases in the S=$\frac{1}{2}$ Heisenberg Model on the Kagome Lattice

The properties of ground state of spin-$\frac{1}{2}$ kagome antiferromagnetic Heisenberg (KAFH) model have attracted considerable interest in the past few decades, and recent numerical simulations reported a spin liquid phase. The nature of the spin liquid phase remains unclear. For instance, the interplay between symmetries and $Z_2$ topological order leads to different types of $Z_2$ spin liquid phases. In this paper, we develop a numerical simulation method based on symmetric projected entangled-pair states (PEPS), which is generally applicable to strongly correlated model systems in two spatial dimensions. We then apply this method to study the nature of the ground state of the KAFH model. Our results are consistent with that the ground state is a $U(1)$ Dirac spin liquid rather than a $Z_2$ spin liquid.

Introduction -Quantum spin liquids (QSL) can be defined as zero-temperature quantum phases of spin systems in the absence of symmetry breaking. In the presence of translational symmetry, and if there are odd number of half-integer spins per unit cell, the Hastings-Oshikawa-Lieb-Schultz-Mattis theorem 1-3 indicates that a QSL is necessarily a nontrivial quantum phase beyond the Landau's paradigm. It has been pointed out that geometric frustration, strong quantum fluctuation and/or strong spin-orbit coupling may be helpful to realize a QSL in realistic spin systems 4,5 .
The nearest neighbour (NN) spin-1 2 kagome antiferromagnetic Heisenberg (KAFH) model is a spin model with a strong geometric frustration. Despite the simple form of the model Hamiltonian, the nature of the ground state of this model has been a long-standing puzzle and attracted considerable interest in the past few decades [6][7][8][9][10][11][12][13][14][15][16][17][18][19] . In particular, recently, a series of numerical simulations find this ground state to be a QSL. However, the nature of the QSL is still under debate. While numerical simulations based on density matrix renormalization group (DMRG) techniques 20,21 report evidences of a gapped Z 2 QSL 13-15 , state-of-the-art variational Monte Carlo simulations find the ground state to be a gapless U(1) Dirac spin liquid 16 . In addition, it is known that there are many different candidate Z 2 QSLs that may be realized in this model [22][23][24][25][26] due to the interplay between the symmetry and the Z 2 topological order -a phenomenon coined symmetry enriched topological(SET) phases. Consequently it is still unclear which one of these candidate Z 2 QSLs may be realized in this model.
The difficulty of the problem, to a large extent, is due to the lack of suitable theoretical/numerical techniques. In order to simulate even moderate system sizes of frustrated quantum spin systems like the KAFH model, one has to work with certain kinds of variational wavefunctions. The choice of variational wavefunctions often brings up the following dilemma: On the one hand, one would like to work with wavefunctions in specific universality classes so that the analytical understanding of the simulation is available. This is the philosophy behind most variational Monte Carlo simulations. On the other hand, in order to perform an unbiased simulation and to obtain accurate energetics, one hopes that the choice of the variational wavefunctions is as general as possible. For instance, DMRG simulations are based on the matrix product states (MPS) 27,28 -a quite general class of variational wavefunctions.
The problem is that the two desired features of the variational wavefunctions usually do not come together. For example, different candidate Z 2 QSLs in the KAFH model are characterized by the different symmetry fractionalization patterns on the anyon quasiparticle excitations [29][30][31][32] . It is highly nontrivial to extract such analytical understandings from a MPS 33 , although the DMRG simulations based on MPS provide very good energetics. At the same time, the variational Monte Carlo simulations based on the U (1) Dirac spin liquid state 16 , although having very clear analytical understanding, may be questioned about their generality.
It would be very interesting to develop new variational simulation schemes, hopefully capturing both desired features. This is indeed one of the motivations of an earlier piece of work by us, where we particularly pay attention to symmetric tensor-network wavefunctions 25,34 . In two spatial dimensions, we are focusing on the symmetric projected entangled pair states (PEPS) [35][36][37] , which are natural generalizations of MPS. It turns out that one can systematically classify general PEPS wavefunctions according to symmetry. Consequently one can obtain a finite number of classes of symmetric PEPS wavefunctions, and perform a variational simulation within each class separately. On the one hand, the analytical understanding of each class of symmetric PEPS wavefunction is available, which is related to, but not limited to, the symmetry fractionalization phenomenon. On the other hand, because the classification of symmetric PEPS is quite general, after variationally simulating different classes of symmetric PEPS, one is expected to have rather good energetics and nearly unbiased understanding of the quantum phase diagram.
In this work, we further develop the numerical simulation scheme based on symmetric PEPS, and apply it attempting to determine the nature of ground state of the KAFH model. The classification of symmetric PEPS 25,34 allows us to construct classes of generic PEPS wavefunctions for different Z 2 QSLs. We particularly focus on symmetric PEPS wavefunctions with bond dimension D = 6 and D = 7, and study the four Z 2 QSLs that can be realized by these bond dimensions. These four classes of symmetric PEPS correspond to Sachdev's Q 1 = Q 2 state, Q 1 = −Q 2 state and two other π-flux states. We perform variational simulations for the KAFH model for each class separately, and obtain four optimal energy densities. If the ground state is one of the four Z 2 QSLs, these optimal energy densities are expected to be significantly different, and the ground state is the Z 2 QSL with the lowest energy density.
However, surprisingly, we find that the optimal energy densities for both Q 1 = Q 2 state and Q 1 = −Q 2 state are nearly degenerate and comparable with the previously reported ground energy density of this model, while the two π-flux states have energy densities significantly higher. In fact, the most natural explanation for such a nearly degenerate energy density between the Q 1 = Q 2 state and Q 1 = −Q 2 state, without resorting to fine-tuning, is that the ground state is actually a U (1) Dirac QSL. This is because both the Q 1 = Q 2 state and Q 1 = −Q 2 state can be viewed as descendent states from the same parent U (1) Dirac QSL, and therefore can both be used to approximate the parent state. Consequently although we use Z 2 QSLs as trial wavefunctions, our results can be viewed as a supporting evidence of the U (1) Dirac QSL.
Spin-1 2 symmetric PEPS on kagome lattice -The kagome PEPS and various notations for sites and bonds are shown in Fig. 1(a). To construct a spin-1 2 kagome PEPS, we associated every site/bond of the kagome lattice with a site/bond tensor. As shown in Fig. 1(b), a site tensor is formed by a physical leg which support a physical spin-1 2 , and four virtual legs, while a bond tensor is formed by two virtual legs. Every leg is associated with a specific local Hilbert space, and a tensor can be viewed as a quantum state in the Hilbert space of the tensor product of all its leg Hilbert spaces. The physical wavefunction is obtained by contracting all connected virtual legs of site tensors and bond tensors.
The classification of symmetric spin liquid phases on the kagome PEPS was obtained in Ref. 25. Here, we briefly review the procedure and the result. The symmetry group of the spin-1 2 kagome system can be generated by translation symmetries T 1(2) , six-fold rotations about the center of the hexagon C 6 , mirror reflection σ along the dashed line in Fig.1(a), time-reversal symmetry T , and spin rotation symmetry U θ n . A global symmetry transformation g induces a gauge transformation W g (x, y, s, i) on all internal legs of tensors. Here (x, y, s) denotes the site position and i labels the leg, as shown in Fig.1(b). Different spin liquid phases are characterized by gauge inequivalent symmetry transform rules W g on internal legs of the tensor network.
In our case, physical legs are spin-1 2 's, while internal legs support virtual spin representations. We can label a virtual Hilbert space as 25, the global 2π spin rotation induces a special pure gauge transformation {J}, which leaves every single tensor invariant up to some phase factor. For instance, when V = 0 ⊕ 1 2 , J = diag [1, −1, −1] on every internal leg. {J} together with the identity action form a Z 2 invariant gauge group (IGG), which is related to the Z 2 toric code topological order.
The Z 2 IGG will enter tensor equations for symmetries and enrich the classification. Briefly speaking, W g , which is the symmetry action on internal legs, satisfies group multiplication rules up to an IGG element (either trivial or nontrivial) as well as a phase factor. Given the Z 2 IGG and global symmetries of the model, one obtains 32 inequivalent classes, which are characterized by five Z 2 indices: η 12 , η C6 , η σ and χ σ , χ T . Here η's label Z 2 IGG elements, which characterize symmetry fractionalizations of spinon e-particles, while χ's are phase factor ±1, which are related to "weak SPT" indices. For example, η 12 = I/J corresponds to zero-flux/π-flux spin liquids in the Schwinger boson language.
As listed in Appendix A, for all classes, we solve the symmetry transformation rules W g for arbitrary D by fixing gauge. The fact that tensors are invariant under symmetry actions on both physical legs and internal legs imposes constraints on the Hilbert space of local tensors. Tensors of different classes live in different constraint sub-Hilbert spaces. Here, we focus on two cases: D = 6 with virtual spins 0 ⊕ 1 2 ⊕ 1 and D = 7 with virtual spins 0 ⊕ 0 ⊕ 1 2 ⊕ 1. Only 4 of the 32 classes can be realized in these two cases, which are fully characterized by two indices η 12 and η C6 while other indices are fixed as Symmetric iPEPS algorithm -Given generic tensor wavefunctions for all classes, our goal is to find the optimal PEPS wavefunction for each class, which minimize h ij = Ψ|h ij |Ψ , where h ij is the local Hamiltonian acting on two neighbouring sites i, j. In NN KAFH, As shown in Fig.2(b), h ij is calculated by contracting all legs of a double layer PEPS. Notice that we have already absorbed bond tensors to neighbouring site tensors for convenience. The bond dimension of the double layer PEPS is D 2 , and it is generally impossible to get the exact result of tensor contraction. The key point of the iPEPS algorithm is to find a reasonable approximation for the environment tensor around site i, j. The algorithm is divided into two parts: optimization and measurement. In the following, we will describe these two parts separately.
Optimization. In this paper, we apply a modified simple update algorithm 42 to optimize the wavefunction. As shown in Fig.3, for the simple update method, the environment tensors of three sites are approximated by the direct product of matrices E. Here, all sites share the same matrix E due to lattice symmetries.
The algorithm to obtain E is described in the following: 1. First, we define the local wavefunction |ψ as contracting single layer site tensors in one unit cell with initial environment matrix E. As shown in Fig.3(a), we can decompose |ψ as |ψ = α |φ A α ⊗ |φ B α , where α labels virtual states living in the tensor product space of leg uv and vw.

Define
Then, M is a hermitian matrix, and can be decomposed as M = (X T ) † ·X T . The decomposition can be sped up a lot by implementing spin rotation symmetries. Then, |e α ≡ (X −1 ) αα |φ α form an orthonormal set. Similar analysis on B leads to an orthonormal set {|f β }. As shown in Fig.3(b), |ψ = αβ X αγ Y γβ |e α ⊗ |f β . We then perform singular value decomposition: X · Y = U ΛV , where Λ encodes the entanglement information of |ψ .
3. Update one environment matrix E to be E → ΛU X −1 , as shown in Fig.3(c), and then use spatial symmetry transformation rules in Appendix A to generate all environment matrices at different spatial positions.
4. Repeat the above procedure until Λ converges.
Given an arbitrary PEPS wavefunction |Ψ belonging to some spin liquid class, say, class A, we are able to efficiently measure the approximate"energy density" h ij su using the converged environment matrix E. We then implement standard minimization algorithm, for instance, the conjugate gradient method, to search for the optimal wavefunction in the constraint sub-Hilbert space of class A, which minimize h ij su . Notice that the major advantages of this simple-update algorithm are its stability and speed, although the approximation introduced by the direct-product-environment E is not well under control. In order to control the approximation in the environment tensor, other algorithms like full-update 43 need to be used, which we leave as a topic of future studies. Measurement. By implementing the optimization algorithm to all four classes, we obtain optimal wavefunctions for these classes. We then measure the energy density of each optimal wavefunction as accurately as we can. We mainly use variational Monte Carlo combined with tensor entanglement renormalization method (VMC-TERG) 44 to measure the energy density on a 192-site finite-size sample. (The energy density measurement based on iTEBD algorithm 45,46 is also per-  formed as a complementary check. See Appendix B for details.) VMC-TERG is a single-layer algorithm in which one has to approximate the tensor-contraction by keeping a finite bond-dimension D cut during the real-space tensor renormalization. Namely, a finite D cut would introduce approximation and a scaling analysis with respect to D cut is usually necessary. However, we would like to emphasize that despite having approximation for the tensor-contraction, for any given D cut , the energy measurement by VMC-TERG is variational. This sharp variational meaning of the VMC-TERG algorithm is one of its major advantage comparing with other algorithms.
Result -We perform the above algorithm to the four promising classes with virtual bond dimension D = 6 and D = 7. Results measured by VMC-TERG are presented in Fig.4. Energy densities of optimized wavefunctions are measured in the 8 × 8 × 3 kagome lattice, and the scaling over D cut is applied. We fit energy densities as power law functions of D cut : E ∼ 1 D α cut , 1 ≤ α ≤ 2, where fitting parameter α is chosen to have the best fitting quality 47 . Energy densities obtained from extrapolation to infinite D cut are presented in Table I. We warn the readers that this type of extrapolation-schemes is only empirically justified 44,48 .
As shown in Fig.4 and Table I, energy densities of the two zero-flux classes are significantly lower than those of the two π-flux classes, which indicates the ground state of NN KAFH should be a zero-flux spin liquid. However, these two zero-flux classes have degenerate optimal energy density for both D = 6 and D = 7 within error bar. It appears that our method fails to determine the correct class of the Z 2 spin liquid for the NN KAFH model.
In fact, the most natural way to interpret the energy degeneracy is that the ground state is actually a U (1) Dirac spin liquid 9,16,49 rather than a gapped Z 2 spin liquid. To justify this statement, we note that the U (1) Dirac spin liquid is the "parent class" of these two zeroflux Z 2 spin liquids. One way to see this is to go to the Abrikosov fermion language 38,50,51 , in which the U (1) Dirac spin liquid is described by gapless fermionic spinons coupling to the internal U (1) gauge field. By adding pairings of fermionic spinons, the U (1) gauge field will be Higgsed to Z 2 , leading to Z 2 spin liquids. Patterns of pairing are constrained by lattice symmetries, and different pairing patterns give different Z 2 spin liquids.
It turns out that these two zero-flux classes are exactly the neighboring phases of the same U (1) Dirac spin liquid 24 , while the two π-flux states are not neighboring phases of the U (1) Dirac spin liquid. 52 Consequently, any state belonging to the U (1) Dirac spin liquid can be approximated by wavefunctions of these two descendant zero-flux Z 2 spin liquids classes by turning on very small pairing. This would naturally lead to the optimal energy degeneracy obtained by the two zero-flux classes of symmetric PEPS, without having to resort to fine-tuning.
The optimal energy density measured here on the 192-site sample is comparable to the thermodynamiclimit energy density reported in a recent tensor-networkbased work Ref.53, and is slightly higher than the estimated thermodynamic-limit energy density obtained from DMRG 13 . We expect that the optimal variational energy can be further improved by implementing more accurate optimization methods, such as the fast full update algorithm 43 .
Discussion and Conclusion -We demonstrate a new variational numerical simulation scheme based on symmetric PEPS wavefunctions. Although we study the particular KAFH model in this paper, the classification and simulation of symmetric PEPS wavefunctions are generally applicable to other correlated quantum systems. The main advantage of this scheme is that two desired features of variational simulations are both realized. First, the systematic classifications and constructions of generic PEPS wavefunctions allow one to simulate the quantum phase without losing generality and obtain accurate energetics, which can be comparable with the energetics of other state-of-the-art variational methods. Second, despite being general, sharp analytical understandings for each class of symmetric PEPS wavefunctions are available.
In particular, we simulate four promising candidate spin liquids on the KAFH model. Two distinct zero-flux Z 2 QSLs give nearly degenerate optimal energy density which is comparable with the ground state energy density reported using other methods. The most natural explanation for this degeneracy is that the ground state is actually the U (1) Dirac spin liquid, which is the parent phase of both classes.
It is also informative to compare our simulation with previous parton-based variational studies on the two zero-flux Z 2 QSLs. Ref. 54 reported the variational Monte Carlo simulations based on Guzwiller-projected Schwinger-boson states, and it was found that the Q 1 = Q 2 state has an energy density significantly lower than that of the Q 1 = −Q 2 state. Although the Guzwillerprojected Schwinger-boson states are in the same universality classes as the two symmetric PEPS classes studied here, the energetics performance of the PEPS wavefunctions are much better. This can be intuitively understood as follows. The tunable variational parameters in parton-based wavefunctions quickly become long-ranged in the real space as one increases the number of parameters, which would not improve energetics -a shortrange property of the wavefunctions. However, the tunable variational parameters in symmetric PEPS wave- functions are directly enlarging the local Hilbert space for a local tensor, which can significantly improve energetics.
We thank Ling Wang, E. Miles Stoudenmire and Patrick Lee for helpful discussions, and Yin-Chen He for sharing his unpublished results on this model using DMRG techniques. The PEPS calculations were based on the ITensor library, http://itensor.org/. This study is supported by the Alfred P. Sloan fellowship and National Science Foundation under Grant No. DMR-1151440. We thank Boston College Research Service for providing the computing facilities where the numerical simulations were performed.