Competing Spin Liquid Phases in the S= Heisenberg Model on the Kagome Lattice
Abstract
The properties of ground state of spin 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 topological order leads to different types of spin liquid phases. In this paper, we develop a numerical simulation method based on symmetric projected entangledpair 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 Dirac spin liquid rather than a spin liquid.
Introduction  Quantum spin liquids (QSL) can be defined as zerotemperature quantum phases of spin systems in the absence of symmetry breaking. In the presence of translational symmetry, and if there are odd number of halfinteger spins per unit cell, the HastingsOshikawaLiebSchultzMattis theoremLieb et al. (1961); Oshikawa (2000); Hastings (2005) 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 spinorbit coupling may be helpful to realize a QSL in realistic spin systemsBalents (2010); WitczakKrempa et al. (2014).
The nearest neighbour (NN) spin 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 longstanding puzzle and attracted considerable interest in the past few decadesElser (1989); Marston and Zeng (1991); Nikolic and Senthil (2003); Ran et al. (2007); Singh and Huse (2007); Evenbly and Vidal (2010); Jiang et al. (2008a); Yan et al. (2011); Depenbrock et al. (2012); Jiang et al. (2012); Iqbal et al. (2013); He and Chen (2015); Mei et al. (2016); Liao et al. (2016). 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) techniquesWhite (1992, 1993) report evidences of a gapped QSLYan et al. (2011); Depenbrock et al. (2012); Jiang et al. (2012), stateoftheart variational Monte Carlo simulations find the ground state to be a gapless U(1) Dirac spin liquidIqbal et al. (2013). In addition, it is known that there are many different candidate QSLs that may be realized in this modelSachdev (1992); Wang and Vishwanath (2006); Lu et al. (2011); Jiang and Ran (2015); Qi and Cheng (2016) due to the interplay between the symmetry and the topological order — a phenomenon coined symmetry enriched topological(SET) phases. Consequently it is still unclear which one of these candidate 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)Östlund and Rommer (1995); Dukelsky et al. (1998) — 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 QSLs in the KAFH model are characterized by the different symmetry fractionalization patterns on the anyon quasiparticle excitationsWen (2002); Mesaros and Ran (2013); Essin and Hermele (2013); Barkeshli et al. (2014). It is highly nontrivial to extract such analytical understandings from a MPSZaletel et al. (2015), although the DMRG simulations based on MPS provide very good energetics. At the same time, the variational Monte Carlo simulations based on the Dirac spin liquid stateIqbal et al. (2013), 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 tensornetwork wavefunctionsJiang and Ran (2015); Mambrini et al. (2016). In two spatial dimensions, we are focusing on the symmetric projected entangled pair states (PEPS)Vidal (2008); Verstraete and Cirac (2004); Cirac and Verstraete (2009), 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 PEPSJiang and Ran (2015); Mambrini et al. (2016) allows us to construct classes of generic PEPS wavefunctions for different QSLs. We particularly focus on symmetric PEPS wavefunctions with bond dimension and , and study the four QSLs that can be realized by these bond dimensions. These four classes of symmetric PEPS correspond to Sachdev’s state, 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 QSLs, these optimal energy densities are expected to be significantly different, and the ground state is the QSL with the lowest energy density.
However, surprisingly, we find that the optimal energy densities for both state and 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 state and state, without resorting to finetuning, is that the ground state is actually a Dirac QSL. This is because both the state and state can be viewed as descendent states from the same parent Dirac QSL, and therefore can both be used to approximate the parent state. Consequently although we use QSLs as trial wavefunctions, our results can be viewed as a supporting evidence of the Dirac QSL.
Spin 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 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, 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.Jiang and Ran, 2015. Here, we briefly review the procedure and the result. The symmetry group of the spin kagome system can be generated by translation symmetries , sixfold rotations about the center of the hexagon , mirror reflection along the dashed line in Fig.1(a), timereversal symmetry , and spin rotation symmetry . A global symmetry transformation induces a gauge transformation on all internal legs of tensors. Here denotes the site position and labels the leg, as shown in Fig.1(b). Different spin liquid phases are characterized by gauge inequivalent symmetry transform rules on internal legs of the tensor network.
In our case, physical legs are spin’s, while internal legs support virtual spin representations. We can label a virtual Hilbert space as , where supports spin and denotes number of spin species, while is the “flavor” space. The dimension of is . As shown in Ref.Jiang and Ran, 2015, the global spin rotation induces a special pure gauge transformation , which leaves every single tensor invariant up to some phase factor. For instance, when , on every internal leg. together with the identity action form a invariant gauge group (IGG), which is related to the toric code topological order.
The IGG will enter tensor equations for symmetries and enrich the classification. Briefly speaking, , 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 IGG and global symmetries of the model, one obtains 32 inequivalent classes, which are characterized by five indices: and . Here ’s label IGG elements, which characterize symmetry fractionalizations of spinon particles, while ’s are phase factor , which are related to “weak SPT” indices. For example, corresponds to zeroflux/flux spin liquids in the Schwinger boson language.
As listed in Appendix A, for all classes, we solve the symmetry transformation rules for arbitrary 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 subHilbert spaces. Here, we focus on two cases: with virtual spins and with virtual spins . Only 4 of the 32 classes can be realized in these two cases, which are fully characterized by two indices and while other indices are fixed as and . These four classes happen to contain four types of NN RVB statesWang and Vishwanath (2006); Yang and Yao (2012); Schuch et al. (2012); Poilblanc et al. (2012); Poilblanc and Schuch (2013) as representative wavefunctions. For , in the absence of any symmetry, the Hilbert space of the site tensor would be dimensional. After implementing all symmetries, the constrained subHilbert space of a site tensor turns out to be only dimensional for all four classes, which significantly reduces the number of variational parameters. (For , the symmetryconstrained subHilbert space is dimensional for all four classes.)
Symmetric iPEPS algorithm  Given generic tensor wavefunctions for all classes, our goal is to find the optimal PEPS wavefunction for each class, which minimize , where is the local Hamiltonian acting on two neighbouring sites . In NN KAFH, .
As shown in Fig.2(b), 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 , 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 . 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 algorithmJiang et al. (2008b) 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 . Here, all sites share the same matrix due to lattice symmetries.
The algorithm to obtain is described in the following:

