Entanglement Perturbation Theory for Infinite Quasi-1D Quantum Systems
aa r X i v : . [ c ond - m a t . s t r- e l ] O c t APS/123-QED
Entanglement Perturbation Theory for Infinite Quasi-1D Quantum Systems
Lihua Wang ∗ and Sung Gong Chung , Computational Condensed-Matter Physics Laboratory, RIKEN, Wako, Saitama 351-0198, Japan Department of Physics and Nanotechnology Research and Computation Center,Western Michigan University, Kalamazoo, Michigan 49008, USA and Asia Pacific Center for Theoretical Physics, Pohang, Gyeonbuk 790-784, South Korea (Dated: June 18, 2018)We develop Entanglement Perturbation Theory (EPT) for infinite Quasi-1D quantum systems.The spin 1 / ,094430(2010)] that dimer and N´eel orders appear alternately as the XXZ anisotropy ∆ approachesthe isotropic limit ∆ = 1. The first and second transitions (across dimer, N´eel, and dimer phases)are detected with improved accuracy at ∆ ≈ .
722 and 0 . ≈ .
98, is not detected at this ∆ in our method,raising the possibility that the second N´eel phase is absent.
PACS numbers: 75.10.Pq , 75.10.Jm , 75.40.Mg
In a recent article, Furukawa et al. [1] reported exoticN´eel and dimer orders in a spin 1 / H = X δ =1 L X i =1 J δ (cid:0) S xi S xi + δ + S yi S yi + δ + ∆ S zi S zi + δ (cid:1) (1)characterized by a ferromagnetic NN exchange coupling J < J >
0, with an easy-plane anisotropy ∆ <
1. Let us de-scribe the effect of the NNN exchange by µ ≡ J /J .Fig.1 shows an equivalent coupled 2 spin-chains (zigzagchain) suited for EPT. The theoretical tools employedare the bosonization-Sine-Gordon (SG) theory [2], theRG (renormalization group)-level spectroscopy analysisaugmented by exact diagonalization [3, 4], and numeri-cal RG analysis based on iTEBD (infinite time evolvingblock decimation) [5], a variance of DMRG (density ma-trix renormalization group) method [6] for infinite 1Dsystems. The key finding is that the N´eel order anddimer order alternate numerous times when approachingthe ferromagnetic transition point ( J / | J | , ∆) = (1 / , Js are both positive, the ground state phase di-agram of the model is well understood [2, 7, 8]. Particu-larly, for small µ and with easy-plane anisotropy ∆ <
1, itis a gapless Tomonaga-Luttinger liquid (TLL). The gap-less TLL phase has instabilities toward ordered states[7]. N´eel order will be induced by easy-axis anisotropy∆ >
1, whereas dimer order appears when µ is largerthan some critical ∆-related value. For instance, transi- FIG. 1: (color online) A zigzag chain suitable for EPT. Cir-cles refer to spins. Long dashed lines (blue) indicate theNN coupling J , while short ones (orange) the NNN coupling J = µJ . tion takes place when µ ≈ .
24 at ∆ = 1. The quantumphase transition to these ordered phases can be under-stood within a bosonization-SG theory [2].Now in the case J < J = 0 isknown exact, and the ”infinite” alternation of energiesof the N´eel and dimer states was observed in the pa-rameter space 0 < ∆ <
1, and based on a perturbativetreatment, from the Gaussian Hamiltonian, of the cosineterm in the SG theory, the authors of [1, 9] claimed thatthe same should occur deep inside the frustrated regime | µ | ∼ /
4, and thus the realization of the true alternatingappearance of the N´eel and dimer phases. This predic-tion was first checked by the level spectroscopy [3, 4],showing 5 times of phase alternation between the twophases. This level spectroscopy analysis was, however,still based on the effective SG theory whose validity in theparameter regime | µ | ∼ / √ π ). Thus, authors of [1] also carried out an in-dependent iTEBD calculation of the original model (1),confirming only one dimer to N´eel transition at ∆ ∼ . ∼ .
93. It was re-ported that iTEBD experienced a poor convergence for0 . ≤ ∆ < . < | µ | < . ≤ ∆ < ≥ .
93, only the dimer phaseexists. This means two points. (2) First, the perturba-tive SG analysis when applied to the deeply frustrativeregime | µ | ∼ /
4, is missing some physics, and (3) sec-ond, iTEBD, even with a large entanglement χ ∼ χ ∼
80. This conclusionis not a surprise. Indeed, TEBD operates only on a fewlocal spins, and therefore the deterioration of accuracydue to loss of information arising from successive Hilbertspace truncation, typical of real-space RG, is expectedworse than DMRG.The EPT algorithms EPT-g1 and EPT-g2 both solvefor ground state properties of systems with translationalsymmetry (a recent development of EPT handling inho-mogeneity will be reported elsewhere). EPT-g1 is es-pecially suitable for infinite systems and is our startingpoint. It was previously applied to the J = 0 case[14, 15], and its superb accuracy is demonstrated in thefirst confirmation of a prediction made by the confor-mal field theory [16, 17] on a long-range spin-spin cor-relation which can only be seen for thousands of latticeseparation. It also confirmed yet another field theoryprediction, by bosonaization, for the 4-spin correlationfunctions with better precision than DMRG [18].The central idea of EPT-g1 is as follows. First is thesystem wave function in the matrix product state (MPS)[10, 19–23] representation | ψ i = X { σ } Tr [ ξ σ · ξ σ · · · ξ σ N N ] | σ i ⊗ | σ i ⊗ · · · ⊗ | σ N i (2)where ξ σ i i is the local wave function matrix at site i , and N is the number of sites. σ i (= ↑ , ↓ ) refers to local states.The dimension p of the square matrices ξ σ i i , the entan-glement, controls the accuracy of the wave function | ψ i .In principle | ψ i with optimized ξ σ i i becomes exact with p → ∞ . The matrices ξ σ i i are optimized to maximize thevariational energy: ∂∂ξ i h ψ ( ξ , ξ , · · · , ξ N ) | e − βH | ψ ( ξ , ξ , · · · , ξ N ) ih ψ ( ξ , ξ , · · · , ξ N ) | ψ ( ξ , ξ , · · · , ξ N ) i = 0(3)where β → X i ( ξ , · · · , ξ i − , ξ i +1 , · · · , ξ N ) ξ i = ǫY i ( ξ , · · · , ξ i − , ξ i +1 , · · · , ξ N ) ξ i ( i = 1 , · · · , N ) , (4)where X i and Y i are symmetric matrices depending on { ξ i } , and ǫ is the eigenvalue. The eigenfunction { ξ i } withthe largest ǫ describes the ground state wave function | ψ i .The density matrix K ≡ e − βH can also be writtenin matrix product form, the matrix product operator(MPO) [10, 15]. A generalization to quasi-1D naturallyinherits this MPO of density matrix along with MPS forthe wave function. But now, cf. Fig.2, MPS is defined ona 2-spin composite (grey and black vertical two circles)and MPO is in a complex shell by shell (inner product di-rection) and layer by layer structure (vertical direction).We first deform a zigzag chain into two chains shown inFig.1. The Hamiltonian is rewritten as a summation oflocal bond Hamiltonians, H = X bond H bond (5)There are four types of bonds. Bond { , } , { , } etcform the first series. { , } , { , } etc belong to thesecond. These two series account for NN interactions.The third is { , } , { , } , { , } , { , } etc and the forthis { , } , { , } , { , } , { , } etc. They account for NNNinteractions. The way to group is not unique, but itdoes not affect the result. It is immediately seen that (cid:2) H bond i , H bond j (cid:3) = 0 if they are in the same series. Thisis the criteria to group the bonds. Small β safely sepa-rates series in the density matrix as follows e − βH = Y bond e − βH bond + O (cid:0) β (cid:1) ≈ · · · e − β P − st series H bond e − β P − nd series H bond · · · (6)And for each series, we have e − β P series H bond = · · · f α ⊗ g α ⊗ f γ ⊗ g γ · · · + O (cid:0) β (cid:1) (7)where repeating indexes follow Einstein summation con-vention. f , for instance for the XXX model, takes four2 × , √ βS x , √ βS y , √ βS z on a site. And g likewise on the other site of a bond. We use shell to de-note a series in (7). There are four shells in total. Theyare coupled in the direction of in and out of the paper inFig.2, by the inner product of the operators on the samesite. Also note that MPO contains two layers , composedof two vertically aligned spins, one in the top chain andthe other in the bottom chain. We finally have an explicit FIG. 2: Schematic figure of the MPO of the density matrixfor a zigzag chain. Short dashed lines (orange) indicate NNinteractions and long dashed lines (blue) NNN interactions.Solid lines refer to the inner prodoct.FIG. 3: Simplified schematic figure of the MPO of the densitymatrix for a zigzag chain. The two spins enclosed by an ellipseform an effective site. The dashed lines (orange) indicate thedirect product and solid lines the inner product. expression of the whole density matrix as follows K ≡ · · · Γ αβγ,α ′ β ′ γ ′ ⊗ Γ α ′ β ′ γ ′ ,α ′′ β ′′ γ ′′ (8)where Γ αβγ,α ′ β ′ γ ′ = f α ′ ( µJ ) f i ( J ) g β ( J ) g α ( µJ ) ⊗ (9) f γ ′ ( µJ ) f β ′ ( J ) g i ( J ) g γ ( µJ )and Γ likewise. α , β and γ etc run from 1 to 4. Theyare the indexes involved in the direct product. If theyare absorbed into an effective index, both Γ and Γ be-come 4-leg tensors. The two effective legs run from 1 to4 = 64. And the other two legs run from 1 to 4, respon-sible for the inner product with a new local wave function ξ now defined on a 2-spin composite, denoted by an el-lipse in Fig.3. It is also noted that, reflecting the Γ , repetition structure, the local wave functions are also oftwo species ξ , . The generalized eiegnvalue problem (4) is solved iteratively for the local wave functions, and af-ter convergence, energy etc can be calculated as before[11, 15].The calculation on zigzag chains is more expensivecompared with a single spin chain because one now hasto deal with the larger number of local states and largerindexes in MPO. Let us discuss the major time consump-tion. There are two. The first are generalized eigenvalueequations whose size is np with n being the number ofeffective local states. The second is an eigenvalue decom-position of a transfer matrix in the horizontal directionin Fig.3, a basis unit involving Γ and Γ and in total4 2-spin composites connected to Γ , . Its size is mp where m is the size of the 3 bonds in the direct productbetween Γ and Γ . In contrast to n = 2 and m = 4 ina spin chain, they are 4 and 64 now. It is seen that theeigenvalue decomposition of mp x mp transfer matrixbecomes especially large.Fortunately, there is a simple way to overcome thisdifficulty. The key point is in the original idea of EPT-g1, namely convert the Hamiltonian eigenvalue prob-lem to that for the density matrix K ≡ e − βH , noting e − βH → − βH , leading to an MPO representation ofthe density matrix. A simple example is the product oftwo local density matrixes, (1 − βH ij ) ⊗ (1 − βH kl ) ∼ − βH ij − βH kl + β H ij ⊗ H kl , where the last term isunimportant in the β → β = 10 − ∼ − isused in this work like before [11]. Likewise, Γ , αβγ,α ′ β ′ γ ′ should contain such higher order terms in β which wecan consistently discard. A rule to do so is in the factthat, in f α and g α matrixes, α = 1 is identity matrixwhile other three terms are all order √ β . We shouldkeep in mind that whenever there is a term of order √ β ,there should be a companion order √ β term at the otherend of the bond. This amounts to keep, among 64 indexcombinations in αβγ , those terms with only one of αβγ not being 1, namely, 111 , , , , , ,
311 and114 , , d ≡ D ~S i (cid:16) ~S i − − ~S i +1 (cid:17)E for µ = 0 . d is maximum when both Js are positive [8]. TheEPT result is 0 . p = 35 vs. 0 . µ = 2 . p -equivalent χ = 200 ∼ h S zi i and the dimer order parameter D xy is N ee l o r de r pa r a m e t e r D xy || d i m e r o r de r pa r a m e t e r FIG. 4: The phase diagram along the line connecting( | µ | , ∆) = (0 . , .
65) and (0 . , . defined by D xy ≡ (cid:0) S xi − S xi + S yi − S yi (cid:1) − (cid:0) S xi S xi +1 + S yi S yi +1 (cid:1) (10)The main difference between our phase diagram andthat of [1] is that our calculation converged solidly forall ∆ up to 1, whereas iTEBD failed for 0 . < ∆ < > .
93, disproving the prediction of [1] thatthere should be an infinite number of such transitions.The failed prediction, we argue, is most likely due to thefailure of perturbative analysis of the Sine-Gordon theorydeep in the frustrated regime | µ | ∼ /
4. Technically, theLanczos diagonalization not only fails to converge due tohighly degenerate nature for 0 . ≤ ∆ ≤ > .
95 so that the level spectroscopyanalysis concludes incorrectly for this parameter regime.As for the nature of the N´eel to dimer transition, ourcalculation, Fig.4, does not agree with [1] that it is firstorder. Rather it appears second order.While our calculation converged solidly, there is a clearevidence of a hard competition between the N´eel anddimer orders near the transition point. In Fig.5, we haveplotted the N´eel and dimer order parameters as a func-tion of entanglement for ∆ = 0 . .
930 with even moresubtle competition between the two orders since theirmagnitudes are both small.To summarize, we have developed EPT for quasi-1Dstrongly correlated quantum systems. There are two im-mediate applications of interest. First is the 1D Hubbardmodel with long-range hopping, long range Coulomb in-teraction, orbital degeneracy and Hund-coupling in thepursuit of itinerant ferromagnetism and triplet supercon-ductivity. Second is a spin tube, namely the effective
10 20 30 40 50 60 70 80-0.050.000.050.100.150.200.25 p Dxy
FIG. 5: Order parameters vs entanglement p at ∆ = 0 . p . sites now contain the third spin, taking into account thirdnearest neighbor exchange interactions. Works are cur-rently under way, and will be reported elsewhere.Note added – After the completion of the paper, wehave heard from the authors of [1] that the Lifshitz linealong which we and they have done the calculations isactually not precise, and that their prediction of theN´eel phase could still be observed close to the ferro-magnetic point, ∆ = 1 .
0. We have explored ∆ = 0 . J / | J | = 0 . J / | J | = 0 . S x S in this parameter region to be strictly zero. Inthe meantime, we have observed a finite chiral order pa-rameter for the point ∆ = 0 . J / | J | = 0 .
35 whichis deeply in the chiral phase. We thus conclude that ourmain claim in the paper is unchanged.The numerical calculation utilized the RIKEN Inte-grated Cluster of Clusters at Advanced Center for Com-puting and Communication, RIKEN. We thank the au-thors of [1] for useful discussions. ∗ Electronic address: [email protected][1] S. Furukawa, M. Sato, and A. Furusaki, Phys. Rev. B ,094430 (2010).[2] F. D. M. Haldane, Phys. Rev. B , 4925 (1982).[3] K. Okamoto and K. Nomura, Phys. Lett. A , 433(1992).[4] K. Nomura, J. Phys. A , 5451 (1995).[5] G. Vidal, Phys. Rev. Lett. , 070201 (2007). [6] S. R. White, Phys. Rev. B , 10345 (1993).[7] C. K. Majumdar and D. K. Ghosh, J. Math. Phys. ,1399 (1969).[8] S. R. White and I. Affleck, Phys. Rev. B , 9862 (1996).[9] S. Furukawa, M. Sato, and S. Onoda, arXiv:1003.3940v1.[10] S.G.Chung, Phys. Lett. A , 707 (2006).[11] S.G.Chung, Phys. Lett. A , 396 (2007).[12] S.G.Chung and K.Ueda, Phys. Lett. A , 4845 (2008).[13] S.G.Chung and L.Wang, Phys. Lett. A , 2277 (2009).[14] L. Wang, Entanglement perturbation theory for antiferro-magnetic spin chains (ProQuest, UMI Dissertation Pub-lishing, 2011), ISBN 978-1243624109.[15] L.Wang and S.G.Chung, arXiv:1109.6758v1.[16] I. Affeck, J. Phys. A: Math. Gen , 4573 (1998).[17] R. R. P. Singh, M. E. Fisher, and R. Shankar, Phys. Rev.B , 2562 (1989).[18] T.Hikihara and A.Furusaki, Phys. Rev. B , 064427 (2004).[19] S. Rommer and S. ¨Ostlund, Phys. Rev. B. , 2164(1997).[20] D. Perez-Garcia, F. Verstraete, M. M. Wolf, and J. I.Cirac, Quantum Inf. Comput. , 401 (2007).[21] S. R. Clark, J. Prior, M. J. Hartmann, D. Jaksch, andM. B. Plenio, New Journal of Physics , 025005 (2010).[22] J. I. Latorre and V. Picˇo, J. Phys. A: Math. Theor. ,265302 (2009).[23] P. Pippan, S. R. White, and H. G. Evertz, Phys. Rev. B , 081103(R) (2010).[24] Z. L. Berezinskii, Zh. Eksp. Teer. Fiz. , 907 (1970).[25] J. M. Kosterlitz and D. J. Thouless, J. Phys. C , 1181(1973).[26] S. G. Chung, Phys. Rev. B60