The spin1/2 Kagome XXZ model in a field:
competition
between lattice nematic and solid orders
Abstract
We study numerically the spin1/2 XXZ model in a field on an infinite Kagome lattice. We use different algorithms based on infinite Projected Entangled Pair States (iPEPS) for this, namely: (i) an approach with simplex tensors and 9site unit cell, and (ii) an approach based on coarsegraining three spins in the Kagome lattice and mapping it to a squarelattice model with local and nearestneighbor interactions, with usual PEPS tensors, 6 and 12site unit cells. Similarly to our previous calculation at the SU(2)symmetric point (Heisenberg Hamiltonian), for any anisotropy from the Ising limit to the XY limit, we also observe the emergence of magnetization plateaus as a function of the magnetic field, at using 6 9 and 12site PEPS unit cells, and at and using a 9site PEPS unit cell, the later setup being able to accommodate solid order. We also find that, at , (lattice) nematic and VBCorder states are degenerate within the accuracy of the 9site simplexmethod, for all anisotropy. The 6 and 12site coarsegrained PEPS methods produce almostdegenerate nematic and VBCSolid orders. We also find that, within our accuracy, the 6site coarsegrained PEPS method gives slightly lower energies, which can be explained by the larger amount of entanglement this approach can handle, even in the cases where the PEPS unitcell is not commensurate with the expected ground state unit cell. Furthermore, we do not observe chiral spin liquid behaviors at and close to the XY point, as has been recently proposed. Our results are the first tensor network investigations of the XXZ model in a field, and reveal the subtle competition between nearby magnetic orders in numerical simulations of frustrated quantum antiferromagnets, as well as the delicate interplay between energy optimization and symmetry in tensor network numerical simulations.
pacs:
03.65.Ud, 71.27.+a, 75.10.bÊ
I Introduction
Frustrated quantum antiferromagnets pose some of the hardest numerical challenges in condensed matter physics. This is so because, given frustration, quantum Monte Carlo is hampered by the infamous sign problem. In this context, the Kagome lattice is one of the most famous examples of frustrated lattices given the presence of triangles, where competing antiferromagnetic interactions between nearestneighbour spins lead to many eigenstates with very different structure but with very similar and low energy, all of them competing to be the true ground state of the system. ÊThis makes the Kagome lattice an ideal experimental and theoretical setup to have very large quantum fluctuations at low temperatures, specially for spin1/2 systems with Heisenberglike interactions, but also for larger spins fru ().
An important model in this context is the the spin1/2 XXZ model in a field on the Kagome lattice. In such model, an anisotropy angle is introduced that distinguishes the interactions between nearest neighbours from the and ones. For generic values of anisotropy and field the model has symmetry, and one recovers the full symmetry of the Heisenberg model at the isotropic point when the field is absent. The anisotropy allows us to study the crossover from the XY model () to the Ising model (). As such, the model exhibits frustration, and in the presence of a magnetic field it gives rise to quantized plateaus in the longitudinal (i.e. along the field) magnetization , being the most prominent at , Êbut also possible at and plateaus (). The existence, nature and properties of such plateaus have been a hot topic of discussion, with some works proposing that in the XY point (i.e. for purely and interactions) the plateau may give rise to a chiral topological quantum spin liquid chirxy (). In fact, even the existence of the plateau at the XY point is also a matter of controversy. Moreover, it is not clear if the plateau undergoes a phase transition as one tunes the anisotropy between the XY and the Ising points.
In this paper we study the zerotemperature phase diagram of the above model with several methods based on infinite Projected Entangled Pair States (iPEPS) ipeps (), similarly to previous studies of the Heisenberg point () heisipeps (). In particular, we use:

An approach with simplex tensors and 9site unit cell pess ().

