Critical Correlations for Short-Range Valence-Bond Wave Functions on the Square Lattice
aa r X i v : . [ c ond - m a t . s t r- e l ] N ov Critical Correlations for Short-Range Valence-Bond Wave Functions on the Square Lattice
A. Fabricio Albuquerque and Fabien Alet Laboratoire de Physique Th´eorique, Universit´e de Toulouse and CNRS, UPS (IRSAMC), F-31062 Toulouse, France (Dated: October 17, 2018)We investigate the arguably simplest SU (2) -invariant wave functions capable of accounting for spin-liquidbehavior, expressed in terms of nearest-neighbor valence-bond states on the square lattice and characterized bydifferent topological invariants. While such wave-functions are known to exhibit short-range spin correlations,we perform Monte Carlo simulations and show that four-point correlations decay algebraically with an exponent . . This is reminiscent of the classical dimer problem, albeit with a slower decay. Furthermore, thesecorrelators are found to be spatially modulated according to a wave-vector related to the topological invariants.We conclude that a recently proposed spin Hamiltonian that stabilizes the here considered wave-function(s) asits (degenerate) ground-state(s) should exhibit gapped spin and gapless non-magnetic excitations. PACS numbers: 75.10.Kt,75.10.Jm,75.40.Mg
Introduction —
The quest for quantum spin-liquid (QSL)states of matter is a longstanding research topic that canbe traced back to Anderson’s proposal. Building on earlierwork, he conjectured that strong quantum fluctuations, en-hanced by frustration and/or low coordination, would weaken SU (2) -broken order and occasionally drive an antiferromag-net towards a “disordered” state with exponentially decayingspin correlations, describable in terms of short-ranged spin-singlet, or valence-bond (VB), degrees of freedom. Interestin Anderson’s insight was further triggered in connection withthe cuprates since, soon after their discovery, spin-singlets ina QSL were interpreted as “pre-formed Cooper pairs” thatwould superconduct upon doping. Major advances (reviewed in Refs. 1, 5, and 6) have takenplace since the original proposal, including a classification ofpossible QSL states and explicit realizations in lattice mod-els in dimension d > . Also, theoretical ideas put for-ward in the early days of high- T c have been considerablydeveloped and resulted in a full-fledged formalism as wellas efficient numerical approaches for handling VB states.On the experimental side, a number of compounds have beenshown not to display magnetic order down to the lowest ac-cessible temperatures, much below the energy scale set byexchange interactions, and are thus candidates for the realiza-tion of QSLs. However, in spite of such advances, a completecharacterization of QSL states is still missing, precluding un-ambiguous identification of experimental realizations, sinceabsence of magnetic order does not exclude the occurrence of,for instance, more conventional valence-bond crystals (VBC)that break lattice symmetries (see e.g. Ref. 21).Within this context, we investigate a family of nearest-neighbor VB (NN-VB) states on the square lattice by per-forming Monte Carlo (MC) simulations based on a recentlyintroduced algorithm. We revisit the pioneering work bySutherland, where a closely related NN-VB state was in-troduced, and provide a thorough characterization of the ar-guably simplest SU (2) -invariant wave functions capable ofaccounting for QSL behavior. Although it has long beenknown that NN-VB states on the square lattice are non-magnetic, the possibility of other types of order, such asVBC, has not yet been excluded. Despite their simplicity,and consequent theoretical appeal, the NN-VB states possess highly non-trivial properties, that we explore in what follows. Wave functions —
NN-VB states are obtained by contract-ing spins attached to NN sites i and j of a lattice into a singletstate, [ i, j ] = √ ( | ↑ i ↓ j i − | ↓ i ↑ j i ) . Since each spin onlypairs with one of its neighbors at a time, there is a one-to-onecorrespondence between NN-VB configurations and those ofhard-core classical dimers on the same lattice. A crucial property of VB states is their non-orthogonality.The overlap between two VB configurations is given by h ψ | ψ i = ± N L − N , where N = L is the number ofsites and N L the number of loops in the transition graph ob-tained by superposing the dimer configurations associated to | ψ i and | ψ i [Fig. 1(a)]. For the square and other bipartitelattices, that can be split into two sublattices A and B , overlapsbetween arbitrary VB states can be ensured to be always posi-tive , so that stochastic methods apply (see below), by choosing i ∈ A and j ∈ B in the anti-symmetric singlet [ i, j ] .An additional important point concerns the fact that dimerconfigurations can be split into topological sectors . On atorus, the transition graph for two dimer coverings belongingto different topological sectors displays non-local loops thatwind around the system [Fig. 1(a)], so that one dimer config-uration can not be continuously deformed onto the other via local dimer rearrangements. For bipartite lattices the num-ber of topological sectors is extensive and each sector can belabelled by topological invariants termed winding numbers , w = ( w x , w y ) : w x ( w y ) is defined as the difference betweenthe number of B ← A and
A → B dimers along a referenceline in the y ( x ) direction [Fig. 1(a)]. VB configurations char-acterized by different w are orthogonal in the thermodynamiclimit: the transition graph between two such configurations, | c w i and | c w i , contains at least one winding loop that com-prises a minimum of L dimers, so that N L ≤ ( N − L ) / , im-plying that h c w | c w i ≤ − L and vanishes when L → ∞ . Having introduced the ingredients, we are able to writedown the NN-VB wave functions we wish to investigate: | ψ w i = X c w | c w i . (1)In contrast to the wave function analyzed by Sutherland, who did not take the existence of topological sectors into ac- L -1 -0.35-0.3-0.25-0.2 ! S . S e " (0, 1): ! e x - (1, 0): ! e y (0, 1): -e y - (1, 0): +e x (0, 1): +e y - (1, 0): -e x (1, 1): -e x, y (1, 1): +e x, y (0, 0): ! e x, y w x w y (a) (b) FIG. 1. (Color online) (a) Transition graph between two NN-VBconfigurations on the square lattice (sublattices A / B are indicated byfilled/open circles). Reference lines for the winding numbers w =( w x , w y ) are indicated by dashed lines: configurations with w =(1 , (red lines) and w = (0 , (blue lines) are shown. Trivialloops with coinciding VBs are depicted as thick-black lines and aloop winding in both directions is evident. (b) NN spin correlationsversus L − (from MC simulations) for winding sectors w = (0 , , (0 , , (1 , and (1 , . For w = (0 , , correlations along ± e x and ± e y are discriminated. Lines are linear-quadratic fits. count and considered an equal amplitude superposition of all NN-VB states, each | ψ w i is an equal amplitude superposi-tion of VB configurations | c w i with fixed winding numbers w = ( w x , w y ) . Our motivation for doing so is our previousremark that h ψ w | ψ w i = 0 for w = w in the thermo-dynamic limit. Interestingly, this implies that each | ψ w i is a(degenerate) ground-state on a torus of the local spin Hamil-tonian recently proposed by Cano and Fendley. Althoughthe number of winding sectors is extensive on a torus, we willshow that it is possible to infer the properties of arbitrary states | ψ w i by focusing on a few sectors with low w . Algorithms —
The expectation value of an observable O inEq. (1) can be measured by evaluating hOi w = 1 Z X c w , ,c w , h c w , |O| c w , ih c w , | c w , i h c w , | c w , i , (2)where Z = h ψ w | ψ w i is the normalization. As first pointedout in Ref. 17, hOi w can be efficiently computed in a stochas-tic manner: the estimator h c w , |O| c w , i / h c w , | c w , i , thatfor most observables of interest is readily evaluated by an-alyzing the loop structure in the transition graphs, is sam-pled by generating pairs of NN-VB configurations | c w , i and | c w , i according to the statistical weight given by their overlap h c w , | c w , i = 2 N L ( w ;1 , − N [ N L ( w ; 1 , denotes the num-ber of loops in the transition graph h c w , | c w , i ].Major advances in sampling techniques have been achievedsince the work by Liang et al. and particularly well suited toour purposes is a recently introduced algorithm. Basically(we refer to Ref. 22 for details), one combines non-local up-dates for the underlying dimer configurations, so to ensuresmall auto-correlations times, with spin updates that allowfor an efficient sampling of the overlap weight, with unitaryacceptance rate. By relying on this algorithm, we simulatesystems with periodic boundary conditions (PBC) of sizesof up to L = 128 . Topological symmetry is easily imple-mented in the simulations by starting from a configuration in l og [ (- ) x + y 〈 S . S r 〉 ] (a) -12 -8 -4 0 4 8 12 x -12-8-4 0 4 8 12 y -8-6-4-2 0 -9 -8 -7 -6 -5 -4 -3 -2 -1 (- ) x + y 〈 S . S r 〉 (b) ξ = 1.35(1) FIG. 2. (Color online) (a) Spin correlation ( − r h S · S r i betweenthe spin at the origin and spins located at r = x e x + y e y . (b) ( − r h S · S r i versus distance r . An exponential fit, ( − r h S · S r i ∼ exp( − r/ξ ) , yields the correlation length ξ = 1 . . Datafor L = 128 and w = (0 , . a given winding sector and discarding measurements for allMC moves that change w : for large L winding updates areexponentially rare and this only causes small efficiency losses. Short-range spin order —
We start by analyzing the spintexture in wavefunction Eq. (1). NN spin correlations, h S r · S r ± e α i ( e α is the unit vector in the α = x, y direction) areplotted as a function of inverse system size L − in Fig. 1(b)for w = (0 , , (1 , , (0 , and (1 , . Deviations amongresults for different w are observed for small systems and,additionally, for w = (0 , vertical/horizontal and A → B / B → A correlations differ. However, all data converge tothe same value in the thermodynamic limit, according to alinear-quadratic best-fit analysis. Although we have no rigor-ous justification for such scaling, we obtain h S r · S r + e α i = − . when L → ∞ for all sectors. Spin correla-tions h S · S r i as a function of distance are plotted in Fig. 2(a-b) for L = 128 and w = (0 , (virtually identical resultsare obtained for other L and w ) and display a perfect stag-gered pattern consistent with (isotropic) short-range N´eel or-der. h S · S r i decays very fast with | r | and an exponential fit[Fig. 2(b)] yields ξ = 1 . for the spin correlation length. Critical correlations —
We proceed to the characterizationof “dimer order” by analyzing the four-point connected cor-relators C ijkl = h ( S i · S j )( S k · S l ) i − h S i · S j ih S k · S l i ,where both i, j and k, l are NN sites on the square lattice. InFig. 3 we show MC data for the spatial dependence of rC ijkl ( r is the distance between dimers) for L = 16 and sectors w = (0 , , (1 , , (0 , and (1 , . We first notice that both C ijkl k (correlations for parallel dimers i , j and k , l ) and C ijkl ⊥ (perpendicular dimers i , j and k , l ) are spatially modulatedfor w = (0 , . Inspection of the results in Fig. 3, and similarones for higher w (not shown), allows us to deduce that mod-ulation for C ijkl k [ C ijkl ⊥ ] is entirely accounted for by a phasefactor cos( Q · r ) [ sin( Q · r ) ], with a wave-vector given in termsof the winding numbers, Q = πL ( w y , w x ) . This inference isconfirmed by our quantitative analysis below. Furthermore,we notice that no clear spatial dependence is noticeable for rC ijkl k in Fig. 3(a), suggesting that four-point correlations inEq. (1) decay algebraically with r with an exponent close to (b) - w = (1, 0)(a) - w = (0, 0) (d) - w = (1, 1)(c) - w = (0, 1) FIG. 3. (Color online) MC results for rC ijkl , for L = 16 (PBC) andindicated w . In all panels, the reference bond is indicated by a thick-black line and the thickness of the remaining ones is proportionalto rC ijkl : blue (pale-red) lines indicate positive (negative) values.Lines along which C ijkl k is strongest are indicated by black circles. r -3 -2 -1 C || ij k l ( r ) L = 16L = 32L = 48L = 64L = 80L = 96L = 112L = 128 r -3 -2 -1 | C || ij k l ( r ) | (1, 1)(0, 1)(1, 0)(0, 0) (a) - = (0, 0) w (b) - L = 128 r -3 -2 -1 | C || ij k l ( r ) / c o s ( Q . r) | , | C ⊥ ij k l ( r ) / s i n ( Q . r) | |C ||ijkl (r) / cos (Q. r) ||C ⊥ ijkl (r) / sin (Q . r) | r -6 -5 -4 -3 -2 -1 | C ⊥ ij k l ( r ) | NegativePositive
L = 128(c) - = (0, 1) w (d) - = (0, 0) w L = 128
FIG. 4. (Color online) MC results for C ijkl as a function of distance.(a) Longitudinal correlations C ijkl k for w = (0 , and various sys-tem sizes. (b) C ijkl k along zero-phase anti-nodal directions (filled cir-cles in Fig. 3) for L = 128 and w = (0 , , (1 , , (0 , and (1 , .(c) C ijkl k ( r ) / cos( Q · r ) and C ijkl ⊥ ( r ) / sin( Q · r ) for w = (0 , and L = 128 (data separated from a nodal line by a parallel displacement dx < are excluded). (d) Absolute value of transverse correlations C ijkl ⊥ for w = (0 , and L = 128 , along the line highlighted inFig. 3(a). In (a-c) the line indicates our best fit yielding the exponent α = 1 . and in (d) the fit yielding α ′ = 2 . . unity (see below).In Fig. 4(a-b) we plot C ijkl k ( r ) along the anti-nodal lineswith strongest correlations (hence smallest relative errors),highlighted in Fig. 3, for: (a) w = (0 , and various sys-tem sizes L and (b) L = 128 and w = (0 , , (1 , , (0 , and (1 , . Results are consistent with power-law de-cay, C ijkl k ∼ r − α , as conjectured in Refs. 18 and 31. Fittingthe data in Fig. 4(b) we arrive at α = 1 . . Although smalldeviations from algebraic behavior are seen for large distancesin Fig. 4(a-c), the value of r at which they start to occur in-creases linearly with L (not shown) and we thus conclude thatthis “upturn” in Fig. 4(a-c) is merely a finite-size effect.We analyze the spatial modulations for correlations and inFig. 4(c) we plot C ijkl k ( r ) and C ijkl ⊥ ( r ) for L = 128 and w = (0 , . By respectively dividing C ijkl k ( r ) , C ijkl ⊥ ( r ) by the phase factors cos( Q x dx ) , sin( Q x dx ) [ dx = 0 alongthe anti-nodal line for C ijkl k ( r ) highlighted in Fig. 3(c) and Q = πL (1 , in this case], we notice that all curves collapse,confirming that C ijkl k ( r ) and C ijkl ⊥ ( r ) decay with the sameexponent and are indeed modulated according to the phasefactors cos( Q · r ) and sin( Q · r ) , with Q = πL ( w y , w x ) .Sub-leading corrections to the scaling exponent can be ob-tained by analyzing C ijkl ⊥ in the w = (0 , sector [seeFig. 3(a)]. Data for C ijkl ⊥ are plotted as a function of r in Fig. 4(d), for L = 128 and w = (0 , . A fit yields α ′ = 2 . for the sub-leading exponent.Finally, we address the point of what specific type of quasi-long-range dimer order is encoded in Eq. (1). In doing so,we analyze the dimer order parameter D defined by D α = N − P r ( − r α S r · S r + e α . Due to the absence of long-rangeorder, h D i is expected to vanish when L → ∞ . However,information concerning the symmetry of the quasi-orderedstate is obtainable by analyzing the angular dependence on φ = arctan( D y /D x ) in the histogram P ( D x , D y ) for oc-currences of D x and D y in the simulations. In Fig. 5 weplot P ( D x , D y ) for L = 96 and w = (0 , . Commonlyobserved VBCs on the square lattice, “columnar” and “pla-quette” states (for a detailed account see Ref. 21), would dis-play, respectively, peaks located at φ = { , ± π/ , π } and φ = {± π/ , ± π/ } . However, no angular structure is ev-ident in Fig. 5 and data are in favor of an continuous U (1) symmetry, a priori different from the U (1) symmetry associ-ated to topological degeneracy. Conclusions —
We have investigated NN-VB wave func-tions on the square lattice [Eq. (1)], characterized by wind-ing numbers w , by performing MC simulations. We con-firm earlier findings in favor of short-ranged spin order and, more interestingly, we find that dimer-dimer correla-tions are critical , a situation reminiscent of the one encoun-tered for classical dimers and thus for the ground-state ofthe quantum dimer model (QDM) on the square lattice at theRokhsar-Kivelson point. However, such correlations decayconsiderably slower for Eq. (1) than in the classical case, sug-gesting increased tendency towards VBC order: an exponent α = 1 . accounts for the decay of both longitudinal andtransverse correlation for all w studied, while in the classical L = w = (0, 0)-0.05 0 0.05 D x -0.05 0 0.05 D y ! -5 ! -4 FIG. 5. (Color online) Histogram P ( D x , D y ) for L = 96 and w =(0 , . case one has α class . = 2 . In this context, it would be in-teresting to analyze how exponents evolve by considering theoverlap as a tunable parameter, between the herein studiedcase and the limit of orthogonal dimer configurations of theQDM. From a broader perspective, we analyze how the wave-function Eq. (1) fits into the general classification of QSLstates. While we stress that our work concerns wave-functions and not the full spectrum of a given Hamiltonian,we notice that each | ψ w i [Eq. (1)] is a (degenerate) ground-state of the local model of Ref. 28 on a torus. We thus predictsuch SU (2) -invariant Hamiltonian to display gapped spin ex-citations, due to short-ranged spin order, and gapless non-magnetic excitations , since four-point correlations are criti- cal and a theorem by Hastings applies. Adopting the termi-nology of Ref. 7, the latter excitations correspond to gaplessgauge modes. Given its extensive degeneracy on a torus, ourresults altogether suggest that the ground-state of the modelof Ref. 28 is a gapped U (1) or SU (2) spin liquid, a state be-lieved to be unstable for a generic local spin model. In thiscontext, the complete characterization of the spectrum of themodel of Ref. 28 and of perturbations thereof would be of highinterest.Directions for further research opened up by our work in-clude the study of wave-functions similar to Eq. (1) in dif-ferent geometries. In d = 3 , we expect NN-VB states onbipartite lattices to display spin order, as it happens for thesimple cubic lattice, but it would be interesting to searchfor traces of the Coulomb phase of dimer models in d = 3 . NN-VB wave-functions on d = 2 geometrically frustrated lat-tices also deserve investigation. However, a sign problem precludes MC simulations as we perform here and alternativeapproaches are called for in this case.
Note Added —
While preparing this manuscript we becameaware of related work by Tang et al. . Acknowledgments —
We thank P. Fendley, M. Mambriniand G. Misguich for useful discussions. Our MC codesare based upon the ALPS libraries.
This work was per-formed using HPC resources from GENCI-CCRT (Grant2010-x2010050225) and CALMIP, and is supported by theFrench ANR program ANR-08-JCJC-0056-01. L. Balents, Nature, , 199 (2010). P. W. Anderson, Mater. Res. Bull., , 153 (1973). P. Fazekas and P. W. Anderson, Philos. Mag., , 474 (1974). P. W. Anderson, Science, , 1196 (1987). G. Misguich, in
Exact Methods in Low-dimensional StatisticalPhysics and Quantum Computing , edited by J. Jacobsen et al. (Oxford Univ. Press, Oxford, 2010). G. Misguich and C. Lhuillier, in
Frustrated spin systems , editedby H. T. Diep (World-Scientific, Singapore, 2005). X.-G. Wen, Phys. Rev. B, , 165113 (2002). G. Misguich et al. , Phys. Rev. B, , 1064 (1999). L. Balents, M. P. A. Fisher, and S. M. Girvin, Phys. Rev. B, ,224412 (2002). M. Hermele, M. P. A. Fisher, and L. Balents, Phys. Rev. B, ,064404 (2004). S. Fujimoto, Phys. Rev. B, , 024429 (2005). K. S. Raman, R. Moessner, and S. L. Sondhi, Phys. Rev. B, ,064413 (2005). A. Seidel, Phys. Rev. B, , 165131 (2009). A. Kitaev, Ann. Phys., , 2 (2006). H. Yao and S. A. Kivelson, Phys. Rev. Lett., , 247203 (2007). Z. Y. Meng et al. , Nature, , 847 (2010). S. Liang, B. Douc¸ot, and P. W. Anderson, Phys. Rev. Lett., ,365 (1988). B. Sutherland, Phys. Rev. B, , 3786 (1988). K. S. D. Beach and A. W. Sandvik, Nucl. Phys. B, , 142(2006). J. Lou and A. W. Sandvik, Phys. Rev. B, , 104432 (2007). M. Mambrini et al. , Phys. Rev. B, , 144422 (2006). A. W. Sandvik and H. G. Evertz, Phys. Rev. B, , 024407 (2010). P. W. Kasteleyn, Physica, , 1209 (1961). H. N. V. Temperley and M. E. Fisher, Philos. Mag., , 1061(1961). M. E. Fisher, Phys. Rev., , 1664 (1961). M. E. Fisher and J. Stephenson, Phys. Rev., , 1411 (1963). N. E. Bonesteel, Phys. Rev. B, , 8954 (1989). J. Cano and P. Fendley, Phys. Rev. Lett., , 067205 (2010). K. S. D. Beach, (2007), cond-mat: 0707.0297. Similar effects are seen in correlators for classical dimers withfixed w , yet with a different spatial structure (unpublished). H. Tasaki, Phys. Rev. B, , 9183 (1989). D. S. Rokhsar and S. A. Kivelson, Phys. Rev. Lett., , 2376(1988). M. Kohmoto and Y. Shapir, Phys. Rev. B, , 9439 (1988). M. B. Hastings, Phys. Rev. B, , 104431 (2004). D. A. Huse, W. Krauth, R. Moessner, and S. L. Sondhi, Phys.Rev. Lett., , 167004 (2003). Y. Tang, A. W. Sandvik, and C. L. Henley, (2010),arXiv:1010:6146. Preliminary results were announced in: Y. Tangand A. W. Sandvik, Bull. Am. Phys. Soc. , P38.12 (2010). M. Troyer, B. Ammon, and E. Heeb, Lecture Notes in Comput.Sci., , 191 (1998). A. F. Albuquerque et al. , J. Magn. Magn. Mater.,310