First, we define the local wavefunction as contracting single layer site tensors in one unit cell with initial environment matrix . As shown in Fig.3(a), we can decompose as , where labels virtual states living in the tensor product space of leg and .

Define . Then, is a hermitian matrix, and can be decomposed as . The decomposition can be sped up a lot by implementing spin rotation symmetries. Then, form an orthonormal set. Similar analysis on leads to an orthonormal set . As shown in Fig.3(b), . We then perform singular value decomposition: , where encodes the entanglement information of .

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” using the converged environment matrix . We then implement standard minimization algorithm, for instance, the conjugate gradient method, to search for the optimal wavefunction in the constraint subHilbert space of class A, which minimize . Notice that the major advantages of this simpleupdate algorithm are its stability and speed, although the approximation introduced by the directproductenvironment is not well under control. In order to control the approximation in the environment tensor, other algorithms like fullupdatePhien et al. (2015) 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 (VMCTERG)Wang et al. (2011) to measure the energy density on a 192site finitesize sample. (The energy density measurement based on iTEBD algorithmJordan et al. (2008); Orús and Vidal (2008) is also performed as a complementary check. See Appendix B for details.) VMCTERG is a singlelayer algorithm in which one has to approximate the tensorcontraction by keeping a finite bonddimension during the realspace tensor renormalization. Namely, a finite would introduce approximation and a scaling analysis with respect to is usually necessary. However, we would like to emphasize that despite having approximation for the tensorcontraction, for any given , the energy measurement by VMCTERG is variational. This sharp variational meaning of the VMCTERG algorithm is one of its major advantage comparing with other algorithms.
Classes  