An approach based on a coarsegraining three spins in the Kagome lattice and mapping it to a squarelattice model with local and nearestneighbour interactions, with usual PEPS tensors, 6 and 12site unit cells heisipeps ().
As we shall see, the use of different unit cells produces different results in the phase diagrams, depending on whether the magnetic structure of the phase at some plateaus is commensurate with the used cell. Moreover, we shall also see that methods for which the structure of some plateaus is not commensurate may produce, however, slightly lower variational energies because they are able to handle a larger amount of entanglement. In particular, we do a detailed analysis of this situation for the plateau. Our different methods reveal a subtle competition between symmetry and energy in the model: while the 9site unit cell seems bettersuited for the expected symmetry properties of the ground state in some parameter regimes such as the and plateaus, the approach based on the 6site unit cell, even if not favourable from the symmetry perspective, turns out to be able to handle more entanglement in the wavefunction and may thus produce slightly better variational energies ^{1}^{1}1In fact, the property that symmetric tensor networks do not necessarily produce better variational energies is wellknown for tensor network methodspess ().. We also observe a competition between different types of order: (i) Valence Bond Crystal solid (VBCSolid) order breaking lattice translation symmetry and (ii) nematic order only breaking rotation symmetry heisipeps0 (). More specifically, we see that the 9site simplexmethod produces degenerate nematic and VBCSolid states ^{2}^{2}2Notice that in Ref.heisipeps () this order was called “Solid 1”, in order to distinguish it from other orders with a unit cell at other values of the magnetization., within our accuracy, for all values of the anisotropy up to the Ising point, whereas the 6 and 12site coarsegrained PEPS methods produce almost degenerate nematic and VBCSolid states, Êwith a slightly lower (but almost degenerate) energy than the 9site approach (see Fig.1 for schematic diagrams of these orders). We also do not find any of the characteristic signatures of a chiral spin liquid phase in the XY point chirxy (). This paper is also the first tensor network study of the XXZ model in a field on a Kagome lattice.
Ii Model and methods
Ê
ii.1 The Hamiltonian
Ê Here we consider the antiferromagnetic XXZ model for spin1/2 in the presence of an external magnetic field along the direction on the Kagome lattice. Its Hamiltonian is given by
(1)  
where is the spin1/2 th operator at site , is the magnetic field, is the anisotropy angle, and the interaction is for nearestneighbour spins. For generic values of and this Hamiltonian is symmetric under spin rotations around the axis. Different interesting points can be accessed by tuning parameter , namely:

The Heisenberg point for .

The XY point for .

