Ground-state phase diagram of the quantum J1-J2 model on the honeycomb lattice
aa r X i v : . [ c ond - m a t . s t r- e l ] F e b Ground-state phase diagram of the quantum J − J model on the honeycomb lattice Fabio Mezzacapo and Massimo Boninsegni Max-Planck-Institut f¨ur Quantenoptik, Hans-Kopfermann-Str.1, D-85748, Garching, Germany Department of Physics, University of Alberta, Edmonton, Alberta, Canada T6G 2E1 (Dated: October 2, 2018)We study the ground-state phase diagram of the quantum J − J model on the honeycomblattice by means of an entangled-plaquette variational ansatz. Values of energy and relevant orderparameters are computed in the range 0 ≤ J /J ≤
1. The system displays classical order for J /J < ∼ . J /J > ∼ . PACS numbers: 02.70.-c, 71.10.Fd, 75.10.Jm.
Frustrated quantum antiferromagnets are a subject ofcurrent intense research, whose study is rendered evenmore timely by recent, intriguing proposals of possibleexperimental realizations with ultracold atoms. Frustra-tion can arise either from the geometry of the system, orfrom competing interactions, and can lead to the sta-bilization of novel, exotic (magnetic and non-magnetic)phases of matter.A prototypical spin model, featuring interaction-induced frustration, is the spin-1/2 antiferromagnetic(AF) Heisenberg Hamiltonian in presence of next-nearest-neighbor (NNN) coupling (usually referred to asthe J − J model): J X h i,j i S i · S j + J X hh i,j ii S i · S j (1)where S i is a spin-1/2 operator associated to the i th lat-tice site and the first (second) summation runs over NN(NNN) sites. Periodic boundary conditions (PBC) areassumed. Model (1) has been extensively studied on the squarelattice, and has more recently elicited great intereston the honeycomb one, . For J = 0 (i.e., with noNNN interaction), the ground state (GS) of (1) featuresAF long-range (N´eel) order on both these two bipartitelattices. When J > J , but as J grows different phases, including disordered ones, becomeenergetically competitive. For example, on the squarelattice the GS displays magnetic long-range order (al-beit of different types) for small and large value of J ,while the nature of the GS in the intermediate region(i.e., J ∼ .
5) is still under debate. On the honeycomblattice, due to its lower coordination number with respectto the square one, the disruptive effect of quantum fluc-tuations on magnetic order is enhanced, as confirmed bythe value of the sublattice magnetization, which in theunfrustrated case is approximatively 10% smaller on thehoneycomb than on the square lattice .Studies of the magnetic phase diagram of the J − J model on the honeycomb lattice, based on different com-putational approaches, have yielded conflicting physi-cal scenarios. For example, it is unclear if the GS of FIG. 1: (color on line). Schematic representation of the hon-eycomb lattice, with the two different definitions of sublattices( A and B ) corresponding to N´eel (left) and collinear (right)orders. the system is disordered (“spin liquid”) for any valueof J , and even the nature of ordered phases remainscontroversial. Exact diagonalization (ED) of (1) onsmall lattices (up to 42 sites), yields evidence of AFN´eel order for J < ∼ .
2. For 0 . < ∼ J < ∼ . J , a collinearphase (not easy to characterize on small lattices) has beensuggested (see Fig. 1). The physical picture emergingfrom ED, whose predictive power is clearly affected byfinite-size limitations, must be assessed by calculationsperformed on systems of much larger size, in order tocarry out a reliable extrapolation of the relevant physicalproperties to the thermodynamic limit. Series expansionsyield, for instance, a similar physical scenario, the maindifference being a GS with spiral magnetic order foundin the intermediate region. Quantum Monte Carlo (QMC) approaches yield nu-merically exact estimates for very large system sizes at J = 0, but are not applicable in presence of frustrationdue to the well-known “sign” problem. Conversely, vari-ational Monte Carlo (VMC) methods are by definitionfree of such problem, and suffer from no significant sizelimitation. Of course, they are also approximate, theiraccuracy depending on the choice of the variational wavefunction (WF). A recent variational investigation of the J − J Hamiltonian on the honeycomb lattice, predictsthe disappearance of N´eel order at J ≃ .
08, a dimerizedrotational symmetry-breaking phase for J > ∼ . It should benoted that, while in studies based on ED the presenceof order of different kinds is assessed by a computationof relevant order parameters, the conclusions of Ref. 4are mainly based on the comparison of energy estimatesyielded by variational WF’s corresponding to differentlyordered, or disordered GS’s.Given the qualitative disagreement between ED andthe only existing variational study, two basic aspects thatstill need be clarified are a ) the type of magnetic order (ifany) of phase(s) that occur as N´eel order is suppressed onincreasing J , and b ) the nature of the GS for J > ∼ . J at which the N´eel or-der vanishes, which according to Ref. 4 is in quantitativeagreement with that found for the half-filled Hubbardmodel. In this paper, we report results of a variational studyof the GS phase diagram of (1) on the honeycomb lat-tice, using the variational family of entangled-plaquettestates (EPS). This is a general ansatz, which has beenshown to yield accurate estimates of GS observables forquantum spin lattice Hamiltonians, including frustratedones.
The goal of this work is to determine the phasediagram of the model of our interest by direct estima-tion of the values of various order parameters, extrapo-lated to the thermodynamic limit. Specifically, besidesGS energies we also compute sublattice magnetizationsand dimer-dimer correlation functions. We find AF N´eelorder up to J < ∼ .
2, whereas for J > ∼ . . < ∼ J < ∼ . N = 324 sites. For each lattice sizeand value of J , we independently and systematically op-timize the EPS WF, by progressively increasing the sizeof our plaquettes (i.e., the number l of lattice sites thata plaquette comprises). Both the WF optimization andthe evaluation of various physical quantities are achievedvia the variational Monte Carlo method. We extrapo-late to the thermodynamic limit the estimates obtainedon finite lattices of different sizes, for a given plaquettesize, and then, asses the dependence of our results on theplaquette size.Table I shows the GS energy per site e ≡ E/N for a lat-tice comprising N = 100 sites, for different values of J ;estimates are shown corresponding to different plaquettesizes. Also shown for comparison are the variational es-timates reported in Ref. 4 by Clark et al. For all valuesof J , our energy estimates are lower than those of Ref. TABLE I: GS energy per site as a function of J for a system ofsize N = 100 with PBC. Data for different plaquette-size l areshown. Variational estimates from Ref. 4 are also reportedfor comparison. Statistical errors (in parentheses) are on thelast digit. J l = 9 l = 16 l = 18 Ref. 40.0 -0.54155(7) -0.54418(6) -0.537901(9)0.1 -0.49049(5) -0.49349(6) -0.488209(3)0.15 -0.46744(5) -0.47087(4) -0.469448(2)0.2 -0.44599(4) -0.45004(6) -0.45145(4) -0.450687(4)0.3 -0.40936(6) -0.41580(7) -0.41751(4) -0.41688(1)0.4 -0.40650(5) -0.41185(5) -0.41125(1)0.5 -0.42432(5) -0.42940(6) -0.42808(2)0.6 -0.45248(5) -0.45712(4)0.7 -0.48519(4) -0.48989(5) -0.47377(3)0.8 -0.52112(6) -0.52610(5)0.9 -0.55928(4) -0.56489(5)1.0 -0.59917(5) -0.60534(5)
4. For J ≤ .
1, as well as for large J , we obtain lowervariational estimates than those of Ref. 4 already withplaquettes of size l =9. Plaquettes of larger size, typically l = 16 (and up to l = 18 for values of 0 . < ∼ J < ∼ . J and different plaquette size l are shown inFig. 2 (left part), as a function of the system size. Ex-trapolation to the thermodynamic limit of numerical es-timates obtained with the same l has been performedby assuming the functional scaling form e ( N, l, J ) = e ( ∞ , l, J ) + α ( l, J ) N − / . Extrapolated values, there-fore, depend on the size l of the plaquette utilized. Sinceby construction expectation values (of any observable)must approach the exact GS values in the limit of largeplaquette size, the question is how to obtain a reliableestimate for such a limit, based on those obtained withrelatively small plaquettes. While in principle calcula-tions with ever increasing plaquette size should be car-ried out, in order to observe numerical convergence ofresults, currently available computer resources limit thelargest plaquette size for which this procedure can be im-plemented in practice. As it turns out, however, resultsfor both the energy and for the relevant magnetic orderparameters obtained with plaquettes of size up to l =16,lend themselves to simple numerical extrapolations. An example of this procedure is illustrated in the rightpart of Fig. 2, showing energy estimates extrapolated tothe thermodynamic limit i.e., e ( ∞ , l, J ), as a functionof the plaquette size. We assume that e ( ∞ , l, J ) canbe expanded in powers of l − / and fit the data usingthe smallest number of powers. As shown in the fig-ure, the values of e ( ∞ , l, J ) fall on a straight line when −0.55−0.54 0.000 0.002 0.004 e ( N , l , J = ) N −3/2 −0.47−0.46 e ( N , l , J = . ) −0.45−0.44−0.43 e ( N , l , J = . ) e ( ¥ , l , J = ) l −3/2 e ( ¥ , l , J = . ) e ( ¥ , l , J = . ) FIG. 2: (color on line). Left: GS energy per site of the J − J model on the honeycomb lattice with PBC, as a function ofthe system size. Estimates are obtained with N plaquettesof l = 4 (boxes), l = 9 (triangles) and l = 16 (circles) sites.Each dashed line is a fit to numerical data with the same l value. Right: GS energy per site, extrapolated to the ther-modynamic limit, as a function of the plaquette-size. Eachdashed line is a fitting function linear in l − / (see text). TheQMC estimate in the J = 0 case is also shown for comparison(star). plotted versus l − / , i.e., an excellent fit to the data isobtained with a single power (dashed lines); in any case,the extrapolated value does not change appreciably onincluding additional terms. For J = 0, this procedureleads to an estimate in agreement with that obtained byQMC (essentially exact in this case).The cogent observable, in order to assess the pres-ence of magnetic order in the J − J Hamiltonian, is thesquared sublattice magnetization (SSM), defined as m ( N ) = (cid:28) N ( X i ∈ A S i − X j ∈ B S j ) (cid:29) (2)where the two summations are taken over lattice sitesbelonging to different sublattices, as shown in Fig. 1.The dependence of the values of the above observableon J , allows one to identify any range of parameters inwhich possibly magnetically disordered phases may exist.It is worth restating that in this work we consider onlytwo types of magnetic long-range order, schematicallyillustrated in Fig. 1, namely N´eel and collinear.Figure 3 shows estimates of the N´eel SSM extrapolatedto the thermodynamic limit (left part), based on the func-tional form m ( N, l, J ) = m ( ∞ , l, J ) + β ( l, J ) N − / , m ( N , l , J = ) N −1/2 m ( N , l , J = . ) m ( N , l , J = . ) m ( ¥ , l , J = ) l −1/2 m ( ¥ , l , J = . ) m ( ¥ , l , J = . ) FIG. 3: (color on line). Left: N´eel SSM of the J − J modelon the honeycomb lattice with PBC, as a function of the sys-tem size. Estimates are obtained with N plaquettes of l = 4(boxes), l = 9 (triangles) and l = 16 (circles) sites. Eachdashed line is a fit to numerical data with the same l value(see text). Right: N´eel SSM, extrapolated to the thermo-dynamic limit, as a function of the plaquette-size. Lines arefunctions bult to infer the l dependence of our data (see text).The QMC estimate in the J = 0 case is also shown for com-parison (star). for different plaquette sizes. Just like for the energy, theright part of the figure shows the extrapolated values m ( ∞ , l, J ), this time plotted versus l − / . The fit tothe data is performed using a quadratic functional de-pendence (dashed lines).For J = 0, our estimates in the large l limit isagain in quantitative agreement with QMC ones . At J = 0 .
15 the N´eel SSM is finite; we make such a state-ment based on the observation that assuming a linear de-pendence of m ( ∞ , l, J ) on l − / (the fit is acceptable),we obtain a finite extrapolated value [dotted line in Fig.3 (right part)]. On including a term proportional to l − the visible upward bend of the data is better described.In turn, the extrapolated value of the order parameterincreases.For J = 0 . J = 0 . . < ∼ J < ∼ . m ( N , l , J = . ) N −1/2 m ( ¥ , l , J = . ) l −1/2 FIG. 4: (color on line). Same of Fig. 3 for J = 0 .
5. Here A and B sublattices are chosen as in the right part of Fig.1. The dependence of the collinear order parameter on theplaquette size is described by using a function which includesbesides the zero-th order, orders higher than the quadratic (in l − / ) one (dashed line) and only the linear one, excluding thepoint at l = 4 (dotted line). Both descriptions yield, when l is saturated, a finite collinear SSM. s vb c ( N ) N −1 J =0.3 J =0.2 FIG. 5: (color on line). Finite size scaling of the plaquettevalence-bond crystal order parameter for two values of J atwhich the system is found in a magnetically disordered GS.Dashed lines are fits to the numerical data with the same J . This order parameter clearly vanishes, for all the l valuesconsidered in this work, in the thermodynamic limit. Datashown refer to the l = 16 case. larger J , the GS of the system orders collinearly.Figure 4 shows the collinear SSM as a function of N (left part), and its values extrapolated to the thermody-namic limit versus l (right part), for J = 0 .
5. The rightpanel shows that the collinear order parameter is finite.Collinear order has been found up to J = 1 . J value considered in this study). Our prediction of N´eel order persisting up to J ∼ . That thecollinearly ordered phase appears at J ∼ . As shown by our energyestimates (Tab. I), the collinear phase is favored with re-spect to the dimerized rotational-symmetry breaking one(described via a resonating valence-bond WF ) proposedfor J > ∼ . J − J model on the honeycomblattice has been investigated by means of variational cal-culations based on the EPS WF. The phase diagram ofthe model displays three distinct regions: magnetic or-der of N´eel and collinear type is found for J < ∼ . J > ∼ . J at which N´eelorder vanishes, according to both ED and our calcula-tions, is considerably larger than that ( J ≃ .
08) foundin Ref. 4. In Ref. 4, a striking similarity was suggestedbetween the phase boundary of the J − J and that ofthe half-filled Hubbard model at the Neel-to-disordertransition; the validity of such a suggestion may have tobe reconsidered, in the light of the results presented herewhich indicate that the simple reduction of the physicsof the Hubbard model to that of an effective spin-1/2system, such as the J − J , might not be achievable.We acknowledge discussions with J. I. Cirac and H.-H.Tu, and thank B. K. Clark for providing us with his en-ergy estimates. This work has been supported by the EUproject QUEVADIS, and the Canadian NSERC throughthe grant G121210893. J. Struck, C. ¨Olschl¨ager, R. Le Targat, P. Soltan-Panahi,A. Eckardt, M. Lewenstein, P. Windpassinger and K. Sen-gstock, Science , 996 (2011). Henceforth, we set the value of the NN AF coupling con-stant J to unity, and take it as our energy unit. V. Murg, F. Verstraete and J. I. Cirac, Phys. Rev. B ,195119 (2009) and references therein. B. K. Clark, D. A. Abanin and S. L. Sondhi, Phys. Rev.Lett. , 087204 (2011). H. Mosadeq, F. Shabazi and S. A. Jafary, J. Phys. Con-dens. Matter , 226006 (2011). A. F. Albuquerque D. Schwandt, B. Het´enyi, S. Capponi,M. Mambrini and A. M. L¨auchli, Phys. Rev. B , 024406(2011). A. Mulder, R. Ganesh, L. Capriotti and A. Paramekanti,Phys. Rev. B , 214419 (2010). J. Reuther, D. A. Abanin and R. Thomale, Phys. Rev. B , 014417 (2011). D. C. Cabra, C. A. Lamas and H. D. Rosales, Phys. Rev.B , 094506 (2011). J. Oitmaa and R. R. P. Singh, Phys. Rev. B , 094424(2011). Y.-M. Lu and Y. Ran, Phys. Rev. B , 024420 (2011). H. Y. Yang and K. P. Schmidt, Eur. Phys. Lett. , 17004(2011). A. W. Sandvik, Phys. Rev. B , 11678 (1997); E. V.Castro, N. M. R. Peres, K. S. D. Beach, and A. W. SandvikPhys. Rev. B , 054422 (2006). Z. Y. Meng, T. C. Lang, S. Wessel, F. F. Assaad and A.Muramatsu, Nature , 847 (2010). F. Mezzacapo, N. Schuch, M. Boninsegni and J. I. CiracNew J. Phys. , 083026 (2009). F. Mezzacapo and J. I. Cirac, New J. Phys. , 103039(2010); F. Mezzacapo, Phys. Rev. B , 115111 (2011); H.J. Changlani, J. M. Kinder, C. J. Umrigar and G. Kin-LicChan, Phys. Rev. B , 245116 (2009); S. Al-Assam, S. R. Clark, C. J. Foot and D. Jaksch, Phys. Rev. B , 205108(2011). The large l limit of our estimates is achieved by using,starting from plaquettes of a given size (in our case foursites), data obtained with plaquettes whose size is system-atically and consistently increased (i.e., by considering pla-quettes with 2 l / + 1 more sites each time that the pla-quette size l is changed). J. D. Reger, J. A. Riera and A. P. Young, J. Phys. Condens.Matter , 1855 (1989). For J → ∞ the Hamiltonian studied here reduces to thatof two decoupled antiferromagnets on the triangular latticewhich features 2 π/ C. N. Varney, K. Sun, V. Galitski and M. Rigol, Phys. Rev.Lett.107