zeroflux I  0.4354(2)  0.4366(3) 
zeroflux II  0.4351(6)  0.4365(5) 
flux I  0.4293(5)  0.4313(7) 
flux II  0.4296(8)  0.4227(4) 
Result  We perform the above algorithm to the four promising classes with virtual bond dimension and . Results measured by VMCTERG are presented in Fig.4. Energy densities of optimized wavefunctions are measured in the kagome lattice, and the scaling over is applied. We fit energy densities as power law functions of : , where fitting parameter is chosen to have the best fitting quality ^{1}^{1}1i.e., the minimal sum of square of the vertical distances from data points to the fitted curve. Energy densities obtained from extrapolation to infinite are presented in Table 1. We warn the readers that this type of extrapolationschemes is only empirically justifiedWang et al. (2011, 2013).
As shown in Fig.4 and Table 1, energy densities of the two zeroflux classes are significantly lower than those of the two flux classes, which indicates the ground state of NN KAFH should be a zeroflux spin liquid. However, these two zeroflux classes have degenerate optimal energy density for both and within error bar. It appears that our method fails to determine the correct class of the 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 Dirac spin liquidRan et al. (2007); Hermele et al. (2008); Iqbal et al. (2013) rather than a gapped spin liquid. To justify this statement, we note that the Dirac spin liquid is the “parent class” of these two zeroflux spin liquids. One way to see this is to go to the Abrikosov fermion languageYang and Yao (2012); Lu et al. (2014); Qi and Fu (2015), in which the Dirac spin liquid is described by gapless fermionic spinons coupling to the internal gauge field. By adding pairings of fermionic spinons, the gauge field will be Higgsed to , leading to spin liquids. Patterns of pairing are constrained by lattice symmetries, and different pairing patterns give different spin liquids.
It turns out that these two zeroflux classes are exactly the neighboring phases of the same Dirac spin liquidLu et al. (2011), while the two flux states are not neighboring phases of the Dirac spin liquid. ^{2}^{2}2The corresponding classes in the fermion language for zeroflux I/II classes are labeled as / respectively in Ref.Lu et al., 2011. Consequently, any state belonging to the Dirac spin liquid can be approximated by wavefunctions of these two descendant zeroflux spin liquids classes by turning on very small pairing. This would naturally lead to the optimal energy degeneracy obtained by the two zeroflux classes of symmetric PEPS, without having to resort to finetuning.
The optimal energy density measured here on the 192site sample is comparable to the thermodynamiclimit energy density reported in a recent tensornetworkbased work Ref.Xie et al., 2014, and is slightly higher than the estimated thermodynamiclimit energy density obtained from DMRGYan et al. (2011). We expect that the optimal variational energy can be further improved by implementing more accurate optimization methods, such as the fast full update algorithmPhien et al. (2015).
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 stateoftheart 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 zeroflux 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 Dirac spin liquid, which is the parent phase of both classes.
It is also informative to compare our simulation with previous partonbased variational studies on the two zeroflux QSLs. Ref.Tay and Motrunich, 2011 reported the variational Monte Carlo simulations based on Guzwillerprojected Schwingerboson states, and it was found that the state has an energy density significantly lower than that of the state. Although the Guzwillerprojected Schwingerboson 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 partonbased wavefunctions quickly become longranged 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 wavefunctions 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 YinChen 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. DMR1151440. We thank Boston College Research Service for providing the computing facilities where the numerical simulations were performed.
References
 Lieb et al. (1961) E. Lieb, T. Schultz, and D. Mattis, Annals of Physics 16, 407 (1961).
 Oshikawa (2000) M. Oshikawa, Phys. Rev. Lett. 84, 1535 (2000).
 Hastings (2005) M. B. Hastings, EPL (Europhysics Letters) 70, 824 (2005).
 Balents (2010) L. Balents, Nature 464, 199 (2010).
 WitczakKrempa et al. (2014) W. WitczakKrempa, G. Chen, Y. B. Kim, and L. Balents, Annual Review of Condensed Matter Physics 5, 57 (2014).
 Elser (1989) V. Elser, Phys. Rev. Lett. 62, 2405 (1989).
 Marston and Zeng (1991) J. Marston and C. Zeng, Journal of Applied Physics 69, 5962 (1991).
 Nikolic and Senthil (2003) P. Nikolic and T. Senthil, Phys. Rev. B 68, 214415 (2003).
 Ran et al. (2007) Y. Ran, M. Hermele, P. A. Lee, and X.G. Wen, Phys. Rev. Lett. 98, 117205 (2007).
 Singh and Huse (2007) R. R. P. Singh and D. A. Huse, Phys. Rev. B 76, 180407 (2007).
 Evenbly and Vidal (2010) G. Evenbly and G. Vidal, Phys. Rev. Lett. 104, 187203 (2010).
 Jiang et al. (2008a) H. C. Jiang, Z. Y. Weng, and D. N. Sheng, Phys. Rev. Lett. 101, 117203 (2008a).
 Yan et al. (2011) S. Yan, D. A. Huse, and S. R. White, Science 332, 1173 (2011).
 Depenbrock et al. (2012) S. Depenbrock, I. P. McCulloch, and U. Schollwöck, Phys. Rev. Lett. 109, 067201 (2012).
 Jiang et al. (2012) H.C. Jiang, Z. Wang, and L. Balents, Nature Physics 8, 902 (2012).
 Iqbal et al. (2013) Y. Iqbal, F. Becca, S. Sorella, and D. Poilblanc, Phys. Rev. B 87, 060405 (2013).
 He and Chen (2015) Y.C. He and Y. Chen, Phys. Rev. Lett. 114, 037201 (2015).
 Mei et al. (2016) J.W. Mei, J.Y. Chen, H. He, and X.G. Wen, arXiv preprint arXiv:1606.09639 (2016).
 Liao et al. (2016) H. J. Liao, Z. Y. Xie, J. Chen, Z. Y. Liu, H. D. Xie, R. Z. Huang, B. Normand, and T. Xiang, arXiv preprint arXiv:1610.04727 (2016).
 White (1992) S. R. White, Phys. Rev. Lett. 69, 2863 (1992).
 White (1993) S. R. White, Phys. Rev. B 48, 10345 (1993).
 Sachdev (1992) S. Sachdev, Phys. Rev. B 45, 12377 (1992).
 Wang and Vishwanath (2006) F. Wang and A. Vishwanath, Phys. Rev. B 74, 174423 (2006).
 Lu et al. (2011) Y.M. Lu, Y. Ran, and P. A. Lee, Phys. Rev. B 83, 224413 (2011).
 Jiang and Ran (2015) S. Jiang and Y. Ran, Phys. Rev. B 92, 104414 (2015).
 Qi and Cheng (2016) Y. Qi and M. Cheng, arXiv preprint arXiv:1606.04544 (2016).
 Östlund and Rommer (1995) S. Östlund and S. Rommer, Phys. Rev. Lett. 75, 3537 (1995).
 Dukelsky et al. (1998) J. Dukelsky, M. A. MartínDelgado, T. Nishino, and G. Sierra, EPL (Europhysics Letters) 43, 457 (1998).
 Wen (2002) X.G. Wen, Phys. Rev. B 65, 165113 (2002).
 Mesaros and Ran (2013) A. Mesaros and Y. Ran, Phys. Rev. B 87, 155115 (2013).
 Essin and Hermele (2013) A. M. Essin and M. Hermele, Phys. Rev. B 87, 104406 (2013).
 Barkeshli et al. (2014) M. Barkeshli, P. Bonderson, M. Cheng, and Z. Wang, arXiv preprint arXiv:1410.4540 (2014).
 Zaletel et al. (2015) M. Zaletel, Y.M. Lu, and A. Vishwanath, arXiv preprint arXiv:1501.01395 (2015).
 Mambrini et al. (2016) M. Mambrini, R. Orus, and D. Poilblanc, arXiv preprint arXiv:1608.06003 (2016).
 Vidal (2008) G. Vidal, Phys. Rev. Lett. 101, 110501 (2008).
 Verstraete and Cirac (2004) F. Verstraete and J. I. Cirac, arXiv preprint condmat/0407066 (2004).
 Cirac and Verstraete (2009) J. I. Cirac and F. Verstraete, Journal of Physics A: Mathematical and Theoretical 42, 504004 (2009).
 Yang and Yao (2012) F. Yang and H. Yao, Phys. Rev. Lett. 109, 147209 (2012).
 Schuch et al. (2012) N. Schuch, D. Poilblanc, J. I. Cirac, and D. PerezGarcia, Physical Review B 86, 115108 (2012).
 Poilblanc et al. (2012) D. Poilblanc, N. Schuch, D. PérezGarcía, and J. I. Cirac, Phys. Rev. B 86, 014404 (2012).
 Poilblanc and Schuch (2013) D. Poilblanc and N. Schuch, Phys. Rev. B 87, 140407 (2013).
 Jiang et al. (2008b) H. C. Jiang, Z. Y. Weng, and T. Xiang, Phys. Rev. Lett. 101, 090603 (2008b).
 Phien et al. (2015) H. N. Phien, J. A. Bengua, H. D. Tuan, P. Corboz, and R. Orús, Phys. Rev. B 92, 035142 (2015).
 Wang et al. (2011) L. Wang, I. Pižorn, and F. Verstraete, Phys. Rev. B 83, 134421 (2011).
 Jordan et al. (2008) J. Jordan, R. Orús, G. Vidal, F. Verstraete, and J. I. Cirac, Phys. Rev. Lett. 101, 250602 (2008).
 Orús and Vidal (2008) R. Orús and G. Vidal, Phys. Rev. B 78, 155117 (2008).
 (47) I.e., the minimal sum of square of the vertical distances from data points to the fitted curve.
 Wang et al. (2013) L. Wang, D. Poilblanc, Z.C. Gu, X.G. Wen, and F. Verstraete, Phys. Rev. Lett. 111, 037202 (2013).
 Hermele et al. (2008) M. Hermele, Y. Ran, P. A. Lee, and X.G. Wen, Phys. Rev. B 77, 224413 (2008).
 Lu et al. (2014) Y.M. Lu, G. Y. Cho, and A. Vishwanath, arXiv preprint arXiv:1403.0575 (2014).
 Qi and Fu (2015) Y. Qi and L. Fu, Phys. Rev. B 91, 100401 (2015).
 (52) The corresponding classes in the fermion language for zeroflux I/II classes are labeled as / respectively in Ref.\[email protected]Lu:2011p224413.
 Xie et al. (2014) Z. Y. Xie, J. Chen, J. F. Yu, X. Kong, B. Normand, and T. Xiang, Phys. Rev. X 4, 011025 (2014).
 Tay and Motrunich (2011) T. Tay and O. I. Motrunich, Phys. Rev. B 84, 020404 (2011).