The Ising point for .
ÊAfter doing an overall analysis of the phase diagram, we shall focus on the incompressible phase at reduced magnetization where, in the classical Ising limit, a macroscopic degeneracy occurs between all configurations with two spin up and one spin down (for ) in every triangle wannier (). When a small XY (ferromagnetic or antiferromagnetic) anisotropy is added, i.e. , Eq. (1) can be mapped to an effective quantum dimer model on the dual hexagonal lattice cabra2005 () whose groundstate was argued to be a VBCSolid moessner2001 (). This result is supported by direct Quantum Monte Carlo simulations on the (nonfrustrated) hardcore boson model isakov2006 () equivalent to the ferromagnetic XXZ chain given by . Here we wish to investigate what happens away from the Ising point in this case.
ii.2 Numerical approaches
Following our approach in Ref.heisipeps (), we use imaginarytime evolution in order to obtain approximations of the ground state of the system in the thermodynamic limit. To implement this, here we use several approaches based on infinitePEPS (or iPEPS), also similarly to what we did in Ref.heisipeps (). On the one hand, we use an approach based on introducing simplex tensors, also called “Projected Entangled Simplex States” (PESS) pess (), where we use a 9site unit cell and the socalled “simple update” su (). On the other hand, we also use an approach where three spins in the Kagome lattice are mapped to a single coarsegrained 8dimensional spin on the square lattice while keeping all interactions local, see Ref.heisipeps (), and then use a standard PEPS algorithm for the square lattice with 2site and 4site squarelattice unit cells, amounting to 6site or 12site unit cells in the original Kagome lattice. We have implemented this second approach both with the “simple update” and also the socalled “fast full update” ffu (), but we saw no significant difference in our results with these two update schemes, so we stick to the simple update because it is more efficient. The approach based on the coarsegrained Kagome lattice is less efficient but, however, has more variational parameters, Êwhich allows to handle a larger amount of entanglement in the ansatz wavefunction. For more specific details about these approaches, we refer the reader to Ref.heisipeps ().
In all cases, our refining parameter is the PEPS bond dimension , which controls the amount of entanglement in the tensor network, and which we consider up to . Depending on the choice of algorithm, translation symmetry may be broken in different ways, and also the number of variational parameters may be different. For instance, for the PESS approach with a 9site unit cell, translation symmetry may be broken in, e.g., a superstructure. On the other hand, the coarsegrained lattice approach with 6site and 12site unit cells may break translation invariance in other ways, such as a structure. Concerning the number of variational parameters, the coarsegrained approach easily has many more than the PESS approach (e.g. a factor of 4 for and the 6site unit cell), which means that the calculations need more CPU time, but they Êare also able to handle more entanglement in the wavefunction. As in Ref.heisipeps (), the characterization of the phases is possible by checking the local magnetizations as a function of the field, as well as the link energy terms , with as in Eq.(1). Expectation values and environment contractions are computed using Corner Transfer Matrices, see, e.g., Refs.ctm (); tn (). Moreover, in order to check the possible chirality of the wavefunction, we compute the entanglement spectra of the state wrapped around half an infinite cylinder, by following the procedure in, e.g, Ref.chirpeps ().
Iii Results
Ê
In Fig.2 and Fig.3 we show the different plateaus that arise as a function of the external magnetic field and the anisotropy angle , as computed with the 9site PESS and the 6site coarsegrained PEPS. We find that the 9site PESS reproduces different plateaus at and , whereas the 6site coarsegrained PEPS reproduces only the , due to the incommensurability of this unit cell with the expected magnetic structure of the other plateaus. Concerning the ground state energy, the approach based on the 6site cell seems to produce slightly lower values than the 9site cell. For instance, at the plateau region for we find that for , the ground state energies per site are (where in fact there is no plateau) and . This slight energetic advantage may be due to the larger number of variational parameters involved in the 6site wavefunction for a given . However, variations beyond the third significant digit should be taken with care, since the relative difference is still very small specially in the context of the several approximations involved in both algorithms, and which we control to the best of our possibilities.
For concreteness, in Fig.4 we show the magnetization curve, i.e. as a function of the external field obtained with the 6site algorithm, for different values of the anisotropy angle . We observe the presence of a very prominent magnetization plateau at , for all the scanned values of , including the XY, Heisenberg and Ising points. Away from the Ising point we also see a jump towards the saturation value Êfor several values of theta jump () – a more finegrained analysis (not displayed) shows that this is indeed the case –. The width of the plateau increases as we tune from the XY towards the Ising point, see Fig.5. Results for the 12site algorithm are essentially identical. Notice also that the 9site PESS algorithm also finds a plateau at the XY point. ^{3}^{3}3Note that it was suggested in Ref.cabra2005 () that the plateau width vanishes at .
In Fig.6 we show the ground state energy per site of in the plateau phase as a function of the anisotropy angle , as computed with our three methods. We see that the 6site approach is the one that produces a lower variational energy. The difference with the 9site PESS is around the third digit, roughly 1%, whereas the difference with the 12site PEPS is even smaller, around 0.1%. As said before, these small differences should be taken with caution, given the different approximations involved in the methods used. The fact that the 9site PESS approach produces higher energy may be a consequence of the smaller number of variational parameters in this ansatz. The quasidegeneracy between the 6site and the 12site PEPS energies shows that a doubling of the unit cell is preferred over a quadrupling. We have also observed that the 9site PESS produces degenerate nematic and VBCSolid states (within our energy resolution), similar to what we found in Ref.heisipeps () at the Heisenberg point. This degeneracy is very slightly lifted in the 6 and 12site PEPS approaches, in favor of the nematic order for XY anisotropy, in the vicinity of the Heisenberg point (with a typically relative energy difference of at the Heisenberg point) and for a small Ising anisotropy, and in favor of VBCSolid order close to the Ising point. This is true even for the largest bond dimensions that we could reach. It may be possible that the energy difference between these states becomes narrower in the largesize limit, thus explaining this observation. The stability of one phase w.r.t. the other is assessed when systematically obtaining the same one in different runs of our algorithms with different initial conditions.
We next perform an analysis of the plateau Êwith the different algorithms. In Fig.7 we show diagrammatically the longitudinal magnetizations at every site as well as the link energy terms at every different link, for the 6site PEPS (left) and the 9site PESS (right) calculations, for points in the plateau and four different anisotropy angles . In our simulations we have seen that, for the 6site PEPS approach, for the phase corresponds, according to the classification in Ref.heisipeps (), to a nematic state, where rotation symmetry is broken down to (3fold degenerate). For we obtain a VBCSolid phase (again according to the classification in Ref.heisipeps ()), where translation symmetry is broken down to a cell, and the discrete symmetries (lattice 6fold rotation) and (reflexion about one axis) are fully broken. Note that the Heisenberg point lies then within the nematic region. Comparison with the results obtained from the 9site PESS approach is useful. In the right panel of Fig. 7 we show different orders obtained in this way. In practice we have seen that, for every possible , nematic and VBCsolid orders are degenerate in energy within our accuracy. The structures in the figure are examples of what we obtain in some runs of the algorithm. By running the algorithm with different random initial conditions, we get either a nematic or VBCSolid state, both with the same energy. Such a degeneracy between nematic and VBCSolid seems to be lifted slightly in the 6 and 12site PEPS approaches.
Finally, in order to see possible signs of chirality, we have computed the entanglement spectrum of the PEPS obtained from the 6site PESS, for half an infinite cylinder of width 12, within the plateau at the XY point. Such state has been conjectured to be chiral chirxy (). In our simulation, however, the entanglement spectrum is perfectly symmetric under timereversal and therefore not chiral, see Fig.8 ^{4}^{4}4ÊAs a word of caution, it is wellknown that finding a gapped chiral topological state with PEPS may actually be quite difficult because of fundamental constraints of the PEPS tensor network, and therefore it is also difficult to get an intuition about from a numerical optimization.
Iv Conclusions
Ê Here we have studied the zerotemperature phase diagram of the Kagome XXZ model in a field and in the thermodynamic limit using different tensor network approaches. We find different plateau structures for the longitudinal magnetization as a function of the field and the anisotropy, depending on the unit cell being used. We also find that a 6site coarsegrained PEPS approach seems to produce slightly lower ground state energies, probably because the larger number of available parameters which allows for more entanglement in the ansatz wavefunction, and in spite of the fact that the unit cell cannot accommodate the superstructures found in other recent work huerga2016 (). Such small energy differences should be considered with caution, though, given the different approximation schemes involved in each method. From our results, we cannot discern if our 6site ansatz, which contains more adjustable parameters, slightly bias the results towards a 6site superstructure, or if this superstructure is indeed a genuine feature of the system. In any case, our results show a tight competition between nematic and VBCSolid orders, which we find for all values of the anisotropy, Êas well as the delicate interplay between the implementation of (lattice) symmetries and the optimization of the energy in tensor network algorithms.
Acknowledgements.
A. K. and R. O. acknowledge financial support from JGU, as well as the University of Toulouse and the CCBPP, where part of this work was performed. Discussions with Y. Iqbal and M. Ziegler are also acknowledged.References
 (1) S. Yan, D. A. Huse, and S. R. White, Science, 332 (6034):11731176 (2011); S. Depenbrock, I. P. McCulloch, and U. Schollwöck, Phys. Rev. Lett. 109, 067201 (2012); Y. Iqbal, F. Becca, S. Sorella, and D. Poilblanc, Phys. Rev. B 87, 060405(R) (2013); D. Poilblanc and N. Schuch, Phys. Rev. B 87, 140407 (2013); D. Poilblanc, N. Schuch, D. PérezGarcía, and J. I. Cirac, Phys. Rev. B 86, 014404 (2012); K. Hida, J. Phys. Soc. Jpn, 70:3673 (2001); A. Honecker, J. Schulenburg, and J. Richter, J. Phys.: Condens. Matter 16, S749S758 (2004); Y. Okamoto, M. Tokunaga, H. Yoshida, A. Matsuo, K. Kindo, and Z. Hiroi, Phys. Rev. B 83, 180407(R) (2011); Krishna Kumar, Kai Sun, and Eduardo Fradkin, Phys. Rev. B 92, 094433 (2015); K. Hida, J. Phys. Soc. Jpn. 69(12):40034007, (2000); T. Liu, W. Li, A. Weichselbaum, J. von Delft, and G. Su, Phys. Rev. B 91, 060403(R) (2015); H. J. Changlani and A. M. Läuchli, Phys. Rev. B 91, 100407 (2015); W. Li, A. Weichselbaum, J. von Delft, and H.Hao Tu, Phys. Rev. B 91, 224414 (2015); D. Ixert, T. Tischler, and K. P. Schmidt, Phys. Rev. B 92, 174422 (2015); O. Götze, D. J. J. Farnell, R. F. Bishop, P. H. Y. Li, and J. Richter, Phys. Rev. B 84, 224428 (2011); W. Li, S. Yang, M. Cheng, Z.Xin Liu, and H.Hao Tu, Phys. Rev. B 89, 174411 (2014); S. Nishimoto and M. Nakamura, Phys. Rev. B 92, 140412(R) (2015); K. Okuta, S. Hara, H. Sato, Y. Narumi, and K. Kindo, J. Phys. Soc. Jpn. 80(6):063703 (2011); T. Heng Han, J. S. Helton, S. Chu, D. G. Nocera, J. A. RodriguezRivera, C. Broholm, and Y. S. Lee, Nature, 492: 406410 (2012); D. E. Freedman, R. Chisnell, T. M. McQueen, Y. S. Lee, C. Payen, and D. G. Nocera, Chem. Commun., 48:6466 (2012); S. Hara, H. Sato, Y. Narumi, J. Phys. Soc. Jpn. 81(7):073707 (2012). J. B. Marston and C. Zeng, J. Appl. Phys. 69 (8):59625964, (1991); K. Yang, L. K. Warman, and S. M. Girvin, Phys. Rev. Lett. 70, 2641 (1993); M. B. Hastings, Phys. Rev. B 63, 014413 (2000); F. Wang and A. Vishwanath, Phys. Rev. B 74, 174423 (2006); Y. Ran, M. Hermele, P. A. Lee, and XiaoGang Wen, Phys. Rev. Lett. 98, 117205 (2007); R. R. P. Singh and D. A. Huse, Phys. Rev. B 77, 144415 (2008); G. Evenbly and G. Vidal, Phys. Rev. Lett. 104, 187203 (2010). J.Wei Mei, J.Yao Chen, H. He, X.Gang Wen, arXiv:1606.09639.
 (2) S. Capponi, O. Derzhko, A. Honecker, A. M. Lauchli, and J. Richter, Phys. Rev. B 88, 144416 (2013); S. Nishimoto, N. Shibata, and C. Hotta, Nat. Commun. 4, 2287 (2013).
 (3) K. Kumar, K. Sun, and Eduardo Fradkin, Phys. Rev. B 90, 174409 (2014).
 (4) J. Jordan, R. Orús, and G. Vidal, F. Verstraete, and J. I. Cirac, Phys. Rev. Lett. 101, 250602 (2008); J. Jordan, R. Orús, and G. Vidal, Phys. Rev. B 79, 174515 (2009).
 (5) T. Picot, M. Ziegler, R. Orús and D. Poilblanc, Phys. Rev. B 93, 060407 (2016).
 (6) Z. Y. Xie, J. Chen, J. F. Yu, X. Kong, B. Normand, and T. Xiang, Phys. Rev. X 4, 011025 (2014); H. J. Liao, Z. Y. Xie, J. Chen, Z. Y. Liu, H. D. Xie, R. Z. Huang, B. Normand, and T. Xiang, arXiv:1610.04727.
 (7) Thibaut Picot and Didier Poilblanc, Phys. Rev. B 91, 064415 (2015).
 (8) Ê G. H. Wannier, Phys. Rev. 79, 357 (1950).
 (9) D. C. Cabra, M. D. Grynberg, P. C. W. Holdsworth, A. Honecker, P. Pujol, J. Richter, D. Schmalfuß, and J. Schulenburg, Phys. Rev. B 71, 144420 (2005).
 (10) R. Moessner, S.L. Sondhi, P. Chandra, Phys. Rev. B 64, 144416 (2001).
 (11) S. V. Isakov, S. Wessel, R. G. Melko, K. Sengupta, Yong Baek Kim, Phys. Rev. Lett. 97, 147202 (2006).
 (12) H. C. Jiang, Z. Y. Weng, and T. Xiang, Phys. Rev. Lett. 101, 090603 (2008).
 (13) H. N. Phien, J. A. Bengua, H. D. Tuan, P. Corboz, and R. Orús, Phys. Rev. B 92, 035142 (2015).
 (14) R. Orús and G. Vidal, Phys. Rev. B 80, 094403 (2009);
 (15) F. Verstraete, J. I. Cirac, and V. Murg, Adv. Phys. 57,143 (2008); J. I. Cirac and F. Verstraete, J. Phys. A: Math. Theor. 42, 504004 (2009); R. Augusiak, F. M. Cucchietti, and M. Lewenstein, in Modern Theories of ManyParticle Systems in Condensed Matter Physics, Lect. Not. Phys. 843, 245294 (2012); J. Eisert, Modeling and Simulation 3, 520 (2013); N. Schuch, QIP, Lecture Notes of the 44th IFF Spring School (2013); R. Orús, Eur. Phys. J. B 87, 280 (2014); R. Orús, Ann. Phys.New York 349 117158 (2014).
 (16) J. I. Cirac, D. Poilblanc, N. Schuch, and F. Verstraete, Phys. Rev. B 83, 245134 (2011).
 (17) J. Schulenburg, A. Honecker, J. Schnack, J. Richter, and H.J. Schmidt, Phys. Rev. Lett. 88, 167207 (2002).
 (18) R. Orús, T.Chieh Wei, O. Buerschaper, A. GarcíaSaez, Phys. Rev. Lett. 113, 257202 (2014).
 (19) Daniel Huerga, Sylvain Capponi, Jorge Dukelsky, Gerardo Ortiz, Phys. Rev. B 94, 165124 (2016).