The phase diagram of the antiferromagnetic XXZ model on the triangular lattice
We determine the quantum phase diagram of the antiferromagnetic spin-1/2 XXZ model on the triangular lattice as a function of magnetic field and anisotropic coupling . Using the density matrix renormalization group (DMRG) algorithm in two dimensions we establish the locations of the phase boundaries between a plateau phase with 1/3 Néel order and two distinct coplanar phases. The two coplanar phases are characterized by a simultaneous breaking of both translational and U(1) symmetries, which is reminiscent of supersolidity. A translationally invariant umbrella phase is entered via a first order phase transition at relatively small values of compared to the corresponding case of ferromagnetic hopping and the classical model. The phase transition lines meet at two tricritical points on the tip of the lobe of the plateau state, so that the two coplanar states are completely disconnected. Interestingly, the phase transition between the plateau state and the upper coplanar state changes from second order to first order for large values of .
pacs:75.10.Jm, 67.80.kb, 05.30.Jp
Competing interactions between quantum spins can prevent conventional magnetic order at low temperatures. In the search of interesting and exotic quantum phases frustrated systems are therefore at the center of theoretical and experimental research in different areas of physics balents10 (); anderson (); sorella99 (); ed (); ccm (); tri11 (); tri10 (); tri12 (); tri13 (); watarai01 (); tocchio13 (); balents13 (); white11 (); tocchio14 (); miyashita10 (); seabra11 (); kawamura85 (); classical (); xxz (); yamamoto2 (); starykh14 (); tri (); moessner08 (); imp (); tri_1st (); tri_sc (); tri_af (); metav14 (); yoshikawa04 (); zhang13 (); tri_afn (); jiang08 (). One of the most straight-forward frustrated system is the spin-1/2 antiferromagnet (AF) on the triangular lattice, which was also the first model to be discussed as a potential candidate for spin-liquid behavior without conventional order by Anderson anderson (). It is now known that the isotropic Heisenberg model on the triangular lattice is not a spin liquid and does show order at zero temperature sorella99 (). Nonetheless, the phase diagram as a function of magnetic field is still actively discussed with recent theoretical calculations ed (); ccm () as well as experimental results tri11 (); tri10 (); tri12 (); tri13 () on BaCoSbO, which appears to be very well described by a triangular AF. Interesting phases have also been found for anisotropic triangular lattices tocchio13 (); balents13 (); white11 () and for the triangular extended Hubbard model tocchio14 (). Hard-core bosons with nearest neighbor interaction on a triangular lattice correspond to the xxz model with ferromagnetic exchange in the xy-plane, which has been studied extensively tri (); moessner08 (); imp (); tri_1st (); tri_sc (); tri_af (). In this case a so-called supersolid phase near half-filling has been established for large interactions tri (), which is characterized by two order parameters, namely a superfluid density and a charge density order. Impurity effects show that the two order parameters are competing imp () and the transition to the superfluid state is first order tri_1st (); tri_sc ().
However, surprisingly little attention has been paid to the role of an antiferromagnetic anisotropic exchange interaction away from half-filling classical (); xxz (); yamamoto2 (); starykh14 (), even though the XXZ model on the triangular lattice
is arguable one of the most fundamental examples of frustrated antiferromagnetism. Only very recently the first complete phase diagram as a function of and was published by Yamamoto et al. using the cluster mean-field theory (CMF) xxz (). In this case, three phases with broken sublattice symmetry and simultaneously broken U(1) symmetry were found, which are stable to surprisingly large compared to the corresponding ferromagnetic model. One of those phases – the so-called -coplanar phase – was not expected to exist at all from simple mean field considerations xxz () and therefore deserves special attention.
We now present quantum simulations of this model using the density matrix renormalization group (DMRG) white (); schollwoeck (); xiang () algorithm in two dimensions. The resulting phase diagram as a function of and is summarized in Fig. 1, which first of all confirms several aspects of the previous study in Ref. [xxz, ]: For large values of , we find an umbrella state with spontaneously broken U(1) symmetry, but no broken sublattice symmetry. With increasing and at small magnetic fields a first order transition occurs to a antiferromagnetic coplanar phase where the spins on one sublattice align against the field, while the other two sublattices form a honeycomb structure with spins still partially pointing in the xy-plane, so that all spins lie in a plane. At large fields a ferrimagnetic coplanar phase is found with parallel canted spins on two sublattices and one sublattice pointing in a different direction. A 1/3 Néel phase with fixed magnetization separates the two co-planar phases. The phase transition to the saturated phase occurs exactly at as for the classical triangular antiferromagnet kawamura85 (); xxz (); classical (); seabra11 (); miyashita10 ().
Our results also show several differences to the previous study xxz (): 1.) The so-called -coplanar phase is missing. As shown below this phase exists only for small system sizes or clusters. 2.) Two tri-critical points, which separate the 1/3-Neel phase from the umbrella phase are pushed to much larger values of and become very close in the thermodynamic limit. 3.) The second order phase transition between the 1/3 Néel phase and the ferrimagnetic coplanar phase curiously turns first order for strong interactions at a special bi-critical point, which has since been confirmed yamamoto2 (). Similar bi-critical points where a phase transition changes from 1st to 2nd order were recently under discussion in binary Bose mixtures yamamoto3 ().
We now discuss the detailed numerical DMRG data at selected points in the phase diagram. Frustrated systems are known to be sensitive to boundary induced behavior zhang13 (), so that periodic boundary conditions (PBC) turned out to be necessary in both directions tri_afn (); jiang08 (). Accordingly, the initial truncation error may be as high as which is normal for 2D DMRG with PBC tri_afn (); jiang08 (). In fact, DMRG “sweeping” improves the data significantly (up to 16%), so that the initial truncation error becomes irrelevant as a measure (which is in fact not very sensitive to ). The final energy values after sweeping go to a unique value for large , so that convergence can be ensured. Note, that the DMRG operates in the canonical ensemble, i.e. the data is given as a function of magnetization per site and the corresponding fields can be obtained as the derivative of the ground state energy with respect to , i.e. sq (); anyon (); gap (). The upper tricritical point can be found by the condition . There is no particle-hole symmetry so the kept states we can afford is at most. Technical details about convergence and finite size scaling can be found in the appendix.
The Heisenberg system in a field has previously been considered using exact diagonalization ed (); miyashita86 (); pierre1994 (); honecker1999 (); honecker2004 (), spin waves zheng-06 (); chubukov91 () and coupled cluster methods (CCM) ccm (). It is well known that the uniform magnetization has a plateau at which is characteristic of the 1/3 Néel phase as shown in the inset of Fig. 2.
The structure factors in the z-direction and in the xy-direction at are useful order parameters to measure the diagonal and the off-diagonal order, respectively. If is finite the system has a broken sublattice symmetry (charge order), while a finite indicates a broken U(1) rotational symmetry (superfluidity). As shown in Fig. 2 both order parameters are finite in the ferrimagnetic and antiferromagnetic coplanar phases. At zero magnetization is larger than , but then decreases with and scales to zero with at , which is exactly the point where becomes largest. In the experiments on BaCoSbO an additional cusp in the susceptibility was observed at higher magnetization tri11 (), which could indicate another phase transition. However, our data does not show any other phase for and . Nonetheless, the off-diagonal structure factor does show a broad maximum around , which is due to the fact that the spins on one of the sublattices are able to align along the xy-plane at approximately this magnetization as shown in Fig. 2. Spins that are aligned within the xy-plane have in turn the largest susceptibility in the z-direction, so this could in part explain the observed maximum in Ref. [tri11, ].
We now turn to larger values of , where the magnetization plateau is larger than for as shown in the inset of Fig. 3. The behavior of the order parameters and is qualitatively similar to the isotropic case as a function of magnetization. However, for the phase transition between the 1/3 Néel phase and the ferrimagnetic coplanar there is a subtle, but important difference in the magnetization curve at strong interactions. As shown in Fig. 3, near the upper phase boundary the calculated field decreases with increasing magnetization (which is fixed for each simulation). This behavior indicates an unstable state and in the thermodynamic limit leads to phase separation, which is an obvious indication of a first order phase transition. In a finite system the energy of the phase boundary can prevent phase separation and the unstable state can be found by numerical simulations at a given magnetization, which is the case here and in related systems poilblanc09 (); sq (); yamamoto2 (). The corresponding first order jump in magnetization must then be determined by a Maxwell construction as indicated in Fig. 3. This jump vanishes somewhere between and , so that we predict a bicritical point where the second order phase transition turns first order in the strong coupling limit as shown in Fig. 1. This surprising behavior can in part be explained from the fact that the end of the plateau approaches the saturation field, so that there is only a small field region where the magnetization changes from down to . However, the coplanar spin state has only a limited susceptibility close to saturation, so that a jump in magnetization may be the only way to resolve this contradiction. In other words, starting from the 1/3 Néel state the configuration must make a finite jump to reach the coplanar state if the upper critical field is too large, since the ferrimagnetic coplanar state is already canted significantly towards the field in this case. In any case, the quantum mechanical mechanism for this behavior is an interesting aspect for future studies. The second order phase transition between the antiferromagnetic coplanar phase and the plateau phase is well understood from a strong coupling expansion tri_sc () in terms of holes which start to occupy the honeycomb sublattice at a critical value of , which is consistent with our numerical data.
We should emphasize that order parameters do not have to be used in order to determine the phase transitions from the 1/3-Néel phase to the coplanar states, since the magnetization plateau can be determined directly from the energies . In order to study the phase boundaries to the umbrella phase, on the other hand, order parameters are essential, but this becomes numerically costly for large system sizes. As additional tools, we therefore want to explore here if different measures of entanglement and quantum discord are useful in 2D DMRG, which have been proposed and used for studying quantum phase transitions in recent related systems ee1 (); ee2 (); ee3 (); Wootters (); dillen08 (); QD (). To define suitable quantum information measures it is useful to consider the reduced density matrix of two neighboring spins. The trace over spin gives the reduced density matrix of a single spin . The von-Neumann entropy of a general density matrix can be used to define the entanglement entropy . The concurrence Wootters (); dillen08 ()
is given in terms of the eigenvalues of the matrix , where characterizes the spin-flipped state. The quantum discord QD (); dillen08 () has been proposed as a good indicator for quantum phase transitions
which is calculated in terms of the conditional quantum entropy
where and . The projectors can be defined in terms of a general parametrization
The minimization over the projectors in (3) then corresponds to a minimization over angles and in the wavefunctions.
In Fig. 4 we show the two order parameters, the concurrence, the entanglement entropy, and the quantum discord at two selected points in the phase diagram, which are indicated by black circles in Fig. 1. All measures give the same locations of the phase transition (in this case , and , , respectively). The quantum information measures based on are computationally less demanding than the structure factors since they can be determined from the correlation functions of only two neighboring spins dillen08 (). They are also universal, since no particular order needs to be assumed. In particular, the quantum discord QD (); dillen08 () turns out to be very reliable in detecting the phase transitions and interestingly the corresponding variational angle in Eq. (5) takes on different values on the two sides of the phase transition. It is so far unclear if this jump in a variational parameter is a generic feature, but it may be useful in future studies as well. We find that the phase transition between the ordered states (Néel and coplanar) to the umbrella phase is always first order, except at the isotropic point , where it is known to be second order tri_afn (). At two tricritical points the second order phase transitions between Néel and coplanar phases meet the first order transition. The phase transition to the umbrella phase can be accurately determined for system sizes of up to , so that a systematic finite size scaling becomes feasible. The (interpolated) first order phase transitions to the umbrella phase can be linearly extrapolated in which gives an estimate in the thermodynamic limit (TD) shown in Fig. 1. Extrapolating with a different power also gives a reasonable fit and pushes the phase transition line out even further by up to , which would make an even larger quantitative difference to the coupled cluster study xxz (). The two tri-critical points approach each other with finite size scaling and we cannot rule out that they merge to one single multicritical point in the TD. While finite size scaling works reasonably well for the first order phase transition, the same is not true for the second order phase transition lines, which show a much more irregular behavior with system size, that we cannot explain.
where is the magnetization on each sublattice . The parameter shows three different values in the ferrimagnetic coplanar, the -coplanar and the umbrella phase, respectively as shown in Fig. 5. For small system sizes of or less a -coplanar phase can be identified, but it shrinks fast with increasing system size. With finite size scaling shown in the inset of Fig.5, the -coplanar phase disappears for . Therefore, we predict that there is no such phase in the thermodynamic limit.
In conclusion, we have analyzed the spin-1/2 XXZ model on the triangular lattice using a two dimensional DMRG method with periodic boundary conditions. The phase diagram shows two coplanar phases with different symmetries of the superfluid condensate, which is separated by an ordered plateau 1/3 Néel phase, with fixed magnetization . The transition to the umbrella state is always first order for finite fields and the critical line in Fig. 1 is monotonically increasing, so that a larger field always leads to an extended ordered state. The transition between the coplanar and the 1/3 Néel phase is generically second order but curiously the upper phase transition line turns first order for , which is yet not fully understood.
Acknowledgements.We are thankful for useful discussions with Axel Pelster about his meanfield calculation of the extended Bose-Hubbard model, Shijie Hu about numerical suggestions, Tao Shi and Raoul Dillenschneider about the spin wave calculations and comments from Alexandros Metavitsiadis, Denis Morath and Dominik Straßel. This work was supported by the “Allianz für Hochleistungsrechnen Rheinland-Pfalz” and by the DFG via the SFB/Transregio 49.
Appendix A Appendix: Finite size scaling of the phase diagram
The first order phase transition to the umbrella phase for a given system size can be determined rather accurately as shown in Fig. 4 in the main manuscript. In those simulations the magnetization is fixed and the corresponding magnetic field is determined by the derivative of the ground state energy at the transition point. This yields the phase transition lines for system sizes , , and shown in Fig. 6 below. Since the data is well behaved, it is possible to determine the corresponding continuous curves for all values of by spline interpolation. For each field it is then possible to use a linear fit in reciprocal system size as shown in Fig. 7 (top), which determines the estimate in the thermodynamic limit in Fig. 6. A reliable error estimate is difficult in this case, since the finite size data is only available for three data points and a square root behavior also yielded reasonable fits, which would push out the estimate of the phase transition line to higher values of by up to . Therefore, the phase transition line shown should be taken as a lower estimate.
For each finite size we find two tri-critical points where one coplanar, the 1/3 Néel, and the umbrella phase meet. The location of the tri-critical points change in both and with finite size, but a reasonable estimate of the corresponding values thermodynamic limit can be made as shown in Fig. 7 (lower two plots). Since the two points come quite close with finite size scaling it cannot be ruled out from our data that they merge into one multi-critical point in the thermodynamic limit.
The situation is even more complicated for the second order phase transition lines from the 1/3 Néel phase to the coplanar phase. In this case we could not find any systematic finite size scaling, since the phase transition lines for and are very close but there is a larger change when going to . We cannot explain this behavior, but it is maybe not surprising, that it is more difficult to pinpoint the transition lines for second order phase transitions than for first order transition lines in the thermodynamic limit.
Appendix B Appendix: Convergence of the DMRG data in 2D
In order to study two-dimensional systems with the DMRG algorithm, the sites need to be ordered along a one-dimensional chain with effectively long range interactions. Therefore, neighboring sites may be rather far apart in the DMRG algorithm, so that more effort is required to capture the quantum correlations. This problem is even more severe with periodic boundary conditions in both direction, which were necessary for the xxz model on the triangular lattice. As a result we find that the variational state of the system is rather poorly described after the initial DMRG buildup. Typically, the truncation error is not very small () and is not a reliable measure of the quality of the simulations (in fact it does not depend much on the number of states kept). The finite size algorithm quickly improves this state by “sweeping”. While the correction in energy maybe as large as 20% in the first sweep, the changes become several orders of magnitude smaller after just a few iterations as shown in Fig. 8.
However, the convergence with the number of sweeps is also not a guarantee that the system is approaching the correct ground state, since metastable states are possible. It is therefore essential to vary the number of states kept. It is also possible to change the number of kept states during the sweeping procedure or change the initial buidup geometry. Typically, we find that if the data has a smooth behavior with the number of states kept it also produces sensible and accurate data which fits well into the phase diagram (e.g. relative to neighboring points in parameter space). Hence each data point in the phase diagram has been carefully checked for consistency. In principle it would also be possible to try an extrapolation fit with the number of states, but in practice this only produces tiny corrections to the phase diagram but may in turn produce artifacts.
The following Figs. 9-13 illustrate the typical behavior of the DMRG data for the energy, the structure factors, and quantum information measures at examplary values of the parameters and different system sizes as a function of number of states kept. Note, that the energies alone determine much of the phase diagram since they define the magnetization plateau of the 1/3 Néel phase.
The behavior of the energy and structure factor is shown in Fig. 9-12 for system sizes 6x6 and 6x9. As expected the accuracy improves with number of states kept and is fully sufficient already for ca. 1200 kept states. In particular, the difference between the data for 800 or 1200 kept states would not be noticable in any of the plots. More importantly, there are no big jumps which would be an indicator for metastable states.
Larger system sizes of are more difficult. As shown in Fig. 13 there may be large jumps when going from 600 states to 900 states, which indicates a metastable situtation. However, all paramaters converge to stable values for larger values of the number of states kept. We have checked convergence for up to at selected points.
- (1) L. Balents, Nature (London) 464, 199 (2010).
- (2) P.W. Anderson, Mater. Res. Bull. 8, 153 (1973).
- (3) L. Capriotti, A. E. Trumper, and S. Sorella, Phys. Rev. Lett. 82, 3899 (1999).
- (4) T. Sakai and H. Nakano, Phys. Rev. B 83, 100405(R) (2011).
- (5) D. J. J. Farnell, R. Zinke, J. Schulenburg, and J. Richter, J. Phys. Condens. Matter 21, 406002 (2009).
- (6) T. Susuki, N. Kurita, T. Tanaka, H. Nojiri, A. Matsuo, K. Kindo, and H. Tanaka Phys. Rev. Lett. 110, 267201 (2013);
- (7) Y. Shirata, H. Tanaka, A. Matsuo, and K. Kindo Phys. Rev. Lett. 108, 057205 (2012)
- (8) H.D. Zhou, C. Xu, A.M. Hallas, H.J. Silverstein, C.R. Wiebe, I. Umegaki, J.Q. Yan, T.P. Murphy, J.-H. Park, Y. Qiu, J.R.D. Copley, J. S. Gardner, and Y. Takano, Phys. Rev. Lett. 109, 267206 (2012).
- (9) G. Koutroulakis, T. Zhou, C.D. Batista, Y. Kamiya, J.D. Thompson, S.E. Brown, and H.D. Zhou, preprint arXiv:1308.6331 unpublished (2013).
- (10) S. Watarai, S. Miyashita and H. Shiba, J. Phys. Soc. Jpn. 70, 532 (2001).
- (11) L. F. Tocchio, H. Feldner, F. Becca, R. Valentí, and C. Gros, Phys. Rev. B 87, 035143 (2013).
- (12) R. Chen, H. Ju, H.-C. Jiang, O. A. Starykh, and L. Balents, Phys. Rev. B 87, 165123 (201 3).
- (13) A. Weichselbaum and S.R. White, Phys. Rev. B 84, 245130 (2011); K. Harada, Phys. Rev. B 86, 184421 (2012).
- (14) L.F. Tocchio, C. Gros, X-F. Zhang, S. Eggert, Phys. Rev. Lett. 113, 246405 (2014).
- (15) S. Wessel and M. Troyer, Phys. Rev. Lett. 95, 127205 (2005); D. Heidarian and K. Damle, Phys. Rev. Lett. 95, 127206(2005); R. G. Melko, A. Paramekanti, A. A. Burkov, A.Vishwanath, D. N. Sheng, and L. Balents, Phys. Rev.Lett. 95, 127207 (2005); M. Boninsegni and N. Prokof’ev , Phys. Rev. Lett. 95, 237204 (2005).
- (16) A. Sen, P. Dutt, K. Damle, and R. Moessner, Phys. Rev. Lett. 100, 147204 (2008).
- (17) X. F. Zhang, Y.C. Wen, and S. Eggert, Phys. Rev. B 82, 220501(R) (2010).
- (18) F. Wang, F. Pollmann, and A. Vishwanath, Phys. Rev. Lett. 102, 017203 (2009). D. Heidarian and A. Paramekanti, Phys. Rev. Lett. 104, 015301 (2010).
- (19) D. Yamamoto, I. Danshita, and C.A.R. Sá de Melo Phys. Rev. A 85, 021601(R) (2012); L. Bonnes and S. Wessel, Phys. Rev. B 84, 054510 (2011); D. Yamamoto, T. Ozaki, C.A.R. Sá de Melo, I. Danshita Phys. Rev. A 88, 033624 (2013).
- (20) X.-F. Zhang, R. Dillenschneider, Y. Yu, and S. Eggert Phys. Rev. B 84, 174515 (2011);
- (21) S. Miyashita, Proc. Jpn. Acad., Ser. B 86, 643 (2010).
- (22) L. Seabra, T. Momoi, P. Sindzingre, and N. Shannon, Phys. Rev. B 84, 214418 (2011).
- (23) H. Kawamura and S. Miyashita, J. Phys. Soc. Jpn. 54, 4530 (1985).
- (24) S. Miyashita, J. Phys. Soc. Jpn. 55, 3605 (1986).
- (25) D. Yamamoto, G. Marmorini, I. Danshita, Phys. Rev. Lett. 112 127203 (2014).
- (26) D. Yamamoto, G. Marmorini, I. Danshita, Phys. Rev. Lett. 112, 259901(E) (2014).
- (27) O.A. Starykh, W. Jin, A.V. Chubukov, Phys. Rev. Lett. 113, 087204 (2014).
- (28) A. Metavitsiadis, R. Dillenschneider, and S. Eggert, Phys. Rev. B 89, 155406 (2014).
- (29) S. Yoshikawa, K. Okunishi, M. Senda and S. Miyashita, J. Phys. Soc. Jpn. 73, 1798 (2004).
- (30) X-F. Zhang and S. Eggert Phys. Rev. Lett. 111, 147201 (2013).
- (31) H.C. Jiang, M.Q. Weng, Z.Y. Weng, D.N. Sheng, and L. Balents, Phys. Rev. B 79, 020409(R) (2009).
- (32) H.C. Jiang, W.Y. Weng, and D.N. Sheng, Phys. Rev. Lett. 101, 117203 (2008).
- (33) S. R. White. Phys. Rev. Lett. 69, 2863, (1992); E.M. Stoudenmire and S.R. White Annu. Rev. Condens. Matter Phys. 3, 111 (2012)
- (34) U. Schollwöck. Rev. Mod. Phys., 77, 259 (2005); Annals of Physics, 326, 96 (2011).
- (35) T. Xiang, J. Lou, and Z. Su. Phys. Rev. B, 64, 104414 (2001).
- (36) Y. Kato, D. Yamamoto, and I. Danshita, Phys. Rev. Lett. 112, 055301 (2014).
- (37) For large magnetization exact diagonalization gives more accurate results than DMRG since the Hilbert space is relatively small.
- (38) G. G. Batrouni and R. T. Scalettar, Phys. Rev. Lett. 84, 1599 (2000)
- (39) T. Keilmann, S. Lanzmich, I. McCulloch and M. Roncaglia, Nature Communication, 2, 361 (2011)
- (40) For system size the even magnetization has a better accuracy, so we choose .
- (41) H. Nishimori, and S. Miyashita, J. Phys. Soc. Jpn. 55, 4448 (1986).
- (42) B. Bernu, P. Lecheminant, C. Lhuillier, and L. Pierre, Phys. Rev. B 50, 10048 (1994).
- (43) A. Honecker, J. Phys.: Condens. Matter 11, 4967 (1999).
- (44) A. Honecker, J. Schulenburg, and J. Richter, J. Phys.: Condens. Matter 16, S749 (2004).
- (45) W. Zheng, J. O. Fjærestad, R. R. P. Singh, R. H. McKenzie, and R. Coldea, Phys. Rev. B 74, 224420 (2006).
- (46) A.V. Chubukov and D.I. Golosov, J. Phys.: Condens. Matter 3, 69 (1991).
- (47) M. Raczkowski and D. Poilblanc, Phys. Rev. Lett. 103, 027001 (2009).
- (48) A. Osterloh, L. Amico1, G. Falci and R. Fazio, Nature 416, 608 (2002)
- (49) Ö. Legeza and J. Sólyom, Phys. Rev. Lett. 96, 116401 (2006)
- (50) L. Amico, R. Fazio, A. Osterloh, and V. Vedral, Rev. Mod. Phys. 80, 517 (2008)
- (51) W.K. Wootters, Phys. Rev. Lett. 80, 2245 (1998).
- (52) R. Dillenschneider, Phys. Rev. B78, 224413 (2008).
- (53) H. Ollivier and W.H. Zurek, Phys. Rev. Lett. 88, 017901 (2001).