Appendix A Symmetry transformation rules and constrained subHilbert spaces for the symmetric kagome PEPS classes
a.1 The symmetry group for KAFH
As shown in Fig.(1), we label the three lattice sites in each unit cell with sublattice index . Further, we specify the virtual index of a given site. We choose Bravais unit vector as and . Thus, we are able to specify the virtual degrees of freedom of site tensors as . The symmetry group of such a twodimensional kagome lattice is generated by the following operations
(1) 
together with time reversal . Here,
The symmetry group of a kagome lattice is defined by the following algebraic relations between its generators:
(2) 
where stands for the identity element in the symmetry group.
Further, consider system with spin rotation symmetry operator , which means spin rotation about axis through angle . We mainly consider halfinteger spins ( symmetry) in this paper. The spin rotation symmetry commutes with all lattice symmetries as well as time reversal symmetry:
(3) 
a.2 Symmetry transformation rules on internal legs
There are symmetric PEPS classes, labeled by five indices: , where and . We choose to be the direct sum of for the integer spin subspace and for the halfinteger spin subspace by fixing gauge.
As shown in Ref.Jiang and Ran, 2015, symmetry transformation rules on internal legs can be represented as:
(5) 
For the rotation transformation , we have
(6) 
For the reflection transformation , we have
(7) 
And for the time reversal transformation , we have
(8) 
where
(9) 
Here is dimension of the extra degeneracy associated with spin. Namely, the total degeneracy for spin living on one virtual leg equals . We have the virtual bond dimension
(10) 
And, is a dimensional antisymmetric matrix.
For ’s, we have
(11) 
where
(12) 
Appendix B Optimal energies measured by iTEBD
We use the iTEBD methodJordan et al. (2008) to measure the energy densities of optimal wavefunctions belonging to two zero flux classes with . As shown in Fig. 5, the energy densities are still fluctuating up to , and it is hard to see the trend for larger . However, these results are in agreement with the VMCTERG result if one intuitively treats the fluctuation ranges as error bars.